diff --git a/.claude/settings.json b/.claude/settings.json new file mode 100644 index 0000000..ed05636 --- /dev/null +++ b/.claude/settings.json @@ -0,0 +1,14 @@ +{ + "hooks": { + "SessionStart": [ + { + "hooks": [ + { + "type": "command", + "command": "bash scripts/ci_setup.sh" + } + ] + } + ] + } +} diff --git a/.dockerignore b/.dockerignore new file mode 100644 index 0000000..deb8270 --- /dev/null +++ b/.dockerignore @@ -0,0 +1,15 @@ +.git +.github +venv +.venv +__pycache__ +*.pyc +.pytest_cache +.ruff_cache +htmlcov +.coverage +data/reference +data/samples +*.tiff +*.tif +docs diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..fc8f5fa --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,48 @@ +name: CI + +on: + push: + branches: [main] + pull_request: + +jobs: + test: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.11", "3.12"] + steps: + - uses: actions/checkout@v4 + + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -e ".[dev,api]" + + - name: Lint (ruff) + run: ruff check src/ tests/ + + - name: Run tests + run: pytest tests/ -q + + docker: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - name: Build image + run: docker build -t dendro:ci . + - name: Smoke-test the API container + run: | + docker run -d -p 8000:8000 --name dendro dendro:ci + for i in $(seq 1 30); do + if curl -fs http://localhost:8000/health >/dev/null; then ok=1; break; fi + sleep 1 + done + curl -fs http://localhost:8000/health + test -n "$ok" + docker rm -f dendro diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..f390021 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,33 @@ +# Dendro dating API container. +# +# Builds an image that serves the cross-dating engine over HTTP. By default it +# ships the small synthetic reference fixtures so the API is functional out of +# the box for a demo; mount or copy a real ITRDB corpus and point +# DENDRO_REFERENCE_DIR at it for production use: +# +# docker run -p 8000:8000 \ +# -e DENDRO_REFERENCE_DIR=/data/reference \ +# -v /path/to/itrdb:/data/reference dendro +FROM python:3.11-slim AS base + +ENV PYTHONUNBUFFERED=1 \ + PIP_NO_CACHE_DIR=1 \ + DENDRO_REFERENCE_DIR=/app/data/reference + +WORKDIR /app + +# Install dependencies first for better layer caching. +COPY pyproject.toml setup.py README.md ./ +COPY src ./src +RUN pip install --upgrade pip && pip install ".[api]" + +# Provide the synthetic fixtures as a default (demo) reference set. +COPY tests/fixtures/reference/synth /app/data/reference/synth + +EXPOSE 8000 + +# Simple container healthcheck against the API. +HEALTHCHECK --interval=30s --timeout=5s --start-period=10s --retries=3 \ + CMD python -c "import urllib.request,sys; sys.exit(0 if urllib.request.urlopen('http://localhost:8000/health').status==200 else 1)" + +CMD ["dendro-api"] diff --git a/README.md b/README.md index 34716f8..2c9f0e1 100644 --- a/README.md +++ b/README.md @@ -1,16 +1,59 @@ # Dendrochronology Dating Tool -A Python CLI tool for dating historic timber by cross-dating ring width measurements against ITRDB reference chronologies. +A Python tool and HTTP service for dating historic timber by cross-dating ring width measurements against ITRDB reference chronologies. Usable from the command line or as a cloud-deployable API. ## Overview This tool helps determine when trees were felled by: 1. Extracting ring width measurements from scanned wood cross-sections 2. Cross-dating against publicly available ITRDB chronologies -3. Identifying the calendar year of the outermost ring (felling year if bark edge present) +3. Interpreting the calendar year of the outermost ring as a felling date (exact with bark edge; a range or *terminus post quem* otherwise) Designed for dating historic New England timber (1600s-1800s), but works for any region with ITRDB coverage. +## Running as a cloud service + +The cross-dating engine is exposed as a stateless FastAPI service (the engine is +used as a library; the API does not shell out to the CLI): + +```bash +pip install -e ".[api]" +export DENDRO_REFERENCE_DIR=/path/to/itrdb # directory of .rwl/.crn files +dendro-api # serves on :8000 +``` + +Endpoints: `GET /health`, `GET /references`, `POST /date` (JSON ring widths), +`POST /date/csv` (CSV upload), `POST /parse` (Tucson file upload). Example: + +```bash +curl -X POST http://localhost:8000/date -H 'content-type: application/json' \ + -d '{"widths": [1.2, 0.9, ...], "has_bark_edge": true, "era_start": 1700, "era_end": 1850}' +``` + +Or with Docker (ships the synthetic demo references; mount a real corpus for production): + +```bash +docker build -t dendro . && docker run -p 8000:8000 dendro +``` + +### Felling-year interpretation + +Dating establishes the year of the **last measured ring**. What that implies for +the felling date depends on the sample: + +- **Bark/waney edge present** → exact felling year. +- **Sapwood present, no bark edge** → estimated felling-year *range* (species sapwood model). +- **No sapwood retained** → *terminus post quem* ("felled after"). + +Pass `--bark-edge/--no-bark-edge`, `--has-sapwood`, and `--sapwood-count` to the +`dendro date` command (or the corresponding API fields). + +### Series orientation + +All series are dated in canonical oldest→newest (pith→bark) order. If your +measurements run bark-inward, pass `--orientation bark_to_pith` and they are +reversed automatically. + ## Installation ```bash @@ -197,7 +240,8 @@ dendrochronology/ │ └── cli/main.py # CLI entry point ├── data/ │ └── reference/ # Downloaded ITRDB chronologies -└── tests/ # Test suite (43 tests) +│ └── api/app.py # FastAPI cloud service +└── tests/ # Test suite (83 tests + fixtures) ``` ## Species Codes @@ -268,7 +312,10 @@ source venv/bin/activate pytest tests/ -v ``` -47 tests total, including 4 real-data validation tests using ITRDB chronologies. +83 tests total. Most run offline against committed synthetic fixtures +(`tests/fixtures`), including an end-to-end dating test that recovers a known +felling year. Four additional validation tests use real ITRDB chronologies and +skip automatically unless that data has been downloaded. ## Algorithm Validation diff --git a/docs/TECHNICAL_PLAN.md b/docs/TECHNICAL_PLAN.md new file mode 100644 index 0000000..5445445 --- /dev/null +++ b/docs/TECHNICAL_PLAN.md @@ -0,0 +1,384 @@ +# Dendro: Technical Review & Plan for a Cloud Dendrochronology Platform + +**Status:** Draft for discussion +**Scope:** Whole codebase (`src/dendro`, `tests`, `scripts`, packaging) +**Goal:** Turn the current local CLI prototype into *actually functional* dendrochronology +software that runs as a service in the cloud. + +--- + +## Implementation status (this branch) + +A first implementation pass has landed the science-correctness fixes, the +testing/CI foundation, and a deployable service layer. What is **done**: + +- **§3.1 Orientation** — canonical oldest→newest enforced end-to-end via + `normalize_orientation`; `--orientation` / API field; reversed input recovers + the same date (tested). +- **§3.2 Felling logic** — new `crossdating/felling.py`; reports distinguish + exact felling year (bark edge), estimated range (sapwood), and *terminus post + quem* (tested). +- **§3.3 RWL master chronology** — now detrend-per-series → biweight mean of + indices, cached. On the synthetic fixture this changed the result from a + *confidently wrong* 1739 (r=0.89) to the correct **1789 (r=0.97, HIGH)**. +- **§3.4 Gleichläufigkeit** — formula corrected, significance test added, folded + into the confidence model (tested). +- **§3.5 Missing rings** — first-cut `detect_missing_ring` locates a likely + missing/false ring and the matcher surfaces it as a warning (tested). Full + auto-realignment search remains future work. +- **§3.6 Spline** — replaced with the Cook & Peters cubic smoothing spline; + frequency response verified to be 0.5 at the cutoff wavelength (tested). +- **§3.7/3.8 Parsers/cleanup** — Tucson writer rounding fixed; a writer↔parser + round-trip test added; era filtering rationalized to the felling-year window. +- **§4.2 Imaging** — directional gradient stops double-counting ring boundaries. +- **Phase 0** — committed deterministic fixtures, CI (ruff + pytest + Docker + smoke test), Dockerfile, lockfile, ruff config, SessionStart hook. 79 tests + pass offline; the 4 real-ITRDB tests skip without downloaded data. +- **Phase 2 (spine)** — FastAPI service (`api/app.py`) wrapping the engine as a + library, with input validation, a cached reference index, and a Docker image. + +Still **outstanding** (tracked below): full missing-ring re-alignment search, +the browser measurement workflow (Phase 3), a persistent project/database model +and job queue (the rest of Phase 2), real-corpus reference data management, and +validation against dplR. + +--- + +## 1. Executive summary + +The repository is a well-organized **prototype** with the right module boundaries +(reference data, imaging, cross-dating, visualization, CLI) and good documentation. It +demonstrably recovers known dates on *clean, same-site* ITRDB series in its synthetic and +within-site validation tests. + +However, it is **not yet scientifically reliable** and **cannot run in the cloud as-is**. +The blockers fall into three groups: + +1. **Correctness bugs in the dating science** that will silently produce wrong or + weakened results on real-world samples (sample orientation is unenforced, the + bark-edge flag is ignored when reporting a felling year, the RWL master chronology is + built by averaging *raw* widths instead of detrended indices, the Gleichläufigkeit + formula is a no-op, and there is no missing-ring handling). +2. **A measurement front-end that requires a local GUI** (matplotlib `plt.show()`), + plus ring detection that is too rudimentary for real hardwood scans. +3. **No service architecture** — no HTTP API, no persistence, no job processing, no + container, no CI, and reference data that is scraped from NOAA HTML at runtime and + re-parsed on every invocation. + +This document inventories the specific issues and proposes a phased plan. The +recommended sequencing is: **fix the science first** (Phase 1), because a cloud service +that returns confidently wrong dates is worse than no service; then **build the service +spine** (Phase 2); then **replace the measurement GUI with a browser workflow** (Phase 3); +then **harden and operationalize** (Phase 4). + +--- + +## 2. What works today (keep this) + +- Clean package layout under `src/dendro` with sensible separation of concerns. +- Tucson `.rwl` reading is good enough for the common ITRDB layout and is the backbone + of the working validation tests. +- The sliding-correlation core (`crossdating/correlator.py`) is straightforward and + correct for the rigid-alignment case; t-value computation is right. +- Diagnostic plotting (`visualization/plots.py`) is genuinely useful and is the kind of + "show your evidence" output dendrochronologists expect (overlay, correlation profile, + segment bars, marker years). +- The validation harness concept (`scripts/validate_crossdating.py`, `tests/test_validation.py`) + — "forget the date, recover it" — is exactly the right way to test this domain. + +--- + +## 3. Correctness issues in the dating science (highest priority) + +These are ordered by how badly they corrupt real-world results. + +### 3.1 Sample orientation is ambiguous and unenforced — *critical* +`CrossdateMatcher.date_sample` documents `values` as *"from bark to pith, outer to inner"* +(newest→oldest), but reference chronologies are stored **oldest→newest**, and +`sliding_correlation` compares the two arrays element-by-element. If the sample is in the +opposite temporal order from the reference, every correlation is computed on a +time-reversed series and dating silently fails. + +Meanwhile the measurement export path (`path_sampler.widths_to_csv` / +`widths_to_tucson`) *reverses* widths to oldest-first — the opposite of what the matcher +docstring asks for. So the tool's own output and its dating engine disagree about +orientation. Depending on how a CSV was produced, results are either correct or garbage, +with no error surfaced. + +**Fix:** Define one canonical internal orientation (recommend **pith→bark / oldest→newest**, +matching ITRDB), convert at every boundary (CSV in, RWL in, measurement out), label the +orientation explicitly in file metadata, and add a guard that runs the correlation both +ways during development to detect regressions. + +### 3.2 The `bark_edge` flag never affects the reported felling year — *critical* +`MatchResult.felling_year` simply returns `proposed_end_year`. The `has_bark_edge` +boolean is threaded through the entire stack but is **never used** in the computation. +For a sample *without* bark/waney edge, the last measured ring is **not** the felling +year — there are unknown missing sapwood (and possibly heartwood) rings. Reporting +`proposed_end_year` as "PROPOSED FELLING YEAR" in that case is scientifically wrong and is +the single most consequential output of the whole tool. + +**Fix:** +- With bark edge → felling year = last ring year (state precision: exact, ± nothing). +- Without bark edge but with sapwood → "felling after" (*terminus post quem*) plus a + sapwood-estimate range (e.g. Hollstein / Baillie & Pilcher sapwood models for oak; + species-specific). Report a range, not a point. +- Heartwood-only → report only "felling after last heartwood ring + minimum sapwood." +- Make the report type distinguish *last-ring date*, *felling year*, and *felling-year + range* as separate fields rather than overloading one number. + +### 3.3 RWL master chronologies are built from raw widths — *high* +In `matcher._match_against_reference`, the RWL branch does: +```python +ref_values = df.mean(axis=1).values # mean of RAW ring widths across cores +detrended, _ = detrend_series(ref_values) +``` +Averaging *raw* widths across cores of different ages/growth rates, then detrending the +average, is methodologically wrong. The standard is to **detrend each series to a +dimensionless index first, then average the indices** (optionally with a robust/biweight +mean). The repo already does it correctly in the test helper +`tests/test_validation.py::_build_site_master` and has `detrend.build_chronology` — but the +production matcher doesn't use either. Every RWL-based match is weakened by this. + +**Fix:** Route all RWL references through a single, correct chronology builder +(detrend-per-series → biweight robust mean → standardized master), cache the result. + +### 3.4 Gleichläufigkeit is computed incorrectly and then ignored — *medium* +In `calculate_gleichlauf`: +```python +glk = 100 * (agreements - zeros * 0.5 + zeros * 0.5) / n +``` +`- zeros*0.5 + zeros*0.5` cancels to zero, so the function returns `100 * agreements / n`, +which is not GLK (it also double-counts zero/zero pairs as agreements). The intended +formula is `GLK = (agreements + 0.5 * semi_agreements) / (n-1)`. Worse, the value is never +used in confidence scoring — GLK is a standard, robust cross-dating statistic and should +be. There's also no GLK significance test (Gsl/`G_sl`). + +**Fix:** Implement GLK correctly, add its significance test, and incorporate it (alongside +t and overlap) into the confidence model. + +### 3.5 No missing/false-ring handling — *high (defines "real" vs "toy")* +The core difficulty of dating field samples is **locally absent rings** (a year with no +measurable growth) and **false/double rings**. These shift everything inward of the error +and destroy a single rigid alignment. The current engine performs exactly one rigid slide +and cannot detect, locate, or compensate for a missing ring. Segment correlation +(`segment_correlation`) can *reveal* that something is wrong but offers no remedy. + +**Fix (phased):** Start with COFECHA-style segmented re-correlation that tests segment +offsets of ±1, ±2 years to flag the likely location of a missing/extra ring; later add a +search that evaluates candidate insertions/deletions and re-scores. This is a substantial +algorithmic addition and should be its own milestone. + +### 3.6 Spline detrending is an ad-hoc approximation — *medium* +`_fit_cubic_spline` maps a target wavelength to scipy's `UnivariateSpline` smoothing +factor via `s = n/(period/2)` then `s*var(y)`. This is not the standard dendro cubic +smoothing spline (Cook & Peters 1981), which is defined by a 50%-frequency-response +cutoff at a given wavelength. Results won't match dplR/standard tooling and aren't +reproducible across series of different length. + +**Fix:** Implement the Cook & Peters cubic smoothing spline (or adopt a vetted +implementation), parameterized by cutoff wavelength as a fraction of series length, and +make the detrending method/curve configurable and recorded in the report. + +### 3.7 Tucson parser robustness — *medium* +- `999` is treated as a stop marker, but `999` in 0.01 mm units is a legitimate 9.99 mm + ring; conflating them can truncate or corrupt wide-ringed series. Tucson's true + terminator is position/format-dependent (and `-9999`/`9990` are the missing/end + sentinels in several variants). +- The CRN reader assumes a rigid 7-char value+depth layout and a hardcoded species + regex; ITRDB has multiple CRN and "-noaa" tabular variants it won't parse. +- Header metadata extraction (lat/long/species/state) is brittle regex over the first 20 + lines. + +**Fix:** Replace with a rigorous, format-aware parser (evaluate adopting an existing +maintained RWL/CRN reader rather than hand-rolling), add a corpus of real-file fixtures, +and parse the structured ITRDB metadata/NOAA template where available instead of regex +scraping. + +### 3.8 Smaller correctness/typing items +- Era filtering in `_match_against_reference` (`era_start - len(sample)`, `+ 50` slop) is + arbitrary and can both admit and reject valid dates inconsistently — replace with an + explicit, documented window on the **felling year**. +- `cross_correlation` (FFT path) is defined but unused and its normalization is suspect; + either fix and use it for speed or remove it. +- `prewhiten` is implemented but never used by the pipeline, and its AR estimation is + questionable; decide whether residual chronologies are in scope and wire it in properly + if so. +- Packaging says Python ≥3.9 but `setup_env.sh` checks for 3.10 and some modules rely on + PEP 604 (`X | Y`) annotations; pin one floor (recommend 3.11) and make it consistent. + +--- + +## 4. Imaging / measurement issues + +### 4.1 The interactive viewer cannot run headless — *critical for cloud* +`MeasurementViewer.show()` calls `plt.show()` and depends on a GUI backend and mouse +events. This is the **single biggest cloud blocker** for the measurement feature: it +fundamentally cannot run on a server. Measurement must move to a browser-based workflow +(image served to a canvas; clicks/paths posted back; server samples the profile and +detects rings) or a headless batch detector with a human review step. + +### 4.2 Ring detection is too naive for real samples — *high* +`detect_rings` uses CLAHE + smoothed gradient magnitude + a hand-rolled local-max scan. +Because it uses `abs(gradient)`, it fires on **both** the earlywood→latewood and the +latewood→earlywood transitions, tending to double-count boundaries. It has no model for +ring curvature, ring-porous hardwood vessels (oak earlywood pores devastate gradient +methods), knots, rays, or reflectance gradients. It will not produce trustworthy widths on +the oak/hemlock the README targets. + +**Fix:** Treat automated detection as *assistive*, always human-reviewed. Improve the +signal model (directional latewood-density gradient, per-species tuning, perpendicular +band integration along the true path normal rather than vertical `axvline`s), and validate +against hand-measured fixtures. Consider established techniques/tools for ring detection as +a reference baseline. + +### 4.3 Minor +- `_update_ring_display` draws vertical lines at the path x-coordinate; for diagonal or + curved paths these mislocate boundaries. +- DPI is user-supplied and trusted; `estimate_dpi` exists but isn't wired into `measure`. + +--- + +## 5. Architecture & cloud-readiness gaps + +### 5.1 No service layer +Today it's a CLI plus an interactive GUI. To "run in the cloud" we need: +- **HTTP API** (recommend **FastAPI**): endpoints for upload, measurement, dating jobs, + results, reference catalog. +- **Async job processing** for cross-dating runs (a worker queue — e.g. RQ/Celery/Arq — + since matching against many references is not instantaneous and shouldn't block requests). +- **Object storage** for uploaded scans and generated plots (S3-compatible). +- **A browser UI** for measurement and for reviewing dating evidence. +- Keep the existing engine as a library the API calls — do **not** shell out to the CLI. + +### 5.2 Reference data is scraped and re-parsed every run +`downloader.py` regex-scrapes NOAA's HTML directory listing at runtime; +`ChronologyIndex` re-parses **every** file on every `CrossdateMatcher` construction, and +`load_chronology` re-parses again during matching. `save_index`/`load_index` exist but are +unused. In a service this is slow, fragile (NOAA layout/rate limits/blocked egress), and +wasteful. + +**Fix:** Curate a **versioned reference dataset** as a build/deploy artifact (or a managed +object-store mirror), parse it **once** into a persisted, serialized index + cached +per-site master chronologies, and load that at service start. Make reference-set version +part of every result's provenance. + +### 5.3 No persistence / project model +A real platform needs durable domain objects: **Project → Structure/Building → Timber/Beam +→ Sample → Measurement series → Dating run → Result**, with provenance (who, when, scanner, +DPI, operator notes) and an audit trail. Today everything is loose CSV/JSON on disk. +**Fix:** Introduce a relational store (Postgres) and a schema for the above, with results +linked to the exact reference-set version and algorithm parameters used. + +### 5.4 No container, CI, lockfile, or reproducible build +No Dockerfile, no CI, `uv.lock` is gitignored, deps unpinned. The project couldn't be +installed reproducibly in a fresh environment during this review. **Fix:** pin deps with a +committed lockfile, add a Dockerfile, add CI (lint + type-check + tests on every PR), and +add a SessionStart hook so web sessions can run tests/linters. + +### 5.5 Input safety +Uploaded CSVs go straight into `pd.read_csv`/`np.loadtxt` with no size/shape validation; +downloaded filenames come from remote HTML with only a `startswith("..")` check (weak +path-traversal guard). **Fix:** validate and bound all inputs; sanitize any +externally-derived filenames; cap upload sizes; never trust scraped paths. + +--- + +## 6. Testing & validation gaps + +- Most unit tests correlate a sample that is an **exact subset** of the reference, so + `r ≈ 1.0` trivially — they exercise plumbing, not dating difficulty. +- The realistic tests (`test_validation.py`) **skip** unless reference data is downloaded, + and no reference data is committed, so in CI the only meaningful tests don't run. +- No fixtures for: missing-ring samples, no-bark-edge samples, hardwood scans, malformed + Tucson files. + +**Fix:** Commit a **small, license-clear fixture set** (a few real RWL files + a couple of +hand-measured sample CSVs with known answers, including at least one missing-ring and one +no-bark-edge case). Make the recover-the-date validation run in CI against those fixtures. +Add property tests for the parsers. + +--- + +## 7. Proposed phased plan + +Each phase is independently shippable and leaves the project in a better state. + +### Phase 0 — Foundations (1 short iteration) +- Commit a lockfile; pin deps; settle the Python floor (3.11). +- Add CI (ruff + mypy + pytest) and a Dockerfile. +- Commit a small real-data fixture set; wire the recover-the-date validation into CI. +- Add a SessionStart hook for web sessions. +**Exit:** `pytest` green in CI on a clean checkout, including ≥3 real-data dating +validations. + +### Phase 1 — Make the science correct (the core of "actually functional") +- **1a** Enforce a single canonical series orientation end-to-end; add guards (§3.1). +- **1b** Make `bark_edge` real: distinct last-ring / felling-year / felling-range outputs, + with species sapwood models (§3.2). +- **1c** Correct RWL master-chronology construction (detrend-then-average, cached) (§3.3). +- **1d** Fix Gleichläufigkeit + significance; fold GLK into confidence (§3.4). +- **1e** Replace ad-hoc spline with Cook & Peters spline; record detrending in reports (§3.6). +- **1f** Harden the Tucson/CRN parsers against real-file variants + fixtures (§3.7). +- **1g** Rationalize era filtering; remove/repair dead code paths (§3.8). +**Exit:** Documented, reproducible dating on the fixture corpus; outputs distinguish +felling-year certainty levels; results match a reference implementation (e.g. dplR) within +tolerance on shared inputs. + +### Phase 2 — Service spine +- FastAPI app wrapping the engine as a library; async job queue; Postgres project model; + object storage; serialized reference index + cached masters loaded at startup; + versioned reference dataset as a deploy artifact (§5.1–5.3). +- Input validation and upload limits (§5.5). +**Exit:** Upload a measurement CSV via API → dating job runs → results + diagnostic plot +retrievable; everything persisted with provenance and reference-set version. + +### Phase 3 — Browser measurement workflow +- Replace the matplotlib GUI with a web canvas: serve the scan, capture path + ring + marks, server-side profile sampling and assistive detection, human review/adjust, export + to the canonical series format (§4.1–4.3). +- Improve assistive ring detection and validate against hand-measured fixtures. +**Exit:** A user can upload a scan, measure rings in the browser, and date the result +without touching a CLI or local GUI. + +### Phase 4 — Algorithmic depth & hardening +- Missing/false-ring detection and assisted re-alignment (§3.5). +- Multi-sample replication / site master building in-app; cross-verification UI. +- Calibrated, species-aware confidence; configurable thresholds with documented provenance. +- Observability (structured logs, metrics), authn/authz, rate limiting, backups. +**Exit:** Handles real field samples with missing rings; defensible confidence reporting; +production-operable. + +--- + +## 8. Build vs. reuse + +Dendrochronology has a de-facto standard implementation in **dplR (R)**. We should not +blindly reimplement decades of refined methods. Recommended posture: +- **Reuse the *algorithms* faithfully** (Cook & Peters spline, biweight robust mean, + standard cross-dating statistics) and **validate our outputs against dplR** on shared + inputs as an acceptance gate. +- **Evaluate adopting a maintained Tucson/RWL parser** rather than hand-rolling format + handling, given how many ITRDB variants exist. +- Keep our own thin engine for the parts we need to control (orientation, sapwood/felling + logic, the service API), so we own the product surface while standing on validated math. + +A short spike at the start of Phase 1 should confirm which external pieces we adopt vs. +reimplement; this plan assumes we reimplement the spline and chronology math (small, +well-specified) and lean on existing parsing/validation references where practical. + +--- + +## 9. Risks & open questions + +- **Reference coverage** for the target region/species/era ultimately bounds dating + success more than any code change; curating the right reference set is a project in its + own right. +- **Missing-ring handling (§3.5)** is the hardest algorithmic work and the difference + between "demo" and "usable on field timber"; it deserves its own design doc. +- **Legal/licensing** of any committed reference fixtures and of derived chronologies + needs checking before redistribution. +- **Scope of measurement-in-browser** (full annotation UI vs. assisted batch) is a product + decision that changes Phase 3 size substantially. + + diff --git a/pyproject.toml b/pyproject.toml index b6e0d5e..7a03378 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,10 +23,19 @@ dependencies = [ dev = [ "pytest>=7.0", "pytest-cov>=4.0", + "httpx>=0.24", + "ruff>=0.4", +] +api = [ + "fastapi>=0.110", + "uvicorn[standard]>=0.29", + "python-multipart>=0.0.9", + "pydantic>=2.0", ] [project.scripts] dendro = "dendro.cli.main:cli" +dendro-api = "dendro.api.app:run" [build-system] requires = ["setuptools>=61.0", "wheel"] @@ -38,3 +47,13 @@ where = ["src"] [tool.pytest.ini_options] testpaths = ["tests"] pythonpath = ["src"] + +[tool.ruff] +line-length = 100 +target-version = "py311" + +[tool.ruff.lint] +select = ["E", "F", "I", "W"] +# E501 (line length) is advisory for this scientific code; bare-except style in +# the legacy modules is tolerated rather than churned in this pass. +ignore = ["E501", "E722", "E741"] diff --git a/requirements.lock b/requirements.lock new file mode 100644 index 0000000..4214442 --- /dev/null +++ b/requirements.lock @@ -0,0 +1,48 @@ +annotated-doc==0.0.4 +annotated-types==0.7.0 +anyio==4.13.0 +certifi==2026.5.20 +charset-normalizer==3.4.7 +click==8.4.1 +contourpy==1.3.3 +coverage==7.14.1 +cycler==0.12.1 +fastapi==0.136.3 +fonttools==4.63.0 +h11==0.16.0 +httpcore==1.0.9 +httptools==0.8.0 +httpx==0.28.1 +idna==3.18 +iniconfig==2.3.0 +kiwisolver==1.5.0 +matplotlib==3.10.9 +numpy==2.4.6 +opencv-python==4.13.0.92 +packaging==26.2 +pandas==3.0.3 +pillow==12.2.0 +pluggy==1.6.0 +pydantic==2.13.4 +pydantic_core==2.46.4 +Pygments==2.20.0 +pyparsing==3.3.2 +pytest==9.0.3 +pytest-cov==7.1.0 +python-dateutil==2.9.0.post0 +python-dotenv==1.2.2 +python-multipart==0.0.32 +PyYAML==6.0.3 +requests==2.34.2 +ruff==0.15.17 +scipy==1.17.1 +six==1.17.0 +starlette==1.3.0 +tqdm==4.68.2 +typing-inspection==0.4.2 +typing_extensions==4.15.0 +urllib3==2.7.0 +uvicorn==0.49.0 +uvloop==0.22.1 +watchfiles==1.2.0 +websockets==16.0 diff --git a/scripts/ci_setup.sh b/scripts/ci_setup.sh new file mode 100755 index 0000000..703f318 --- /dev/null +++ b/scripts/ci_setup.sh @@ -0,0 +1,24 @@ +#!/usr/bin/env bash +# Idempotent environment setup for CI and Claude Code web sessions. +# Ensures the package and its dev/api dependencies are importable so that +# `pytest` and `ruff` work. Safe to run repeatedly; never fails the session. +set -u + +cd "$(dirname "$0")/.." || exit 0 + +# Prefer an existing virtualenv; otherwise install into the active interpreter. +if [ -d venv ]; then + # shellcheck disable=SC1091 + . venv/bin/activate 2>/dev/null || true +fi + +if ! python -c "import dendro" >/dev/null 2>&1; then + python -m pip install --upgrade pip >/dev/null 2>&1 || true + pip install -e ".[dev,api]" >/dev/null 2>&1 || true +fi + +python -c "import dendro" >/dev/null 2>&1 \ + && echo "dendro environment ready" \ + || echo "warning: dendro not importable; run 'pip install -e .[dev,api]'" + +exit 0 diff --git a/src/dendro/api/__init__.py b/src/dendro/api/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/dendro/api/app.py b/src/dendro/api/app.py new file mode 100644 index 0000000..8979074 --- /dev/null +++ b/src/dendro/api/app.py @@ -0,0 +1,249 @@ +""" +HTTP API for the dendrochronology dating engine. + +This is the cloud-facing layer. It wraps the existing cross-dating engine as a +library (it does **not** shell out to the CLI) and exposes stateless endpoints +suitable for running behind a load balancer: + + GET /health liveness/readiness probe + GET /references summary of the loaded reference chronologies + POST /date date a ring-width series supplied as JSON + POST /date/csv date a ring-width series uploaded as a CSV file + POST /parse summarize an uploaded Tucson .rwl file + +The reference index is parsed once at startup and cached in memory (its master +chronologies are also cached lazily by the matcher), so requests do not re-read +the reference corpus. The reference directory is taken from the +``DENDRO_REFERENCE_DIR`` environment variable, defaulting to ``data/reference``. + +Run locally with:: + + dendro-api # or: uvicorn dendro.api.app:app --reload +""" + +from __future__ import annotations + +import io +import os +from contextlib import asynccontextmanager +from pathlib import Path +from typing import Optional + +import numpy as np +import pandas as pd +from fastapi import FastAPI, File, HTTPException, UploadFile +from pydantic import BaseModel, Field + +from ..crossdating.matcher import CrossdateMatcher +from ..reference.chronology_index import ChronologyIndex + +# Guardrails for untrusted input. +MAX_RINGS = 5000 +MAX_UPLOAD_BYTES = 5 * 1024 * 1024 # 5 MB + + +def _reference_dir() -> Path: + return Path(os.environ.get("DENDRO_REFERENCE_DIR", "data/reference")) + + +class _State: + """Process-wide, lazily-loaded reference state.""" + + matcher: Optional[CrossdateMatcher] = None + reference_dir: Optional[Path] = None + + +state = _State() + + +def _load_matcher() -> CrossdateMatcher: + ref_dir = _reference_dir() + state.reference_dir = ref_dir + index = ChronologyIndex(ref_dir) if ref_dir.exists() else ChronologyIndex() + state.matcher = CrossdateMatcher(index=index) + return state.matcher + + +@asynccontextmanager +async def lifespan(app: FastAPI): + _load_matcher() + yield + + +app = FastAPI( + title="Dendro Dating API", + version="0.2.0", + summary="Cross-date tree-ring width series against reference chronologies.", + lifespan=lifespan, +) + + +class DateRequest(BaseModel): + widths: list[float] = Field(..., min_length=1, description="Ring widths.") + sample_name: str = "sample" + has_bark_edge: bool = True + orientation: str = Field("pith_to_bark", pattern="^(pith_to_bark|bark_to_pith)$") + has_sapwood: bool = False + sapwood_count: Optional[int] = None + era_start: int = 1600 + era_end: int = 1900 + species: Optional[list[str]] = None + states: Optional[list[str]] = None + min_overlap: int = Field(30, ge=10, le=500) + top: int = Field(10, ge=1, le=100) + + +def _require_matcher() -> CrossdateMatcher: + if state.matcher is None: + _load_matcher() + assert state.matcher is not None + if len(state.matcher.index) == 0: + raise HTTPException( + status_code=503, + detail=( + "No reference chronologies are loaded. Set DENDRO_REFERENCE_DIR " + "to a directory of ITRDB .rwl/.crn files." + ), + ) + return state.matcher + + +def _validate_widths(widths: list[float]) -> np.ndarray: + if len(widths) > MAX_RINGS: + raise HTTPException(status_code=413, detail=f"Too many rings (> {MAX_RINGS}).") + arr = np.asarray(widths, dtype=np.float64) + if not np.all(np.isfinite(arr)): + raise HTTPException(status_code=422, detail="Widths must all be finite numbers.") + if np.any(arr < 0): + raise HTTPException(status_code=422, detail="Widths must be non-negative.") + return arr + + +def _run_dating(matcher: CrossdateMatcher, req: DateRequest, widths: np.ndarray) -> dict: + report = matcher.date_sample( + values=widths, + sample_name=req.sample_name, + has_bark_edge=req.has_bark_edge, + species_filter=[s.upper() for s in req.species] if req.species else None, + state_filter=[s.upper() for s in req.states] if req.states else None, + era_start=req.era_start, + era_end=req.era_end, + min_overlap=req.min_overlap, + orientation=req.orientation, + has_sapwood=req.has_sapwood, + sapwood_count=req.sapwood_count, + ) + payload = report.to_dict() + payload["matches"] = payload["matches"][: req.top] + return payload + + +@app.get("/health") +def health() -> dict: + n = len(state.matcher.index) if state.matcher else 0 + return { + "status": "ok", + "reference_dir": str(state.reference_dir) if state.reference_dir else None, + "references_loaded": n, + } + + +@app.get("/references") +def references() -> dict: + matcher = state.matcher or _load_matcher() + index = matcher.index + starts = [e.start_year for e in index.entries if e.start_year > 0] + ends = [e.end_year for e in index.entries if e.end_year > 0] + return { + "count": len(index), + "species": index.get_species(), + "states": index.get_states(), + "coverage": { + "start": min(starts) if starts else None, + "end": max(ends) if ends else None, + }, + } + + +@app.post("/date") +def date(req: DateRequest) -> dict: + matcher = _require_matcher() + widths = _validate_widths(req.widths) + return _run_dating(matcher, req, widths) + + +@app.post("/date/csv") +async def date_csv( + file: UploadFile = File(...), + has_bark_edge: bool = True, + orientation: str = "pith_to_bark", + era_start: int = 1600, + era_end: int = 1900, +) -> dict: + matcher = _require_matcher() + raw = await file.read() + if len(raw) > MAX_UPLOAD_BYTES: + raise HTTPException(status_code=413, detail="Uploaded file too large.") + try: + df = pd.read_csv(io.BytesIO(raw)) + except Exception as exc: # noqa: BLE001 - surface parse errors to the client + raise HTTPException(status_code=422, detail=f"Could not parse CSV: {exc}") + + df.columns = [str(c).lower().strip() for c in df.columns] + if "width_mm" in df.columns: + widths = df["width_mm"].to_numpy(dtype=float) + elif "width" in df.columns: + widths = df["width"].to_numpy(dtype=float) + else: + numeric = df.select_dtypes("number") + if numeric.shape[1] == 0: + raise HTTPException(status_code=422, detail="No numeric width column found.") + widths = numeric.iloc[:, -1].to_numpy(dtype=float) + + req = DateRequest( + widths=widths.tolist(), + sample_name=Path(file.filename or "sample").stem, + has_bark_edge=has_bark_edge, + orientation=orientation, + era_start=era_start, + era_end=era_end, + ) + return _run_dating(matcher, req, _validate_widths(req.widths)) + + +@app.post("/parse") +async def parse(file: UploadFile = File(...)) -> dict: + from ..reference.tucson_parser import parse_rwl_file + + raw = await file.read() + if len(raw) > MAX_UPLOAD_BYTES: + raise HTTPException(status_code=413, detail="Uploaded file too large.") + tmp = Path("/tmp") / f"upload_{os.getpid()}_{file.filename or 'series.rwl'}" + tmp.write_bytes(raw) + try: + rwl = parse_rwl_file(tmp) + finally: + tmp.unlink(missing_ok=True) + return { + "filename": file.filename, + "series_count": len(rwl.series), + "series": [ + {"id": sid, "start": s.start_year, "end": s.end_year, "length": s.length} + for sid, s in list(rwl.series.items())[:50] + ], + } + + +def run() -> None: + """Console-script entry point: launch uvicorn.""" + import uvicorn + + uvicorn.run( + "dendro.api.app:app", + host=os.environ.get("DENDRO_HOST", "0.0.0.0"), + port=int(os.environ.get("DENDRO_PORT", "8000")), + ) + + +if __name__ == "__main__": + run() diff --git a/src/dendro/cli/main.py b/src/dendro/cli/main.py index b7ec7e5..06ba56e 100644 --- a/src/dendro/cli/main.py +++ b/src/dendro/cli/main.py @@ -16,7 +16,6 @@ import click import numpy as np - # Default data directory DEFAULT_DATA_DIR = Path.cwd() / "data" @@ -72,7 +71,7 @@ def download(states: str, species: str, output: Optional[str], overwrite: bool): state_list = [s.strip().lower() for s in states.split(",")] species_list = [s.strip().upper() for s in species.split(",")] - click.echo(f"Downloading chronologies for:") + click.echo("Downloading chronologies for:") click.echo(f" States: {', '.join(s.upper() for s in state_list)}") click.echo(f" Species: {', '.join(species_list)}") click.echo(f" Output: {output_dir}") @@ -120,8 +119,8 @@ def measure(image: str, dpi: int, output: Optional[str], auto: bool): Example: dendro measure sample.tiff --dpi=2400 --output=sample.csv """ - from ..imaging.viewer import MeasurementViewer from ..imaging.path_sampler import widths_to_csv + from ..imaging.viewer import MeasurementViewer image_path = Path(image) @@ -151,7 +150,7 @@ def on_complete(w): viewer.show() if widths is not None and len(widths) > 0: - csv_output = widths_to_csv(widths, output_path=output_path) + widths_to_csv(widths, output_path=output_path) click.echo(f"\nSaved {len(widths)} ring measurements to {output_path}") else: click.echo("No measurements recorded.") @@ -200,6 +199,23 @@ def on_complete(w): default=True, help="Sample includes bark edge (for exact felling year)." ) +@click.option( + "--orientation", + type=click.Choice(["pith_to_bark", "bark_to_pith"]), + default="pith_to_bark", + help="Order of the ring series (oldest-first vs bark-first)." +) +@click.option( + "--has-sapwood/--no-sapwood", + default=False, + help="Incomplete sapwood present (no bark edge) for a felling-year range." +) +@click.option( + "--sapwood-count", + type=int, + default=None, + help="Number of sapwood rings present, if known." +) @click.option( "--output", "-o", type=click.Path(), @@ -242,6 +258,9 @@ def date( species: Optional[str], states: Optional[str], bark_edge: bool, + orientation: str, + has_sapwood: bool, + sapwood_count: Optional[int], output: Optional[str], top: int, plot: bool, @@ -263,9 +282,8 @@ def date( dendro date beam1.csv beam2.csv beam3.csv --cross-verify """ import pandas as pd + from ..crossdating.matcher import CrossdateMatcher - from ..crossdating.detrend import DetrendMethod, detrend_series, standardize - from ..crossdating.correlator import sliding_correlation reference_dir = Path(reference) if reference else DEFAULT_DATA_DIR / "reference" @@ -338,6 +356,9 @@ def date( state_filter=state_filter, era_start=era_start, era_end=era_end, + orientation=orientation, + has_sapwood=has_sapwood, + sapwood_count=sapwood_count, ) # Display results @@ -350,7 +371,9 @@ def date( click.echo() if report.consensus_year: - click.echo(f"PROPOSED FELLING YEAR: {report.consensus_year}") + click.echo(f"Last ring dated to: {report.consensus_year}") + if report.felling_estimate: + click.echo(report.felling_estimate.summary()) click.echo(f"Confidence: {report.consensus_confidence}") else: click.echo("No confident date could be determined.") @@ -487,7 +510,7 @@ def parse(rwl_file: str): Example: dendro parse data/reference/nh/nh001.rwl """ - from ..reference.tucson_parser import parse_rwl_file, parse_crn_file + from ..reference.tucson_parser import parse_crn_file, parse_rwl_file filepath = Path(rwl_file) @@ -556,7 +579,7 @@ def _display_segment_analysis(match, consensus_year: int): def _display_marker_years(values: np.ndarray, proposed_start_year: int): """Detect and display potential marker year matches.""" - from ..visualization.plots import detect_marker_years, identify_known_markers, MARKER_YEARS + from ..visualization.plots import MARKER_YEARS, detect_marker_years, identify_known_markers click.echo("\nMarker Year Analysis:") click.echo("-" * 60) @@ -660,12 +683,12 @@ def _run_cross_verify( if len(unique_years) == 1: consensus_year = years[0] click.echo(f"✓ STRONG AGREEMENT: All {len(dated_results)} samples date to {consensus_year}") - click.echo(f" This significantly increases confidence in the dating.") + click.echo(" This significantly increases confidence in the dating.") # Calculate combined statistics total_t = sum(r.matches[0].t_value if r.matches else 0 for _, _, r in dated_results) avg_t = total_t / len(dated_results) - click.echo(f"\n Combined statistics:") + click.echo("\n Combined statistics:") click.echo(f" Samples agreeing: {len(dated_results)}/{len(results)}") click.echo(f" Average t-value: {avg_t:.1f}") @@ -722,8 +745,8 @@ def _run_cross_verify( def _generate_plots(report, values: np.ndarray, matcher, output_path: Path): """Generate and save diagnostic plots.""" - from ..crossdating.detrend import detrend_series, standardize from ..crossdating.correlator import sliding_correlation + from ..crossdating.detrend import detrend_series, standardize from ..visualization.plots import save_diagnostic_plots if not report.matches: diff --git a/src/dendro/crossdating/__init__.py b/src/dendro/crossdating/__init__.py index 9aedf90..a1350d9 100644 --- a/src/dendro/crossdating/__init__.py +++ b/src/dendro/crossdating/__init__.py @@ -1,7 +1,7 @@ """Cross-dating algorithms for matching samples to reference chronologies.""" +from .correlator import calculate_tvalue, sliding_correlation from .detrend import detrend_series, standardize -from .correlator import sliding_correlation, calculate_tvalue from .matcher import CrossdateMatcher, MatchResult __all__ = [ diff --git a/src/dendro/crossdating/correlator.py b/src/dendro/crossdating/correlator.py index d8c9692..b161363 100644 --- a/src/dendro/crossdating/correlator.py +++ b/src/dendro/crossdating/correlator.py @@ -155,25 +155,41 @@ def calculate_gleichlauf(series1: np.ndarray, series2: np.ndarray) -> float: if len(series1) != len(series2) or len(series1) < 2: return 0.0 - # Calculate year-to-year changes - diff1 = np.diff(series1) - diff2 = np.diff(series2) + # Year-to-year interval signs (Eckstein & Bauch 1969). + sign1 = np.sign(np.diff(series1)) + sign2 = np.sign(np.diff(series2)) - # Count agreements - # Sign: positive = +1, zero = 0, negative = -1 - sign1 = np.sign(diff1) - sign2 = np.sign(diff2) + n = len(sign1) - n = len(diff1) - agreements = np.sum(sign1 == sign2) + # Full agreement when both intervals share a sign (including both flat); + # half credit when exactly one interval is flat; no credit on opposite signs. + same = sign1 == sign2 + one_zero = (~same) & ((sign1 == 0) | (sign2 == 0)) + score = np.sum(same) + 0.5 * np.sum(one_zero) - # Count cases where one or both are zero (half weight) - zeros = np.sum((sign1 == 0) | (sign2 == 0)) + return 100.0 * score / n - # Adjusted GLK - glk = 100 * (agreements - zeros * 0.5 + zeros * 0.5) / n - return glk +def gleichlauf_significance(glk_percent: float, n: int) -> float: + """ + Approximate two-sided significance (p-value) of a Gleichläufigkeit value. + + Under the null hypothesis of independent series, GLK is ~50% with standard + error ``1 / (2 * sqrt(n-1))`` (Eckstein & Bauch 1969). This returns the + probability of observing a deviation from 50% at least this large by chance. + + Args: + glk_percent: Gleichläufigkeit as a percentage (0-100). + n: Number of overlapping years (intervals used = n - 1). + + Returns: + Two-sided p-value (lower = more significant). + """ + if n < 3: + return 1.0 + se = 1.0 / (2.0 * np.sqrt(n - 1)) + z = abs(glk_percent / 100.0 - 0.5) / se + return float(2.0 * (1.0 - stats.norm.cdf(z))) def cross_correlation( @@ -295,6 +311,88 @@ def find_best_match( return results[:n_best] +@dataclass +class MissingRingHint: + """A heuristic suggestion that a sample has a missing or false ring.""" + + index: int # sample index where the discontinuity is suggested + kind: str # "missing" (sample omits a ring) or "false" (sample has an extra) + baseline_correlation: float + improved_correlation: float + + @property + def gain(self) -> float: + return self.improved_correlation - self.baseline_correlation + + +def _safe_corr(a: np.ndarray, b: np.ndarray) -> float: + if len(a) < 3 or len(b) < 3 or len(a) != len(b): + return 0.0 + if np.std(a) == 0 or np.std(b) == 0: + return 0.0 + return float(np.corrcoef(a, b)[0, 1]) + + +def detect_missing_ring( + sample: np.ndarray, + reference: np.ndarray, + min_segment: int = 15, + min_gain: float = 0.1, +) -> Optional[MissingRingHint]: + """ + Heuristically locate a single missing or false ring causing a dating misfit. + + A locally absent ring (a year with no measurable growth) or a false/double + ring shifts every later ring by one, so a sample can correlate only weakly + even at its true position. This scans candidate break points: at each one it + tests whether re-aligning the later part of the sample by +/-1 ring against + the reference substantially improves correlation over the rigid alignment. + + The reference must be at least one ring longer than the sample and aligned so + that ``sample[0]`` corresponds to ``reference[0]``. + + Args: + sample: Standardized sample series (oldest -> newest). + reference: Standardized reference series aligned to the sample start; + should extend at least ``len(sample) + 1``. + min_segment: Minimum rings on each side of a candidate break. + min_gain: Minimum correlation improvement to report a hint. + + Returns: + A ``MissingRingHint`` for the best-supported location, or ``None``. + """ + sample = np.asarray(sample, dtype=np.float64) + reference = np.asarray(reference, dtype=np.float64) + n = len(sample) + if n < 2 * min_segment + 1 or len(reference) < n + 1: + return None + + baseline = _safe_corr(sample, reference[:n]) + best: Optional[MissingRingHint] = None + + for k in range(min_segment, n - min_segment): + head = _safe_corr(sample[:k], reference[:k]) + + # "missing": sample omits a ring at k, so its tail lines up one ring + # later in the reference. + tail_missing = _safe_corr(sample[k:], reference[k + 1 : k + 1 + (n - k)]) + # "false": sample has an extra ring at k, so its tail lines up one ring + # earlier in the reference. + tail_false = _safe_corr(sample[k:], reference[k - 1 : k - 1 + (n - k)]) + + for kind, tail in (("missing", tail_missing), ("false", tail_false)): + combined = 0.5 * (head + tail) + if combined - baseline >= min_gain and (best is None or combined > best.improved_correlation): + best = MissingRingHint( + index=k, + kind=kind, + baseline_correlation=baseline, + improved_correlation=combined, + ) + + return best + + def dating_confidence(result: CorrelationResult) -> str: """ Assess confidence level of a dating result. @@ -307,12 +405,22 @@ def dating_confidence(result: CorrelationResult) -> str: Returns: Confidence level: "HIGH", "MEDIUM", or "LOW". """ - # High confidence: strong correlation and good overlap - if result.t_value >= 6.0 and result.correlation >= 0.55 and result.overlap >= 50: + # High confidence: strong correlation, good overlap, and parallel variation. + if ( + result.t_value >= 6.0 + and result.correlation >= 0.55 + and result.overlap >= 50 + and result.gleichlauf >= 62.0 + ): return "HIGH" - # Medium confidence: decent statistics - if result.t_value >= 4.0 and result.correlation >= 0.45 and result.overlap >= 30: + # Medium confidence: decent statistics. + if ( + result.t_value >= 4.0 + and result.correlation >= 0.45 + and result.overlap >= 30 + and result.gleichlauf >= 58.0 + ): return "MEDIUM" # Low confidence diff --git a/src/dendro/crossdating/detrend.py b/src/dendro/crossdating/detrend.py index 6f893e0..9b01fc9 100644 --- a/src/dendro/crossdating/detrend.py +++ b/src/dendro/crossdating/detrend.py @@ -18,9 +18,7 @@ from typing import Optional import numpy as np -from scipy import signal from scipy.optimize import curve_fit -from scipy.interpolate import UnivariateSpline class DetrendMethod(Enum): @@ -184,33 +182,64 @@ def _fit_cubic_spline( y: np.ndarray, mask: np.ndarray, period: float, + freq_response: float = 0.5, ) -> np.ndarray: """ - Fit a cubic smoothing spline with specified frequency response. + Fit a Cook & Peters (1981) cubic smoothing spline. + + This is the dendrochronology-standard "n-year spline": a cubic smoothing + spline whose amplitude frequency response equals ``freq_response`` (50% by + convention) at a wavelength of ``period`` years. It is implemented as the + Whittaker/Reinsch smoother ``g = (I + lambda * D'D)^-1 y`` where ``D`` is the + second-difference operator, which is mathematically equivalent to the + natural cubic smoothing spline on unit-spaced data. + + The penalty's transfer magnitude is ``(2 - 2 cos w)`` per the second + difference, so the smoother response is ``H(w) = 1 / (1 + lambda (2-2cos w)^2)``. + Setting ``H = freq_response`` at ``w0 = 2*pi/period`` gives + ``lambda = (1/freq_response - 1) / (2 - 2 cos w0)^2``. Args: - period: Wavelength at which 50% of variance is retained. + x: Sample index positions (unused beyond length; data are unit-spaced). + y: Raw ring-width values. + mask: Boolean mask of valid (positive, non-NaN) values. + period: Cutoff wavelength in years (e.g. 2/3 of the series length). + freq_response: Fraction of variance retained at ``period`` (default 0.5). + + Returns: + The fitted growth curve over the full series length (strictly positive). """ - # Convert period to smoothing parameter - # Approximate relationship between period and smoothing + from scipy.sparse import diags, identity + from scipy.sparse.linalg import spsolve + n = len(y) - s = n / (period / 2) # Smoothing factor + yf = np.asarray(y, dtype=np.float64).copy() + + # Linearly interpolate masked-out points so the smoother stays well-posed. + if not np.all(mask): + if np.sum(mask) < 2: + return np.full(n, max(np.nanmean(yf[mask]) if np.any(mask) else 1.0, 0.01)) + idx = np.arange(n) + yf[~mask] = np.interp(idx[~mask], idx[mask], yf[mask]) + + if n < 4: + return np.maximum(np.full(n, np.mean(yf)), 0.01) + + # Keep the cutoff wavelength physically meaningful for the series length. + period = float(np.clip(period, 4.0, max(4.0, 2.0 * n))) + omega0 = 2.0 * np.pi / period + amp = 2.0 - 2.0 * np.cos(omega0) # second-difference transfer magnitude + lam = (1.0 / freq_response - 1.0) / (amp * amp) try: - spline = UnivariateSpline( - x[mask], - y[mask], - s=s * np.var(y[mask]), - k=3, - ) - curve = spline(x) + e = np.ones(n) + d2 = diags([e[:-2], -2.0 * e[1:-1], e[2:]], [0, 1, 2], shape=(n - 2, n)) + system = (identity(n, format="csc") + lam * (d2.transpose() @ d2)).tocsc() + curve = spsolve(system, yf) return np.maximum(curve, 0.01) - except Exception: - # Fall back to linear coeffs = np.polyfit(x[mask], y[mask], 1) - curve = np.polyval(coeffs, x) - return np.maximum(curve, 0.01) + return np.maximum(np.polyval(coeffs, x), 0.01) def standardize( @@ -271,7 +300,6 @@ def prewhiten(values: np.ndarray, order: int = 1) -> np.ndarray: # Fit AR model using Yule-Walker equations try: - from scipy.signal import lfilter # Estimate AR coefficients r = np.correlate(values[valid_mask] - np.mean(values[valid_mask]), diff --git a/src/dendro/crossdating/felling.py b/src/dendro/crossdating/felling.py new file mode 100644 index 0000000..5a5dae4 --- /dev/null +++ b/src/dendro/crossdating/felling.py @@ -0,0 +1,150 @@ +""" +Felling-year interpretation from a dated last ring. + +Dating a sample establishes the calendar year of its *last measured ring*. What +that implies about the felling (cutting) date depends on what survives on the +sample: + +- **Bark edge / waney edge present** -> the last ring is the final year of + growth, so the felling year is known exactly (to the season). +- **Sapwood present but no bark edge** -> some outer sapwood rings are missing. + The felling year is *after* the last ring, by an amount estimated from + species-specific sapwood statistics. The result is a *range*. +- **Heartwood/sapwood boundary only** -> we know the felling year is no earlier + than the last heartwood ring plus the minimum sapwood count: a *terminus post + quem* (felled after). + +This module turns the dated last-ring year plus those observations into an +explicit ``FellingEstimate`` rather than overloading a single "felling year" +number (the previous behaviour, which silently reported the last ring as the +felling year even with no bark edge). + +Sapwood estimates are intentionally coarse and species-grouped; precise sapwood +models (e.g. Hollstein, Baillie & Pilcher, Sohar et al.) should be selected per +region and species for publication-grade work. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from enum import Enum +from typing import Optional + + +class FellingType(str, Enum): + """How precisely the felling date is known.""" + + EXACT = "exact" # bark edge present + RANGE = "range" # sapwood present, estimated range + AFTER = "after" # terminus post quem only + + +# Typical (min, median, max) number of sapwood rings by broad species group. +# Conifers do not have a reliably countable sapwood band for this purpose, so +# only a nominal minimum is offered. Oak figures are commonly used European/NA +# ranges and should be refined per regional model before publication. +SAPWOOD_ESTIMATES: dict[str, tuple[int, int, int]] = { + "QUERCUS": (9, 15, 50), # white/red oak (QUAL, QURU, ...) + "DEFAULT": (10, 20, 40), +} + +# Species codes that belong to the oak group. +_OAK_CODES = {"QUAL", "QURU", "QUVE", "QUPR", "QUST", "QUSP", "QUMO"} + + +def _sapwood_model(species: Optional[str]) -> tuple[int, int, int]: + if species and species.upper() in _OAK_CODES: + return SAPWOOD_ESTIMATES["QUERCUS"] + return SAPWOOD_ESTIMATES["DEFAULT"] + + +@dataclass +class FellingEstimate: + """An interpretation of the felling date from a dated last ring.""" + + last_ring_year: int + felling_type: FellingType + earliest_felling: int + latest_felling: Optional[int] # None for open-ended "felled after" + note: str + + def summary(self) -> str: + if self.felling_type is FellingType.EXACT: + return f"Felling year: {self.earliest_felling} (exact, bark edge present)" + if self.felling_type is FellingType.RANGE: + return ( + f"Felling year range: {self.earliest_felling}-{self.latest_felling} " + "(estimated from sapwood)" + ) + return f"Felled after {self.earliest_felling} (terminus post quem)" + + def to_dict(self) -> dict: + return { + "last_ring_year": self.last_ring_year, + "felling_type": self.felling_type.value, + "earliest_felling": self.earliest_felling, + "latest_felling": self.latest_felling, + "note": self.note, + } + + +def estimate_felling( + last_ring_year: int, + has_bark_edge: bool, + has_sapwood: bool = False, + sapwood_count: Optional[int] = None, + species: Optional[str] = None, +) -> FellingEstimate: + """ + Interpret a dated last ring as a felling date. + + Args: + last_ring_year: Calendar year of the outermost measured ring. + has_bark_edge: True if bark/waney edge is present (exact felling year). + has_sapwood: True if (incomplete) sapwood is present but no bark edge. + sapwood_count: Number of sapwood rings already measured/present, if known. + species: Species code, used to pick a sapwood model. + + Returns: + A ``FellingEstimate`` describing the felling date and its uncertainty. + """ + if has_bark_edge: + return FellingEstimate( + last_ring_year=last_ring_year, + felling_type=FellingType.EXACT, + earliest_felling=last_ring_year, + latest_felling=last_ring_year, + note="Bark/waney edge present: felling year is exact.", + ) + + sap_min, sap_med, sap_max = _sapwood_model(species) + present = sapwood_count or 0 + + if has_sapwood: + # Some sapwood present: estimate the missing remainder to bark. + missing_min = max(0, sap_min - present) + missing_max = max(missing_min, sap_max - present) + earliest = last_ring_year + missing_min + latest = last_ring_year + missing_max + return FellingEstimate( + last_ring_year=last_ring_year, + felling_type=FellingType.RANGE, + earliest_felling=earliest, + latest_felling=latest, + note=( + f"Sapwood present ({present} rings); estimated total " + f"{sap_min}-{sap_max} sapwood rings for this species group." + ), + ) + + # No bark edge, no sapwood retained: only a lower bound on the felling year. + return FellingEstimate( + last_ring_year=last_ring_year, + felling_type=FellingType.AFTER, + earliest_felling=last_ring_year + sap_min, + latest_felling=None, + note=( + "No bark edge or sapwood retained: felling year is at least the last " + f"ring plus the minimum sapwood allowance ({sap_min} rings)." + ), + ) diff --git a/src/dendro/crossdating/matcher.py b/src/dendro/crossdating/matcher.py index f4786bf..37a9fcb 100644 --- a/src/dendro/crossdating/matcher.py +++ b/src/dendro/crossdating/matcher.py @@ -20,15 +20,42 @@ import numpy as np import pandas as pd +from ..reference.chronology_index import ChronologyIndex, ChronologyMetadata +from ..reference.tucson_parser import Chronology, RWLFile from .correlator import ( - CorrelationResult, + dating_confidence, find_best_match, + gleichlauf_significance, segment_correlation, - dating_confidence, ) -from .detrend import detrend_series, standardize, DetrendMethod -from ..reference.chronology_index import ChronologyIndex, ChronologyMetadata -from ..reference.tucson_parser import parse_rwl_file, parse_crn_file, Chronology, RWLFile +from .detrend import DetrendMethod, build_chronology, detrend_series, standardize +from .felling import FellingEstimate, estimate_felling + + +def normalize_orientation(values: np.ndarray, orientation: str) -> np.ndarray: + """ + Return a ring-width series in canonical order (oldest -> newest / pith -> bark). + + All cross-dating in this package assumes the sample runs oldest ring first, + matching how ITRDB reference chronologies are stored. A sample measured from + the bark inward must be reversed before dating. + + Args: + values: Ring-width series. + orientation: ``"pith_to_bark"`` (already canonical) or ``"bark_to_pith"`` + (will be reversed). + + Returns: + The series in oldest-to-newest order. + """ + values = np.asarray(values, dtype=np.float64) + if orientation == "pith_to_bark": + return values + if orientation == "bark_to_pith": + return values[::-1].copy() + raise ValueError( + f"Unknown orientation {orientation!r}; expected 'pith_to_bark' or 'bark_to_pith'." + ) @dataclass @@ -46,11 +73,22 @@ class MatchResult: overlap: int gleichlauf: float confidence: str + glk_p_value: float = 1.0 segment_correlations: list[tuple[int, float, float]] = field(default_factory=list) + @property + def last_ring_year(self) -> int: + """Calendar year of the outermost (last measured) ring.""" + return self.proposed_end_year + @property def felling_year(self) -> int: - """The proposed felling year (last ring + bark edge).""" + """ + Backwards-compatible alias for the last-ring year. + + Note: this is only the felling year when a bark edge is present. Use the + report's ``felling_estimate`` for the correct interpretation otherwise. + """ return self.proposed_end_year @@ -65,27 +103,32 @@ class CrossdateReport: matches: list[MatchResult] consensus_year: Optional[int] = None consensus_confidence: str = "LOW" + felling_estimate: Optional[FellingEstimate] = None warnings: list[str] = field(default_factory=list) def to_dict(self) -> dict: - """Convert report to dictionary for serialization.""" + """Convert report to a JSON-serializable dictionary (native types only).""" return { "sample_name": self.sample_name, - "sample_length": self.sample_length, - "has_bark_edge": self.has_bark_edge, + "sample_length": int(self.sample_length), + "has_bark_edge": bool(self.has_bark_edge), "detrend_method": self.detrend_method, - "consensus_year": self.consensus_year, + "consensus_year": int(self.consensus_year) if self.consensus_year is not None else None, "consensus_confidence": self.consensus_confidence, - "warnings": self.warnings, + "felling_estimate": self.felling_estimate.to_dict() if self.felling_estimate else None, + "warnings": list(self.warnings), "matches": [ { "reference": m.reference_name, "species": m.reference_species, "state": m.reference_state, - "felling_year": m.felling_year, - "correlation": round(m.correlation, 4), - "t_value": round(m.t_value, 2), - "overlap": m.overlap, + "last_ring_year": int(m.last_ring_year), + "felling_year": int(m.felling_year), + "correlation": round(float(m.correlation), 4), + "t_value": round(float(m.t_value), 2), + "gleichlauf": round(float(m.gleichlauf), 1), + "glk_p_value": float(m.glk_p_value), + "overlap": int(m.overlap), "confidence": m.confidence, } for m in self.matches @@ -117,6 +160,10 @@ def __init__( else: self.index = ChronologyIndex() + # Cache of standardized master chronologies keyed by reference filepath, + # so each reference is parsed and detrended at most once per matcher. + self._master_cache: dict[str, Optional[tuple[np.ndarray, int]]] = {} + def date_sample( self, values: np.ndarray, @@ -129,12 +176,17 @@ def date_sample( detrend_method: DetrendMethod = DetrendMethod.SPLINE, min_overlap: int = 30, max_references: int = 20, + orientation: str = "pith_to_bark", + has_sapwood: bool = False, + sapwood_count: Optional[int] = None, ) -> CrossdateReport: """ Cross-date a sample against available reference chronologies. Args: - values: Ring width measurements (from bark to pith, outer to inner). + values: Ring-width measurements. By default interpreted as + oldest-to-newest (``orientation="pith_to_bark"``); pass + ``orientation="bark_to_pith"`` for a series measured bark inward. sample_name: Identifier for the sample. has_bark_edge: Whether sample includes the outermost ring (bark). species_filter: Only match against these species. @@ -144,11 +196,16 @@ def date_sample( detrend_method: Method for removing growth trend. min_overlap: Minimum overlapping years required. max_references: Maximum number of references to try. + orientation: Order of ``values`` (``"pith_to_bark"`` or + ``"bark_to_pith"``); normalized internally to oldest-first. + has_sapwood: Whether incomplete sapwood is present (no bark edge). + sapwood_count: Number of sapwood rings present, if known. Returns: CrossdateReport with all match results and assessment. """ - values = np.asarray(values, dtype=np.float64) + # Normalize to canonical oldest-to-newest order before anything else. + values = normalize_orientation(values, orientation) # Validate input if len(values) < min_overlap: @@ -217,7 +274,7 @@ def date_sample( ) if match is not None: all_matches.append(match) - except Exception as e: + except Exception: continue # Skip problematic references # Sort by t-value @@ -235,6 +292,22 @@ def date_sample( # Assess consensus self._assess_consensus(report) + # If the best match is imperfect, check whether a single missing/false + # ring would explain the misfit, and point the user at its location. + if report.matches and report.matches[0].correlation < 0.7: + self._add_missing_ring_hint(report, sample_std) + + # Interpret the dated last ring as a felling date. + if report.consensus_year is not None: + best_species = report.matches[0].reference_species if report.matches else None + report.felling_estimate = estimate_felling( + last_ring_year=report.consensus_year, + has_bark_edge=has_bark_edge, + has_sapwood=has_sapwood, + sapwood_count=sapwood_count, + species=best_species, + ) + return report def _match_against_reference( @@ -247,40 +320,10 @@ def _match_against_reference( ) -> Optional[MatchResult]: """Match sample against a single reference chronology.""" - # Load the reference data - data = self.index.load_chronology(meta) - - if data is None: - return None - - # Get reference values and years - if isinstance(data, Chronology): - ref_values = data.values - ref_start = data.start_year - - # Detrend/standardize reference (CRN files are often already indexed) - if np.mean(ref_values) > 500: # Likely scaled to 1000 - ref_std = ref_values / 1000.0 - else: - ref_std = standardize(ref_values) - elif isinstance(data, RWLFile): - # For RWL files, build a master chronology - if not data.series: - return None - - # Get chronology - df = data.to_dataframe() - ref_values = df.mean(axis=1).values - ref_start = int(df.index.min()) - - # Detrend/standardize - try: - detrended, _ = detrend_series(ref_values) - ref_std = standardize(detrended) - except Exception: - ref_std = standardize(ref_values) - else: + master = self._get_master(meta) + if master is None: return None + ref_std, ref_start = master # Find best match best_matches = find_best_match( @@ -292,9 +335,10 @@ def _match_against_reference( best = best_matches[0] - # Filter by era - sample_end = best.position + len(sample) - 1 - if best.position < era_start - len(sample) or sample_end > era_end + 50: + # Filter by era: the dated last ring (felling year candidate) must fall + # within the requested window. + last_ring_year = best.position + len(sample) - 1 + if not (era_start <= last_ring_year <= era_end): return None # Calculate segment correlations for detailed analysis @@ -322,9 +366,92 @@ def _match_against_reference( overlap=best.overlap, gleichlauf=best.gleichlauf, confidence=dating_confidence(best), + glk_p_value=gleichlauf_significance(best.gleichlauf, best.overlap), segment_correlations=seg_corrs, ) + def _get_master(self, meta: ChronologyMetadata) -> Optional[tuple[np.ndarray, int]]: + """ + Return a standardized master chronology for a reference, building it once. + + RWL files are turned into a site master by detrending **each** series to + a dimensionless index and then taking a robust (biweight) mean of the + indices -- not by averaging raw widths, which would leave each tree's + age-related growth trend in the master and produce spurious matches. + + Args: + meta: Index entry for the reference file. + + Returns: + ``(standardized_values, start_year)`` or ``None`` if it cannot be built. + """ + cached = self._master_cache.get(meta.filepath) + if cached is not None or meta.filepath in self._master_cache: + return cached + + result: Optional[tuple[np.ndarray, int]] = None + data = self.index.load_chronology(meta) + + if isinstance(data, Chronology) and data.values is not None and len(data.values): + # CRN files are already standardized indices (often scaled by 1000). + ref_values = np.asarray(data.values, dtype=np.float64) + ref_std = ref_values / 1000.0 if np.nanmean(ref_values) > 500 else ref_values + result = (ref_std, data.start_year) + + elif isinstance(data, RWLFile) and data.series: + series_list: list[np.ndarray] = [] + years_list: list[np.ndarray] = [] + for s in data.series.values(): + if s.length < 10: + continue + try: + detrended, _ = detrend_series(s.values, method=DetrendMethod.SPLINE) + series_list.append(standardize(detrended, method="ratio")) + years_list.append(s.years) + except Exception: + continue + if series_list: + years, chron, depth = build_chronology( + series_list, years_list, method="biweight" + ) + valid = ~np.isnan(chron) + if np.any(valid): + first = int(np.argmax(valid)) + last = len(chron) - int(np.argmax(valid[::-1])) + result = (chron[first:last], int(years[first])) + + self._master_cache[meta.filepath] = result + return result + + def _add_missing_ring_hint(self, report: CrossdateReport, sample_std: np.ndarray): + """Add a warning if a missing/false ring likely explains a weak match.""" + from .correlator import detect_missing_ring + + best = report.matches[0] + meta = next( + (m for m in self.index.entries if m.site_name == best.reference_name), None + ) + if meta is None: + return + master = self._get_master(meta) + if master is None: + return + ref_std, ref_start = master + offset = best.proposed_start_year - ref_start + if offset < 0 or offset + len(sample_std) + 1 > len(ref_std): + return + hint = detect_missing_ring( + sample_std, ref_std[offset : offset + len(sample_std) + 20] + ) + if hint is not None and hint.gain >= 0.1: + year = best.proposed_start_year + hint.index + report.warnings.append( + f"Possible {hint.kind} ring near year {year} (sample ring " + f"{hint.index + 1}): re-aligning there raises correlation from " + f"{hint.baseline_correlation:.2f} to {hint.improved_correlation:.2f}. " + "Re-examine the measurements in that region." + ) + def _assess_consensus(self, report: CrossdateReport): """Assess consensus among top matches.""" @@ -401,11 +528,11 @@ def _assess_consensus(self, report: CrossdateReport): "Check measurements in these regions." ) - # Geographic diversity warning - states = set(m.reference_state for m in good_matches[:5]) + # Geographic diversity warning (only when state metadata is available). + states = {m.reference_state for m in good_matches[:5] if m.reference_state} if len(states) == 1: report.warnings.append( - f"All top matches from single state ({list(states)[0]}). " + f"All top matches from single state ({next(iter(states))}). " "Geographic diversity would strengthen confidence." ) diff --git a/src/dendro/imaging/__init__.py b/src/dendro/imaging/__init__.py index 302ccb5..d46161e 100644 --- a/src/dendro/imaging/__init__.py +++ b/src/dendro/imaging/__init__.py @@ -1,7 +1,7 @@ """Image processing for ring width extraction from scanned samples.""" -from .ring_detector import detect_rings from .path_sampler import sample_along_path +from .ring_detector import detect_rings from .viewer import MeasurementViewer __all__ = [ diff --git a/src/dendro/imaging/path_sampler.py b/src/dendro/imaging/path_sampler.py index 449fd99..33afc7b 100644 --- a/src/dendro/imaging/path_sampler.py +++ b/src/dendro/imaging/path_sampler.py @@ -256,8 +256,9 @@ def widths_to_tucson( Returns: Tucson format string. """ - # Convert mm to 0.01mm units (standard Tucson) - widths_001mm = (widths_mm * 100).astype(int) + # Convert mm to 0.01mm units (standard Tucson). Round rather than truncate + # so a 1.27 mm ring becomes 127, not 126. + widths_001mm = np.rint(np.asarray(widths_mm, dtype=np.float64) * 100).astype(int) # Reverse if needed (Tucson goes from oldest to newest) # Our widths are bark-to-pith (newest to oldest), so reverse @@ -276,10 +277,7 @@ def widths_to_tucson( idx = 0 while idx < len(widths_001mm): - # Start of decade - decade_start = (year // 10) * 10 - - # How many values fit in this decade + # How many values fit in the current decade row values_in_decade = min(10 - (year % 10), len(widths_001mm) - idx) # Build line diff --git a/src/dendro/imaging/ring_detector.py b/src/dendro/imaging/ring_detector.py index c848b5e..d6f2c35 100644 --- a/src/dendro/imaging/ring_detector.py +++ b/src/dendro/imaging/ring_detector.py @@ -248,14 +248,13 @@ def _calculate_gradient(profile: np.ndarray) -> np.ndarray: mode='same' ) - # Calculate gradient (first derivative) - gradient = np.gradient(smoothed) - - # We want to detect transitions from light to dark (ring boundaries) - # Take absolute value to catch both directions - gradient = np.abs(gradient) - - return gradient + # Signed first derivative. A tree-ring (year) boundary is a *single* + # transition -- the abrupt step from dense dark latewood at the end of one + # year to light earlywood at the start of the next. Returning the signed + # gradient lets the boundary finder pick that one polarity, instead of + # taking the absolute value (which fires on both the earlywood->latewood and + # latewood->earlywood edges and double-counts every ring). + return np.gradient(smoothed) def _find_boundaries( @@ -276,25 +275,26 @@ def _find_boundaries( Returns: List of boundary positions (pixel indices). """ - # Normalize gradient - grad_max = np.max(gradient) - if grad_max == 0: - return [] + gradient = np.asarray(gradient, dtype=np.float64) - grad_norm = gradient / grad_max + # Pick the transition polarity that carries more energy and detect peaks of + # that single direction only (avoids counting both edges of each ring). + pos_energy = float(np.sum(gradient[gradient > 0])) + neg_energy = float(-np.sum(gradient[gradient < 0])) + signal = gradient if pos_energy >= neg_energy else -gradient - # Find peaks above threshold - threshold_value = threshold * 0.5 # Scale threshold + grad_max = np.max(signal) + if grad_max <= 0: + return [] - # Simple peak finding - peaks = [] - for i in range(min_distance, len(gradient) - min_distance): - # Check if this is a local maximum - window = gradient[max(0, i - min_distance // 2): - min(len(gradient), i + min_distance // 2 + 1)] + grad_norm = signal / grad_max + threshold_value = threshold * 0.5 # Scale threshold - if gradient[i] == np.max(window) and grad_norm[i] > threshold_value: - # Check distance from last peak + peaks: list[int] = [] + for i in range(min_distance, len(signal) - min_distance): + window = signal[max(0, i - min_distance // 2): + min(len(signal), i + min_distance // 2 + 1)] + if signal[i] == np.max(window) and grad_norm[i] > threshold_value: if not peaks or (i - peaks[-1]) >= min_distance: peaks.append(i) diff --git a/src/dendro/imaging/viewer.py b/src/dendro/imaging/viewer.py index 92ee1c1..37c4324 100644 --- a/src/dendro/imaging/viewer.py +++ b/src/dendro/imaging/viewer.py @@ -12,14 +12,13 @@ from dataclasses import dataclass, field from pathlib import Path -from typing import Optional, Callable +from typing import Callable, Optional import numpy as np try: import matplotlib.pyplot as plt - from matplotlib.widgets import Button, Slider - from matplotlib.backend_bases import MouseEvent, KeyEvent + from matplotlib.backend_bases import KeyEvent, MouseEvent # noqa: F401 HAS_MATPLOTLIB = True except ImportError: HAS_MATPLOTLIB = False @@ -360,7 +359,7 @@ def _finalize(self): self.session.ring_widths_mm = widths self.session.is_finalized = True - print(f"\nMeasurement complete!") + print("\nMeasurement complete!") print(f"Number of rings: {len(widths)}") print(f"Total span: {boundaries_mm[-1] - boundaries_mm[0]:.1f} mm") print(f"Mean ring width: {np.mean(widths):.2f} mm") diff --git a/src/dendro/reference/__init__.py b/src/dendro/reference/__init__.py index a773477..65766f5 100644 --- a/src/dendro/reference/__init__.py +++ b/src/dendro/reference/__init__.py @@ -1,8 +1,8 @@ """Reference chronology handling - ITRDB data download and parsing.""" -from .tucson_parser import parse_rwl_file, parse_crn_file -from .downloader import download_chronologies from .chronology_index import ChronologyIndex +from .downloader import download_chronologies +from .tucson_parser import parse_crn_file, parse_rwl_file __all__ = [ "parse_rwl_file", diff --git a/src/dendro/reference/chronology_index.py b/src/dendro/reference/chronology_index.py index 61f670a..6a96804 100644 --- a/src/dendro/reference/chronology_index.py +++ b/src/dendro/reference/chronology_index.py @@ -6,11 +6,11 @@ import json import re -from dataclasses import dataclass, asdict +from dataclasses import asdict, dataclass from pathlib import Path from typing import Optional -from .tucson_parser import parse_rwl_file, parse_crn_file, Chronology, RWLFile +from .tucson_parser import Chronology, RWLFile, parse_crn_file, parse_rwl_file @dataclass diff --git a/src/dendro/reference/downloader.py b/src/dendro/reference/downloader.py index 95147ca..2d922dd 100644 --- a/src/dendro/reference/downloader.py +++ b/src/dendro/reference/downloader.py @@ -17,7 +17,6 @@ import requests from tqdm import tqdm - # Base URLs for ITRDB data ITRDB_MEASUREMENTS_BASE = "https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/usa/" ITRDB_CHRONOLOGIES_BASE = "https://www.ncei.noaa.gov/pub/data/paleo/treering/chronologies/northamerica/usa/" diff --git a/src/dendro/visualization/__init__.py b/src/dendro/visualization/__init__.py index 4e4498f..1a78421 100644 --- a/src/dendro/visualization/__init__.py +++ b/src/dendro/visualization/__init__.py @@ -9,10 +9,10 @@ """ from .plots import ( - plot_crossdate_results, plot_correlation_profile, - plot_segment_analysis, + plot_crossdate_results, plot_sample_reference_overlay, + plot_segment_analysis, save_diagnostic_plots, ) diff --git a/src/dendro/visualization/plots.py b/src/dendro/visualization/plots.py index d879053..bc27b59 100644 --- a/src/dendro/visualization/plots.py +++ b/src/dendro/visualization/plots.py @@ -8,13 +8,12 @@ from pathlib import Path from typing import Optional -import numpy as np -import matplotlib.pyplot as plt import matplotlib.patches as mpatches +import matplotlib.pyplot as plt +import numpy as np -from ..crossdating.matcher import CrossdateReport, MatchResult -from ..crossdating.correlator import sliding_correlation, CorrelationResult - +from ..crossdating.correlator import CorrelationResult +from ..crossdating.matcher import CrossdateReport # Known marker years for New England region MARKER_YEARS = { @@ -206,7 +205,7 @@ def plot_segment_analysis( # Plot as bars width = (years[1] - years[0]) * 0.8 if len(years) > 1 else 20 - bars = ax.bar(years, correlations, width=width, color=colors, alpha=0.7, edgecolor='black') + ax.bar(years, correlations, width=width, color=colors, alpha=0.7, edgecolor='black') # Add threshold lines ax.axhline(0.55, color='green', linestyle='--', alpha=0.7, label='HIGH (r=0.55)') diff --git a/tests/fixtures/generate_fixtures.py b/tests/fixtures/generate_fixtures.py new file mode 100644 index 0000000..5784561 --- /dev/null +++ b/tests/fixtures/generate_fixtures.py @@ -0,0 +1,175 @@ +""" +Deterministic synthetic ITRDB-style fixtures for testing the dating pipeline. + +NOAA's ITRDB archive is not reachable from CI, and redistributing real +chronologies raises licensing questions, so we generate realistic, +*deterministic* tree-ring data instead. The generated series share a common +"climate" signal (an AR(1) process plus sharp marker years), each modulated by +an individual age-related growth trend and measurement noise -- exactly the +structure real cross-dating relies on. + +Running this module rewrites the committed fixtures under +``tests/fixtures/reference/synth`` and a known-date sample CSV. The output is +fully determined by the seeds below, so regenerating must produce identical +files (this is asserted by ``tests/test_fixtures.py``). +""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np + +FIXTURE_DIR = Path(__file__).parent +REF_DIR = FIXTURE_DIR / "reference" / "synth" +SAMPLE_DIR = FIXTURE_DIR / "samples" + +# Calendar span of the synthetic "climate" +CLIMATE_START = 1500 +CLIMATE_END = 2000 + +# Known felling year of the committed dating sample (with bark edge). +SAMPLE_FELLING_YEAR = 1789 +SAMPLE_N_RINGS = 90 + + +def make_climate_signal(seed: int = 20240101) -> np.ndarray: + """An AR(1) climate signal with a few sharp marker years.""" + rng = np.random.default_rng(seed) + n = CLIMATE_END - CLIMATE_START + 1 + phi = 0.3 + c = np.zeros(n) + for i in range(1, n): + c[i] = phi * c[i - 1] + rng.normal(0, 1) + c = (c - c.mean()) / c.std() + # Sharp marker years (e.g. 1816 "Year Without Summer"): very narrow rings. + for year in (1816, 1709, 1780, 1601): + idx = year - CLIMATE_START + if 0 <= idx < n: + c[idx] = -3.0 + return c + + +def _neg_exp_trend(n: int, a: float, b: float, c: float) -> np.ndarray: + """Age-related growth trend: large when young, decaying with age.""" + t = np.arange(n) + return a * np.exp(-b * t) + c + + +def synth_core( + climate: np.ndarray, + start_year: int, + n_rings: int, + rng: np.random.Generator, + sensitivity: float = 0.45, +) -> np.ndarray: + """Generate one raw ring-width series (in mm) from the shared climate.""" + s = start_year - CLIMATE_START + sig = climate[s : s + n_rings] + trend = _neg_exp_trend( + n_rings, + a=rng.uniform(1.5, 3.0), + b=rng.uniform(0.01, 0.03), + c=rng.uniform(0.4, 0.8), + ) + noise = rng.normal(0, 0.12, n_rings) + widths_mm = trend * (1.0 + sensitivity * sig + noise) + return np.clip(widths_mm, 0.05, None) + + +def widths_to_tucson_lines(series_id: str, start_year: int, widths_mm: np.ndarray) -> list[str]: + """Format a series as standard Tucson .rwl lines (0.01 mm units, 999 stop).""" + vals = np.rint(widths_mm * 100).astype(int) + vals = np.clip(vals, 1, 9989) # keep away from the 999/9990 sentinels + sid = series_id[:8].ljust(8) + lines: list[str] = [] + year = start_year + idx = 0 + n = len(vals) + while idx < n: + count = min(10 - (year % 10), n - idx) + line = f"{sid}{year:4d}" + for i in range(count): + line += f"{int(vals[idx + i]):6d}" + idx += count + year += count + if idx >= n: + line += " 999" + lines.append(line) + return lines + + +def write_site( + path: Path, + site_code: str, + species: str, + climate: np.ndarray, + cores: list[tuple[str, int, int]], + seed: int, +) -> None: + """Write one .rwl site file with a header and several cores.""" + rng = np.random.default_rng(seed) + lines = [ + f"{site_code} 1 Synthetic Site {site_code} {species}", + f"{site_code} 2 Test Region {species} 45 00 N 72 00 W 300M", + f"{site_code} 3 Synthetic generator 1500 2000", + ] + for core_id, start, length in cores: + widths = synth_core(climate, start, length, rng) + lines.extend(widths_to_tucson_lines(core_id, start, widths)) + path.write_text("\n".join(lines) + "\n") + + +def generate_all() -> None: + REF_DIR.mkdir(parents=True, exist_ok=True) + SAMPLE_DIR.mkdir(parents=True, exist_ok=True) + climate = make_climate_signal() + + # Two synthetic sites of eastern-hemlock-like and white-pine-like cores. + write_site( + REF_DIR / "synth01.rwl", + site_code="SYNTH01", + species="TSCA", + climate=climate, + cores=[ + ("SY01A", 1600, 300), + ("SY01B", 1620, 280), + ("SY01C", 1650, 250), + ("SY01D", 1580, 320), + ("SY01E", 1700, 200), + ("SY01F", 1660, 240), + ("SY01G", 1550, 350), + ("SY01H", 1690, 210), + ], + seed=111, + ) + write_site( + REF_DIR / "synth02.rwl", + site_code="SYNTH02", + species="PIST", + climate=climate, + cores=[ + ("SY02A", 1610, 290), + ("SY02B", 1630, 270), + ("SY02C", 1640, 260), + ("SY02D", 1600, 300), + ("SY02E", 1680, 220), + ("SY02F", 1655, 245), + ], + seed=222, + ) + + # A known-date dating sample: a core that ends (bark edge) at the felling year. + rng = np.random.default_rng(999) + sample_start = SAMPLE_FELLING_YEAR - SAMPLE_N_RINGS + 1 + sample_widths = synth_core(climate, sample_start, SAMPLE_N_RINGS, rng) + # CSV is written oldest -> newest (canonical / pith -> bark orientation). + csv_lines = ["year,width_mm"] + for i, w in enumerate(sample_widths): + csv_lines.append(f"{sample_start + i},{w:.3f}") + (SAMPLE_DIR / "known_1789.csv").write_text("\n".join(csv_lines) + "\n") + + +if __name__ == "__main__": + generate_all() + print(f"Wrote fixtures to {REF_DIR} and {SAMPLE_DIR}") diff --git a/tests/fixtures/reference/synth/synth01.rwl b/tests/fixtures/reference/synth/synth01.rwl new file mode 100644 index 0000000..20dca62 --- /dev/null +++ b/tests/fixtures/reference/synth/synth01.rwl @@ -0,0 +1,218 @@ +SYNTH01 1 Synthetic Site SYNTH01 TSCA +SYNTH01 2 Test Region TSCA 45 00 N 72 00 W 300M +SYNTH01 3 Synthetic generator 1500 2000 +SY01A 1600 215 5 526 311 224 243 110 187 129 147 +SY01A 1610 170 98 235 91 170 164 253 316 262 204 +SY01A 1620 94 161 287 232 205 279 201 246 239 90 +SY01A 1630 223 344 243 155 162 136 254 138 139 23 +SY01A 1640 127 156 147 52 70 156 167 116 138 124 +SY01A 1650 116 66 96 181 153 129 100 121 165 224 +SY01A 1660 37 59 145 236 131 251 154 69 133 112 +SY01A 1670 132 177 159 66 96 103 59 186 207 108 +SY01A 1680 169 73 118 180 167 110 154 111 60 138 +SY01A 1690 104 142 92 160 78 28 140 76 81 5 +SY01A 1700 74 85 113 121 105 44 121 119 124 5 +SY01A 1710 114 163 137 134 124 109 78 124 144 126 +SY01A 1720 89 83 23 53 127 181 97 45 133 118 +SY01A 1730 144 99 73 21 43 48 79 134 147 167 +SY01A 1740 125 183 103 110 8 36 5 79 133 127 +SY01A 1750 169 128 158 110 60 19 65 118 126 139 +SY01A 1760 152 55 55 37 133 97 41 94 111 47 +SY01A 1770 15 9 98 107 65 109 146 148 100 5 +SY01A 1780 5 58 64 105 70 86 59 68 86 75 +SY01A 1790 76 50 62 130 55 74 96 78 87 73 +SY01A 1800 92 102 35 84 74 5 23 124 124 84 +SY01A 1810 94 85 91 153 53 69 5 39 11 90 +SY01A 1820 61 44 65 36 64 8 75 92 105 36 +SY01A 1830 41 23 119 129 144 15 50 107 53 57 +SY01A 1840 64 49 98 79 68 73 70 41 107 106 +SY01A 1850 99 71 122 74 84 49 43 57 62 99 +SY01A 1860 75 35 47 115 62 75 22 50 66 70 +SY01A 1870 109 100 37 63 58 48 23 44 44 86 +SY01A 1880 66 87 142 110 70 68 112 66 28 27 +SY01A 1890 74 68 48 99 31 22 36 5 27 41 999 +SY01B 1620 68 215 395 346 272 383 342 358 309 117 +SY01B 1630 265 455 359 272 185 173 339 176 208 44 +SY01B 1640 139 232 155 92 116 151 211 164 172 158 +SY01B 1650 129 89 77 224 176 145 145 134 130 229 +SY01B 1660 31 62 127 264 135 262 138 56 87 109 +SY01B 1670 146 212 164 68 98 68 46 156 219 125 +SY01B 1680 153 58 105 167 133 111 125 92 46 105 +SY01B 1690 107 95 55 127 80 7 132 43 73 5 +SY01B 1700 84 41 94 96 95 33 82 122 101 5 +SY01B 1710 88 116 136 90 98 96 56 111 135 92 +SY01B 1720 63 78 8 85 99 147 96 28 109 97 +SY01B 1730 148 74 78 17 37 43 49 114 135 143 +SY01B 1740 98 138 95 81 28 22 19 54 126 100 +SY01B 1750 152 124 131 82 38 6 61 90 124 118 +SY01B 1760 113 38 55 50 118 87 37 85 95 40 +SY01B 1770 28 41 115 87 46 83 132 130 89 11 +SY01B 1780 5 50 48 78 69 59 58 71 75 47 +SY01B 1790 74 52 73 101 55 79 81 87 79 75 +SY01B 1800 73 59 30 64 52 5 15 98 127 72 +SY01B 1810 69 74 95 139 26 62 5 38 17 81 +SY01B 1820 46 54 60 50 53 6 79 80 99 12 +SY01B 1830 34 27 131 107 135 38 50 109 73 75 +SY01B 1840 51 53 85 60 69 80 71 41 95 87 +SY01B 1850 107 59 120 61 92 42 55 56 56 105 +SY01B 1860 85 32 50 117 50 78 36 36 72 73 +SY01B 1870 100 94 47 60 71 54 36 46 41 89 +SY01B 1880 51 70 136 122 68 75 110 77 36 43 +SY01B 1890 81 69 61 83 20 27 50 5 38 64 999 +SY01C 1650 218 228 141 346 349 312 230 204 407 611 +SY01C 1660 70 67 267 555 251 522 295 79 253 265 +SY01C 1670 309 447 251 183 203 194 208 327 452 243 +SY01C 1680 349 158 238 359 290 212 294 247 115 288 +SY01C 1690 203 209 150 329 179 20 258 124 166 5 +SY01C 1700 141 130 242 227 171 51 216 254 275 5 +SY01C 1710 212 185 231 161 213 140 132 244 228 180 +SY01C 1720 160 126 35 130 190 281 191 92 196 180 +SY01C 1730 251 173 134 28 60 76 100 201 263 277 +SY01C 1740 170 268 140 142 41 32 46 82 195 187 +SY01C 1750 273 188 221 123 101 5 82 147 187 183 +SY01C 1760 198 45 60 74 200 129 57 106 120 68 +SY01C 1770 42 20 116 120 77 119 220 202 131 5 +SY01C 1780 5 64 69 128 56 108 106 93 91 46 +SY01C 1790 79 67 85 141 74 91 127 105 106 111 +SY01C 1800 98 113 49 95 76 9 19 112 129 99 +SY01C 1810 96 95 110 155 41 68 5 34 19 103 +SY01C 1820 72 51 78 58 72 5 89 84 105 14 +SY01C 1830 41 5 142 128 137 34 53 80 63 74 +SY01C 1840 59 22 100 97 96 100 90 46 99 130 +SY01C 1850 118 68 112 69 90 43 65 69 57 114 +SY01C 1860 77 45 59 117 55 90 29 52 75 77 +SY01C 1870 122 96 60 65 74 62 13 37 43 80 +SY01C 1880 71 84 128 118 67 70 111 74 27 68 +SY01C 1890 87 64 62 101 26 25 61 5 30 44 999 +SY01D 1580 261 186 141 328 241 519 582 226 91 232 +SY01D 1590 215 164 291 175 17 154 354 131 235 173 +SY01D 1600 210 5 414 297 170 246 117 133 146 155 +SY01D 1610 134 119 226 90 169 153 237 321 223 216 +SY01D 1620 101 172 246 217 135 210 177 208 180 58 +SY01D 1630 156 277 250 156 129 111 241 147 123 5 +SY01D 1640 71 131 94 56 72 114 112 97 87 95 +SY01D 1650 78 54 50 126 120 100 73 111 159 216 +SY01D 1660 27 31 103 208 122 221 82 53 87 65 +SY01D 1670 117 160 115 60 66 38 63 161 181 81 +SY01D 1680 143 73 67 150 116 85 116 102 55 106 +SY01D 1690 88 98 73 144 72 30 101 55 75 5 +SY01D 1700 46 30 94 93 79 44 74 108 99 5 +SY01D 1710 95 91 102 101 94 77 66 119 110 83 +SY01D 1720 94 67 5 72 119 152 79 41 113 80 +SY01D 1730 109 80 65 6 41 19 69 102 139 164 +SY01D 1740 102 147 83 86 15 25 9 54 126 111 +SY01D 1750 155 120 130 71 47 5 59 105 111 103 +SY01D 1760 122 32 35 40 104 75 29 75 88 63 +SY01D 1770 33 30 94 83 63 88 138 127 89 5 +SY01D 1780 5 53 41 85 53 81 67 66 67 61 +SY01D 1790 63 49 51 91 50 67 104 71 85 70 +SY01D 1800 68 67 24 66 56 7 14 88 107 74 +SY01D 1810 78 65 78 141 25 51 5 24 5 80 +SY01D 1820 70 38 47 48 47 5 62 79 79 27 +SY01D 1830 40 8 122 117 111 20 37 94 48 56 +SY01D 1840 47 51 71 60 78 88 77 38 109 111 +SY01D 1850 86 45 93 48 83 54 55 56 58 82 +SY01D 1860 69 39 51 111 49 70 19 53 73 68 +SY01D 1870 98 86 37 52 60 47 30 43 25 78 +SY01D 1880 63 68 121 102 67 68 96 62 31 54 +SY01D 1890 60 57 38 85 29 11 38 5 21 49 999 +SY01E 1700 125 127 351 310 267 118 235 344 334 5 +SY01E 1710 307 354 347 294 299 249 246 359 353 291 +SY01E 1720 231 211 5 181 322 434 288 93 274 297 +SY01E 1730 421 207 213 45 100 113 184 343 367 385 +SY01E 1740 296 359 253 196 39 60 54 140 292 274 +SY01E 1750 388 287 366 221 74 9 153 208 301 298 +SY01E 1760 278 96 129 81 293 148 65 153 173 106 +SY01E 1770 57 72 165 175 93 179 286 297 180 9 +SY01E 1780 5 121 122 192 129 154 136 130 148 102 +SY01E 1790 152 106 112 259 88 137 182 150 157 140 +SY01E 1800 159 155 57 153 98 45 45 192 221 116 +SY01E 1810 144 143 151 265 24 108 5 65 44 162 +SY01E 1820 113 81 83 69 116 5 108 133 127 26 +SY01E 1830 62 20 199 170 201 23 58 166 107 99 +SY01E 1840 101 77 131 113 126 124 114 69 130 148 +SY01E 1850 163 129 146 109 118 64 85 87 82 152 +SY01E 1860 123 72 63 139 74 111 32 63 93 96 +SY01E 1870 163 131 61 78 97 74 30 44 76 138 +SY01E 1880 82 97 165 139 96 99 141 85 25 62 +SY01E 1890 109 99 63 112 29 41 63 13 47 55 999 +SY01F 1660 110 23 178 444 226 492 276 116 183 160 +SY01F 1670 240 303 281 157 149 69 164 261 371 197 +SY01F 1680 287 86 143 260 230 188 190 159 102 220 +SY01F 1690 189 170 121 287 151 41 194 92 169 5 +SY01F 1700 80 75 167 120 127 49 116 184 205 5 +SY01F 1710 141 192 167 162 155 151 95 190 184 132 +SY01F 1720 108 105 5 105 154 215 127 64 121 128 +SY01F 1730 183 90 91 5 34 84 90 170 184 193 +SY01F 1740 132 176 114 142 5 20 7 99 174 121 +SY01F 1750 189 153 178 105 57 5 65 128 139 150 +SY01F 1760 141 47 34 49 151 98 32 83 98 50 +SY01F 1770 23 5 108 76 69 104 139 165 118 5 +SY01F 1780 5 58 67 86 69 73 58 51 85 60 +SY01F 1790 77 51 69 108 58 84 112 106 66 77 +SY01F 1800 76 96 33 75 52 5 12 102 113 81 +SY01F 1810 75 84 96 150 20 66 5 30 15 92 +SY01F 1820 55 51 60 58 58 5 57 92 96 39 +SY01F 1830 45 35 145 115 146 18 49 102 65 51 +SY01F 1840 59 45 83 61 74 74 73 45 114 109 +SY01F 1850 113 70 104 52 86 41 65 53 55 93 +SY01F 1860 90 36 43 112 45 85 30 33 51 79 +SY01F 1870 123 95 48 60 60 52 32 40 39 84 +SY01F 1880 49 68 114 104 64 70 101 76 17 60 +SY01F 1890 78 75 50 80 16 28 54 5 30 45 999 +SY01G 1550 501 240 189 84 107 381 332 243 214 239 +SY01G 1560 190 171 168 216 127 189 286 122 172 249 +SY01G 1570 207 203 160 309 192 171 203 380 43 124 +SY01G 1580 183 89 76 170 153 275 295 112 12 122 +SY01G 1590 112 54 116 101 18 93 131 113 100 134 +SY01G 1600 117 5 207 142 97 113 81 66 53 100 +SY01G 1610 66 53 110 59 69 113 124 160 113 106 +SY01G 1620 28 78 115 94 87 105 85 96 79 28 +SY01G 1630 76 128 107 83 70 47 121 60 67 5 +SY01G 1640 57 80 72 13 42 46 78 66 64 59 +SY01G 1650 45 24 25 76 72 68 48 47 71 120 +SY01G 1660 14 19 44 93 66 129 54 34 66 43 +SY01G 1670 66 84 67 42 50 24 24 84 103 66 +SY01G 1680 86 34 49 95 69 58 65 70 38 67 +SY01G 1690 72 57 29 77 44 14 69 37 43 5 +SY01G 1700 21 32 58 47 50 28 50 85 73 5 +SY01G 1710 74 77 89 57 71 49 57 78 79 61 +SY01G 1720 51 54 6 39 74 101 63 31 63 74 +SY01G 1730 87 57 50 7 22 32 46 68 101 104 +SY01G 1740 87 110 69 53 7 10 14 51 93 79 +SY01G 1750 115 94 101 62 36 5 50 69 86 102 +SY01G 1760 98 20 29 39 91 74 19 56 64 34 +SY01G 1770 18 13 64 71 43 72 114 107 72 7 +SY01G 1780 5 30 39 70 36 61 54 49 54 43 +SY01G 1790 51 43 52 82 36 63 74 59 76 53 +SY01G 1800 68 67 31 62 34 8 14 73 99 59 +SY01G 1810 58 64 64 122 26 48 5 13 10 66 +SY01G 1820 31 36 37 30 55 5 51 64 66 25 +SY01G 1830 36 16 91 89 103 13 42 81 51 51 +SY01G 1840 50 32 64 55 70 62 58 17 91 96 +SY01G 1850 80 43 93 50 64 37 40 56 51 85 +SY01G 1860 65 31 41 85 37 53 34 50 56 59 +SY01G 1870 98 76 51 46 43 55 10 37 46 70 +SY01G 1880 48 64 98 95 59 51 89 59 24 52 +SY01G 1890 57 55 37 72 14 25 35 5 18 30 999 +SY01H 1690 261 253 184 302 171 5 224 91 103 5 +SY01H 1700 132 96 210 192 190 26 156 231 273 5 +SY01H 1710 194 211 208 134 151 126 140 175 189 149 +SY01H 1720 122 133 14 110 163 221 134 75 167 127 +SY01H 1730 199 113 82 27 43 41 71 154 173 207 +SY01H 1740 132 192 114 110 36 24 15 66 168 124 +SY01H 1750 162 137 161 99 43 5 67 119 140 126 +SY01H 1760 128 36 52 43 119 89 43 86 89 49 +SY01H 1770 6 17 81 84 51 91 133 136 79 5 +SY01H 1780 5 43 67 87 67 61 55 62 84 51 +SY01H 1790 59 53 64 98 52 59 82 81 67 76 +SY01H 1800 82 71 31 63 44 8 18 64 104 56 +SY01H 1810 69 64 63 118 26 57 5 33 13 64 +SY01H 1820 45 24 41 42 40 5 55 69 77 18 +SY01H 1830 24 15 90 105 108 21 37 88 45 54 +SY01H 1840 45 34 67 49 68 64 52 35 83 73 +SY01H 1850 80 48 78 45 50 42 44 50 47 72 +SY01H 1860 71 30 30 85 36 61 17 40 54 61 +SY01H 1870 81 67 37 40 60 37 5 23 34 61 +SY01H 1880 42 73 88 80 51 54 77 46 24 45 +SY01H 1890 41 55 41 60 27 20 32 5 16 38 999 diff --git a/tests/fixtures/reference/synth/synth02.rwl b/tests/fixtures/reference/synth/synth02.rwl new file mode 100644 index 0000000..be22bd9 --- /dev/null +++ b/tests/fixtures/reference/synth/synth02.rwl @@ -0,0 +1,162 @@ +SYNTH02 1 Synthetic Site SYNTH02 PIST +SYNTH02 2 Test Region PIST 45 00 N 72 00 W 300M +SYNTH02 3 Synthetic generator 1500 2000 +SY02A 1610 156 148 268 113 218 246 314 440 287 332 +SY02A 1620 63 180 272 342 245 260 234 317 229 62 +SY02A 1630 209 378 307 213 228 151 273 140 173 29 +SY02A 1640 87 217 146 36 84 138 196 134 165 134 +SY02A 1650 103 68 51 160 164 149 93 136 158 280 +SY02A 1660 63 31 122 218 126 218 120 55 110 94 +SY02A 1670 134 183 129 79 84 70 60 162 183 121 +SY02A 1680 158 37 107 167 112 109 115 102 78 119 +SY02A 1690 96 88 82 144 70 14 125 30 99 5 +SY02A 1700 62 46 115 106 99 28 77 124 125 5 +SY02A 1710 106 113 137 94 110 94 88 125 134 104 +SY02A 1720 83 87 29 77 99 160 104 30 127 111 +SY02A 1730 143 97 73 20 26 44 71 109 145 155 +SY02A 1740 103 146 93 76 9 24 5 71 128 106 +SY02A 1750 169 109 142 99 55 5 72 113 122 114 +SY02A 1760 120 28 55 37 125 87 45 75 102 38 +SY02A 1770 15 7 88 97 69 84 150 146 85 10 +SY02A 1780 5 56 62 91 63 83 55 64 83 71 +SY02A 1790 62 64 55 112 39 92 96 93 69 84 +SY02A 1800 70 66 31 73 65 17 23 106 111 71 +SY02A 1810 93 82 91 151 14 56 5 26 15 101 +SY02A 1820 54 50 56 45 57 5 65 104 89 13 +SY02A 1830 39 11 145 112 139 29 45 91 69 69 +SY02A 1840 61 50 75 81 85 96 86 37 122 131 +SY02A 1850 106 70 124 54 92 39 59 67 48 101 +SY02A 1860 90 46 66 118 77 104 21 61 84 88 +SY02A 1870 125 114 56 68 77 69 25 49 32 86 +SY02A 1880 70 96 138 109 79 70 111 62 36 49 +SY02A 1890 83 72 67 99 28 13 55 5 48 44 999 +SY02B 1630 268 429 340 270 253 116 346 196 207 37 +SY02B 1640 121 265 159 73 159 204 199 177 214 170 +SY02B 1650 81 84 81 223 216 203 163 176 190 298 +SY02B 1660 5 24 155 316 165 318 160 73 157 183 +SY02B 1670 169 266 192 84 80 72 82 181 250 125 +SY02B 1680 216 61 121 190 144 125 144 129 96 122 +SY02B 1690 135 134 81 155 119 28 142 69 99 5 +SY02B 1700 61 74 124 117 123 28 94 146 141 5 +SY02B 1710 110 154 149 113 135 88 75 137 161 114 +SY02B 1720 100 100 11 87 133 151 95 51 122 118 +SY02B 1730 166 84 94 5 18 29 61 125 148 175 +SY02B 1740 131 160 91 91 21 28 5 55 122 97 +SY02B 1750 173 114 164 86 45 5 75 121 112 121 +SY02B 1760 139 32 60 56 123 92 46 77 85 50 +SY02B 1770 5 27 90 75 65 94 150 137 97 5 +SY02B 1780 5 45 53 97 55 69 56 77 80 69 +SY02B 1790 65 49 68 123 60 79 96 84 77 66 +SY02B 1800 81 83 29 69 45 6 13 85 118 77 +SY02B 1810 67 86 84 131 28 64 5 19 18 78 +SY02B 1820 57 32 59 50 60 5 57 78 81 13 +SY02B 1830 42 20 120 113 133 19 45 86 69 59 +SY02B 1840 55 51 74 68 79 85 63 51 106 103 +SY02B 1850 99 48 103 62 79 35 40 62 56 97 +SY02B 1860 82 22 49 108 65 72 27 45 61 70 +SY02B 1870 108 87 53 46 57 49 13 41 37 88 +SY02B 1880 68 70 108 115 66 59 104 56 24 51 +SY02B 1890 66 58 71 89 33 21 52 5 28 52 999 +SY02C 1640 119 243 103 99 100 123 204 138 150 153 +SY02C 1650 119 28 54 185 186 135 114 120 145 269 +SY02C 1660 20 52 145 264 153 219 139 65 110 83 +SY02C 1670 128 168 150 68 90 66 73 135 196 105 +SY02C 1680 149 52 82 162 109 98 113 72 35 109 +SY02C 1690 95 93 59 118 56 18 96 42 65 5 +SY02C 1700 42 37 92 75 63 23 58 77 104 5 +SY02C 1710 81 96 88 70 81 66 56 91 84 63 +SY02C 1720 63 63 13 51 83 111 71 31 74 74 +SY02C 1730 90 57 54 16 27 38 37 81 106 109 +SY02C 1740 67 116 65 59 13 13 5 53 85 77 +SY02C 1750 111 78 101 71 27 5 42 76 82 87 +SY02C 1760 82 27 38 28 88 63 21 61 59 39 +SY02C 1770 18 14 71 51 37 67 97 104 71 9 +SY02C 1780 5 29 44 53 31 50 41 44 49 39 +SY02C 1790 47 41 44 71 33 50 66 60 51 45 +SY02C 1800 53 57 27 50 41 6 8 63 71 56 +SY02C 1810 57 63 64 92 27 46 5 16 18 57 +SY02C 1820 34 32 39 32 50 5 54 61 63 9 +SY02C 1830 18 5 78 78 86 15 27 72 55 36 +SY02C 1840 36 26 52 47 48 51 61 41 77 76 +SY02C 1850 79 43 86 51 59 29 35 43 39 68 +SY02C 1860 60 21 32 72 39 57 15 36 47 52 +SY02C 1870 72 71 22 38 47 38 16 24 26 62 +SY02C 1880 52 57 83 76 50 48 83 44 26 29 +SY02C 1890 47 49 42 66 16 17 36 5 29 27 999 +SY02D 1600 273 5 478 292 246 270 151 187 173 169 +SY02D 1610 174 122 230 157 169 214 293 405 278 249 +SY02D 1620 108 188 324 259 247 269 236 300 240 103 +SY02D 1630 211 370 291 240 196 99 266 118 156 5 +SY02D 1640 140 209 171 51 120 146 191 147 173 102 +SY02D 1650 110 87 78 192 201 151 115 124 215 301 +SY02D 1660 25 92 154 284 152 294 162 82 142 146 +SY02D 1670 151 254 191 74 142 90 99 203 273 107 +SY02D 1680 212 106 155 225 172 122 156 111 74 203 +SY02D 1690 151 133 73 193 100 36 158 74 115 5 +SY02D 1700 82 73 173 126 148 43 140 167 175 5 +SY02D 1710 135 184 183 130 156 129 98 181 180 135 +SY02D 1720 116 108 18 129 176 224 143 47 132 139 +SY02D 1730 202 112 84 19 72 67 80 178 223 234 +SY02D 1740 171 222 148 120 5 33 5 109 190 154 +SY02D 1750 238 157 217 116 74 7 72 158 183 158 +SY02D 1760 199 53 59 54 181 115 34 136 137 79 +SY02D 1770 30 44 136 153 74 132 211 216 135 5 +SY02D 1780 5 58 89 124 89 104 92 83 140 83 +SY02D 1790 97 52 69 151 68 91 118 124 112 115 +SY02D 1800 111 111 38 105 93 17 26 136 160 112 +SY02D 1810 125 127 108 195 36 94 5 42 16 120 +SY02D 1820 85 58 57 82 79 20 97 118 137 28 +SY02D 1830 65 15 166 154 182 40 66 124 101 91 +SY02D 1840 82 50 113 99 102 107 107 53 143 151 +SY02D 1850 139 91 151 85 115 59 73 58 91 122 +SY02D 1860 108 46 63 164 63 109 33 62 97 92 +SY02D 1870 148 126 67 75 109 87 27 81 69 96 +SY02D 1880 77 122 151 153 96 110 135 90 28 58 +SY02D 1890 106 87 60 97 35 11 35 5 30 55 999 +SY02E 1680 418 130 216 427 307 246 305 302 103 293 +SY02E 1690 246 241 180 376 196 31 254 127 188 5 +SY02E 1700 162 106 228 241 205 90 199 203 268 5 +SY02E 1710 185 261 262 187 271 198 126 245 292 217 +SY02E 1720 137 160 42 167 243 293 163 104 216 199 +SY02E 1730 264 179 133 29 72 33 100 208 256 284 +SY02E 1740 192 255 134 141 17 75 14 116 203 185 +SY02E 1750 225 190 231 148 72 5 101 165 188 171 +SY02E 1760 205 68 68 49 182 118 56 100 128 49 +SY02E 1770 31 19 132 117 101 117 181 199 120 13 +SY02E 1780 5 69 88 119 71 102 102 102 118 73 +SY02E 1790 84 71 82 156 53 115 116 84 106 95 +SY02E 1800 87 103 45 83 63 5 25 125 116 99 +SY02E 1810 95 107 108 165 43 58 5 16 5 94 +SY02E 1820 52 60 53 59 81 5 64 82 93 12 +SY02E 1830 30 25 129 122 138 14 45 107 51 56 +SY02E 1840 72 44 73 71 72 81 87 48 110 114 +SY02E 1850 118 90 117 66 76 56 56 55 49 87 +SY02E 1860 81 38 50 105 39 79 32 42 67 75 +SY02E 1870 100 100 57 53 68 60 17 43 54 82 +SY02E 1880 71 65 116 104 59 77 113 81 35 37 +SY02E 1890 70 70 44 86 23 29 33 5 35 50 999 +SY02F 1655 202 161 193 237 350 +SY02F 1660 61 88 173 368 197 376 260 46 174 128 +SY02F 1670 176 267 206 121 94 38 127 233 266 152 +SY02F 1680 237 81 132 228 176 143 167 148 77 176 +SY02F 1690 159 149 105 175 135 18 167 58 93 5 +SY02F 1700 90 88 118 116 138 35 123 136 153 5 +SY02F 1710 125 149 157 126 145 114 86 137 147 124 +SY02F 1720 101 95 6 85 133 176 113 40 137 118 +SY02F 1730 120 101 85 5 49 44 73 132 161 171 +SY02F 1740 123 150 102 100 23 26 5 67 124 114 +SY02F 1750 165 118 141 90 39 5 50 87 110 136 +SY02F 1760 121 42 56 48 133 87 55 81 94 47 +SY02F 1770 15 14 86 77 53 85 124 123 104 5 +SY02F 1780 5 46 61 88 55 67 63 67 92 44 +SY02F 1790 63 35 53 103 43 55 97 82 62 87 +SY02F 1800 76 74 33 69 61 8 19 88 100 72 +SY02F 1810 75 80 72 131 27 73 5 26 14 72 +SY02F 1820 41 48 49 38 59 5 63 82 78 31 +SY02F 1830 30 8 117 112 119 16 37 79 53 53 +SY02F 1840 50 42 77 65 67 73 58 40 83 88 +SY02F 1850 96 68 107 52 74 35 52 51 46 85 +SY02F 1860 69 32 43 102 45 77 33 44 55 71 +SY02F 1870 112 81 33 54 56 56 28 40 37 67 +SY02F 1880 58 65 107 93 66 69 96 60 21 38 +SY02F 1890 70 68 29 69 22 16 38 5 28 34 999 diff --git a/tests/fixtures/samples/known_1789.csv b/tests/fixtures/samples/known_1789.csv new file mode 100644 index 0000000..f9d599f --- /dev/null +++ b/tests/fixtures/samples/known_1789.csv @@ -0,0 +1,91 @@ +year,width_mm +1700,2.003 +1701,1.949 +1702,3.574 +1703,3.249 +1704,3.109 +1705,1.394 +1706,2.602 +1707,4.637 +1708,4.467 +1709,0.050 +1710,3.439 +1711,4.712 +1712,4.522 +1713,3.002 +1714,3.690 +1715,3.071 +1716,2.131 +1717,4.384 +1718,4.148 +1719,3.407 +1720,2.739 +1721,2.061 +1722,0.141 +1723,2.621 +1724,2.987 +1725,4.493 +1726,3.154 +1727,1.266 +1728,3.417 +1729,2.806 +1730,4.278 +1731,2.732 +1732,2.377 +1733,0.050 +1734,1.238 +1735,0.800 +1736,1.608 +1737,3.942 +1738,3.905 +1739,4.672 +1740,3.170 +1741,4.092 +1742,2.784 +1743,2.721 +1744,0.678 +1745,0.770 +1746,0.219 +1747,1.933 +1748,3.037 +1749,2.992 +1750,3.722 +1751,3.143 +1752,4.003 +1753,2.651 +1754,1.352 +1755,0.186 +1756,1.544 +1757,2.580 +1758,3.115 +1759,3.501 +1760,3.212 +1761,0.781 +1762,1.387 +1763,1.315 +1764,2.965 +1765,2.142 +1766,0.690 +1767,1.989 +1768,2.059 +1769,1.624 +1770,0.748 +1771,0.470 +1772,2.078 +1773,2.221 +1774,1.217 +1775,2.152 +1776,3.162 +1777,3.510 +1778,2.131 +1779,0.282 +1780,0.050 +1781,1.363 +1782,1.536 +1783,1.997 +1784,1.461 +1785,1.896 +1786,1.382 +1787,1.283 +1788,1.754 +1789,1.197 diff --git a/tests/test_api.py b/tests/test_api.py new file mode 100644 index 0000000..c9d1f6d --- /dev/null +++ b/tests/test_api.py @@ -0,0 +1,90 @@ +"""Tests for the FastAPI dating service.""" + +import io +from pathlib import Path + +import pandas as pd +import pytest + +pytest.importorskip("fastapi") +from fastapi.testclient import TestClient # noqa: E402 + +from dendro.api import app as api_app # noqa: E402 +from dendro.crossdating.matcher import CrossdateMatcher # noqa: E402 +from dendro.reference.chronology_index import ChronologyIndex # noqa: E402 + +FIX = Path(__file__).parent / "fixtures" +SAMPLE = FIX / "samples" / "known_1789.csv" + + +@pytest.fixture +def client(monkeypatch): + # Point the service at the committed synthetic reference fixtures. + matcher = CrossdateMatcher(index=ChronologyIndex(FIX / "reference")) + api_app.state.matcher = matcher + api_app.state.reference_dir = FIX / "reference" + with TestClient(api_app.app) as c: + # TestClient lifespan reloads the matcher from DENDRO_REFERENCE_DIR; + # re-pin it to the fixtures afterwards. + api_app.state.matcher = matcher + yield c + + +def test_health(client): + r = client.get("/health") + assert r.status_code == 200 + assert r.json()["status"] == "ok" + + +def test_references(client): + r = client.get("/references") + assert r.status_code == 200 + body = r.json() + assert body["count"] == 2 + # Synthetic fixture filenames don't encode species/state (as with many real + # ITRDB files), so just assert the coverage envelope is reported. + assert body["coverage"]["start"] is not None + assert body["coverage"]["end"] is not None + + +def test_date_json_recovers_year(client): + widths = pd.read_csv(SAMPLE)["width_mm"].tolist() + r = client.post( + "/date", + json={ + "widths": widths, + "has_bark_edge": True, + "era_start": 1700, + "era_end": 1850, + "min_overlap": 40, + }, + ) + assert r.status_code == 200 + body = r.json() + assert body["consensus_year"] == 1789 + assert body["felling_estimate"]["felling_type"] == "exact" + + +def test_date_csv_upload(client): + raw = SAMPLE.read_bytes() + r = client.post( + "/date/csv?era_start=1700&era_end=1850", + files={"file": ("known_1789.csv", io.BytesIO(raw), "text/csv")}, + ) + assert r.status_code == 200 + assert r.json()["consensus_year"] == 1789 + + +def test_date_rejects_negative_widths(client): + r = client.post("/date", json={"widths": [1.0, -2.0, 3.0] * 20}) + assert r.status_code == 422 + + +def test_parse_upload(client): + raw = (FIX / "reference" / "synth" / "synth01.rwl").read_bytes() + r = client.post( + "/parse", + files={"file": ("synth01.rwl", io.BytesIO(raw), "text/plain")}, + ) + assert r.status_code == 200 + assert r.json()["series_count"] == 8 diff --git a/tests/test_correlator.py b/tests/test_correlator.py index e14406e..58f9a7f 100644 --- a/tests/test_correlator.py +++ b/tests/test_correlator.py @@ -1,15 +1,14 @@ """Tests for cross-dating correlation algorithms.""" import numpy as np -import pytest from dendro.crossdating.correlator import ( - sliding_correlation, - calculate_tvalue, + CorrelationResult, calculate_gleichlauf, - find_best_match, + calculate_tvalue, dating_confidence, - CorrelationResult, + find_best_match, + sliding_correlation, ) diff --git a/tests/test_detrend.py b/tests/test_detrend.py index 0a494dd..94ae085 100644 --- a/tests/test_detrend.py +++ b/tests/test_detrend.py @@ -4,11 +4,11 @@ import pytest from dendro.crossdating.detrend import ( + DetrendMethod, + build_chronology, detrend_series, - standardize, prewhiten, - build_chronology, - DetrendMethod, + standardize, ) diff --git a/tests/test_felling.py b/tests/test_felling.py new file mode 100644 index 0000000..78f6e00 --- /dev/null +++ b/tests/test_felling.py @@ -0,0 +1,44 @@ +"""Tests for felling-year interpretation.""" + +from dendro.crossdating.felling import FellingType, estimate_felling + + +def test_bark_edge_exact(): + fe = estimate_felling(1789, has_bark_edge=True) + assert fe.felling_type is FellingType.EXACT + assert fe.earliest_felling == 1789 + assert fe.latest_felling == 1789 + + +def test_no_bark_no_sapwood_is_after(): + fe = estimate_felling(1789, has_bark_edge=False) + assert fe.felling_type is FellingType.AFTER + assert fe.earliest_felling > 1789 + assert fe.latest_felling is None + + +def test_sapwood_present_gives_range(): + fe = estimate_felling( + 1789, has_bark_edge=False, has_sapwood=True, sapwood_count=5, species="QUAL" + ) + assert fe.felling_type is FellingType.RANGE + assert fe.earliest_felling >= 1789 + assert fe.latest_felling is not None + assert fe.latest_felling >= fe.earliest_felling + + +def test_oak_uses_oak_sapwood_model(): + oak = estimate_felling(1800, has_bark_edge=False, has_sapwood=True, species="QURU") + other = estimate_felling(1800, has_bark_edge=False, has_sapwood=True, species="PIST") + # The two species groups should not produce identical ranges in general. + assert (oak.earliest_felling, oak.latest_felling) != ( + other.earliest_felling, other.latest_felling + ) or oak.note != other.note + + +def test_summary_and_dict_roundtrip(): + fe = estimate_felling(1789, has_bark_edge=True) + assert "1789" in fe.summary() + d = fe.to_dict() + assert d["felling_type"] == "exact" + assert d["earliest_felling"] == 1789 diff --git a/tests/test_fixtures.py b/tests/test_fixtures.py new file mode 100644 index 0000000..5eb7674 --- /dev/null +++ b/tests/test_fixtures.py @@ -0,0 +1,34 @@ +"""Ensure the committed synthetic fixtures are reproducible and parseable.""" + +from pathlib import Path + +from dendro.reference.tucson_parser import parse_rwl_file + +FIX = Path(__file__).parent / "fixtures" +REF = FIX / "reference" / "synth" + + +def test_fixture_files_exist(): + assert (REF / "synth01.rwl").exists() + assert (REF / "synth02.rwl").exists() + assert (FIX / "samples" / "known_1789.csv").exists() + + +def test_fixtures_parse(): + rwl = parse_rwl_file(REF / "synth01.rwl") + assert len(rwl.series) == 8 + for series in rwl.series.values(): + assert series.length > 100 + assert (series.values > 0).all() + + +def test_fixtures_are_deterministic(tmp_path, monkeypatch): + # Regenerating into a temp dir must reproduce the committed bytes exactly. + import tests.fixtures.generate_fixtures as gen + + committed = (REF / "synth01.rwl").read_bytes() + monkeypatch.setattr(gen, "REF_DIR", tmp_path / "ref") + monkeypatch.setattr(gen, "SAMPLE_DIR", tmp_path / "samples") + gen.generate_all() + regenerated = (tmp_path / "ref" / "synth01.rwl").read_bytes() + assert regenerated == committed diff --git a/tests/test_known_date.py b/tests/test_known_date.py index 5683669..ff34c43 100644 --- a/tests/test_known_date.py +++ b/tests/test_known_date.py @@ -6,7 +6,6 @@ """ import numpy as np -import pytest from dendro.crossdating.correlator import find_best_match from dendro.crossdating.detrend import standardize @@ -22,7 +21,6 @@ def test_exact_subset_recovery(self): # Create a "reference chronology" with distinctive pattern # Mix of periodic signal and random component n_years = 200 - years = np.arange(1700, 1700 + n_years) # Climate-like signal: multi-frequency sine waves signal = ( diff --git a/tests/test_science_fixes.py b/tests/test_science_fixes.py new file mode 100644 index 0000000..40ea5f3 --- /dev/null +++ b/tests/test_science_fixes.py @@ -0,0 +1,98 @@ +""" +Targeted regression tests for the cross-dating correctness fixes: +- Cook & Peters spline frequency response +- Gleichlaeufigkeit significance +- missing/false ring detection +- Tucson writer/parser round-trip +""" + +import tempfile +from pathlib import Path + +import numpy as np + +from dendro.crossdating.correlator import ( + calculate_gleichlauf, + detect_missing_ring, + gleichlauf_significance, +) +from dendro.crossdating.detrend import _fit_cubic_spline, standardize +from dendro.imaging.path_sampler import widths_to_tucson +from dendro.reference.tucson_parser import parse_rwl_file + + +class TestSplineFrequencyResponse: + def _response(self, period, wavelength, n=600): + t = np.arange(n) + y = 10.0 + np.sin(2 * np.pi * t / wavelength) + curve = _fit_cubic_spline(t, y, np.ones(n, bool), period) - 10.0 + k = int(wavelength) + sl = slice(k, n - k) + s = np.sin(2 * np.pi * t / wavelength) + c = np.cos(2 * np.pi * t / wavelength) + a = 2 * np.mean(curve[sl] * s[sl]) + b = 2 * np.mean(curve[sl] * c[sl]) + return np.hypot(a, b) + + def test_fifty_percent_at_cutoff(self): + # By definition the response is 0.5 at the cutoff wavelength. + assert abs(self._response(50.0, 50.0) - 0.5) < 0.03 + + def test_long_wavelengths_retained(self): + assert self._response(50.0, 200.0) > 0.9 + + def test_short_wavelengths_removed(self): + assert self._response(50.0, 12.5) < 0.1 + + +class TestGleichlaufSignificance: + def test_chance_level_not_significant(self): + # 50% GLK over a long series is not significant. + assert gleichlauf_significance(50.0, 100) > 0.5 + + def test_high_glk_significant(self): + assert gleichlauf_significance(80.0, 100) < 0.01 + + def test_partial_agreement_value(self): + s1 = np.array([1, 2, 3, 2, 3]) + s2 = np.array([1, 2, 1, 2, 3]) + assert 40 <= calculate_gleichlauf(s1, s2) <= 60 + + +class TestMissingRing: + def test_locates_deleted_ring(self): + rng = np.random.default_rng(7) + ref = standardize(rng.normal(0, 1, 200) + np.sin(np.linspace(0, 20 * np.pi, 200))) + true = ref[40:141].copy() + sample = np.delete(true, 50)[:100] + hint = detect_missing_ring(sample, ref[40:160]) + assert hint is not None + assert hint.kind == "missing" + assert abs(hint.index - 50) <= 2 + assert hint.gain > 0.1 + + def test_clean_sample_no_hint(self): + rng = np.random.default_rng(8) + ref = standardize(rng.normal(0, 1, 200) + np.sin(np.linspace(0, 20 * np.pi, 200))) + clean = ref[40:140].copy() + assert detect_missing_ring(clean, ref[40:160]) is None + + +class TestTucsonRoundTrip: + def test_round_trip_widths(self): + rng = np.random.default_rng(1) + widths_mm = np.round(rng.uniform(0.5, 3.0, 60), 2) + text = widths_to_tucson(widths_mm, "RT01", end_year=1850) + with tempfile.NamedTemporaryFile("w", suffix=".rwl", delete=False) as f: + f.write(text) + path = f.name + try: + rwl = parse_rwl_file(path) + series = next(iter(rwl.series.values())) + got = series.values / 100.0 # 0.01 mm -> mm + expect = widths_mm[::-1] # writer stores oldest-first + assert series.end_year == 1850 + assert len(got) == len(expect) + assert np.allclose(got, expect, atol=1e-6) + finally: + Path(path).unlink() diff --git a/tests/test_synth_dating.py b/tests/test_synth_dating.py new file mode 100644 index 0000000..e119bfc --- /dev/null +++ b/tests/test_synth_dating.py @@ -0,0 +1,90 @@ +""" +End-to-end dating validation against committed synthetic fixtures. + +Unlike ``test_validation.py`` (which needs downloaded ITRDB data and therefore +skips in CI), these tests run against the deterministic fixtures under +``tests/fixtures`` and exercise the full pipeline: parse -> per-series detrend -> +robust master chronology -> sliding correlation -> felling interpretation. + +The fixtures are built so the known felling year is 1789; recovering it proves +the master-chronology, orientation, and felling fixes work together. +""" + +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from dendro.crossdating.matcher import CrossdateMatcher +from dendro.reference.chronology_index import ChronologyIndex + +FIX = Path(__file__).parent / "fixtures" +REF = FIX / "reference" +SAMPLE = FIX / "samples" / "known_1789.csv" + + +@pytest.fixture(scope="module") +def matcher() -> CrossdateMatcher: + return CrossdateMatcher(index=ChronologyIndex(REF)) + + +@pytest.fixture(scope="module") +def sample_widths() -> np.ndarray: + return pd.read_csv(SAMPLE)["width_mm"].values + + +def test_recovers_known_felling_year(matcher, sample_widths): + report = matcher.date_sample( + values=sample_widths, + sample_name="known_1789", + has_bark_edge=True, + era_start=1700, + era_end=1850, + min_overlap=40, + ) + assert report.consensus_year == 1789 + assert report.consensus_confidence == "HIGH" + # Strong, unambiguous match expected from clean synthetic data. + assert report.matches[0].correlation > 0.8 + assert report.matches[0].t_value > 10 + + +def test_bark_edge_gives_exact_felling(matcher, sample_widths): + report = matcher.date_sample( + values=sample_widths, has_bark_edge=True, + era_start=1700, era_end=1850, min_overlap=40, + ) + fe = report.felling_estimate + assert fe is not None + assert fe.felling_type.value == "exact" + assert fe.earliest_felling == fe.latest_felling == 1789 + + +def test_no_bark_edge_is_terminus_post_quem(matcher, sample_widths): + report = matcher.date_sample( + values=sample_widths, has_bark_edge=False, + era_start=1700, era_end=1850, min_overlap=40, + ) + fe = report.felling_estimate + assert fe is not None + assert fe.felling_type.value == "after" + # Felled strictly after the last dated ring. + assert fe.earliest_felling > 1789 + assert fe.latest_felling is None + + +def test_orientation_reversed_still_dates(matcher, sample_widths): + # The same series measured bark-first must recover the same year. + report = matcher.date_sample( + values=sample_widths[::-1], has_bark_edge=True, + orientation="bark_to_pith", + era_start=1700, era_end=1850, min_overlap=40, + ) + assert report.consensus_year == 1789 + + +def test_master_chronology_is_cached(matcher, sample_widths): + matcher.date_sample(values=sample_widths, era_start=1700, era_end=1850, min_overlap=40) + # Every indexed reference should now have a cache entry. + assert len(matcher._master_cache) >= 1 diff --git a/tests/test_tucson_parser.py b/tests/test_tucson_parser.py index f14b1c5..3684ad5 100644 --- a/tests/test_tucson_parser.py +++ b/tests/test_tucson_parser.py @@ -1,16 +1,15 @@ """Tests for Tucson format parser.""" -import numpy as np -import pytest -from io import StringIO -from pathlib import Path import tempfile +from pathlib import Path + +import numpy as np from dendro.reference.tucson_parser import ( - parse_rwl_file, - _parse_rwl_line, RingWidthSeries, RWLFile, + _parse_rwl_line, + parse_rwl_file, ) diff --git a/tests/test_validation.py b/tests/test_validation.py index 5a27cd6..0c368ec 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -5,10 +5,11 @@ known dates from actual tree ring data. """ -import pytest -import numpy as np from pathlib import Path +import numpy as np +import pytest + # Skip if reference data not downloaded REF_DIR = Path(__file__).parent.parent / "data" / "reference" SKIP_IF_NO_DATA = pytest.mark.skipif( @@ -28,7 +29,7 @@ class TestRealDataValidation: def _build_site_master(self, rwl, exclude_series_id): """Build master chronology from all series except one.""" - from dendro.crossdating.detrend import detrend_series, standardize, DetrendMethod + from dendro.crossdating.detrend import DetrendMethod, detrend_series, standardize other_series = [(sid, s) for sid, s in rwl.series.items() if sid != exclude_series_id and s.length >= 50] @@ -56,9 +57,9 @@ def _build_site_master(self, rwl, exclude_series_id): def _crossdate_within_site(self, rwl_path): """Test cross-dating one core against others from same site.""" - from dendro.reference.tucson_parser import parse_rwl_file from dendro.crossdating.correlator import find_best_match - from dendro.crossdating.detrend import detrend_series, standardize, DetrendMethod + from dendro.crossdating.detrend import DetrendMethod, detrend_series, standardize + from dendro.reference.tucson_parser import parse_rwl_file rwl = parse_rwl_file(rwl_path) diff --git a/tests/test_visualization.py b/tests/test_visualization.py index d9ba6b6..011750f 100644 --- a/tests/test_visualization.py +++ b/tests/test_visualization.py @@ -3,12 +3,11 @@ """ import numpy as np -import pytest from dendro.visualization.plots import ( + MARKER_YEARS, detect_marker_years, identify_known_markers, - MARKER_YEARS, ) @@ -93,10 +92,10 @@ class TestPlotFunctions: def test_import_plot_functions(self): """All plot functions should be importable.""" from dendro.visualization.plots import ( - plot_crossdate_results, plot_correlation_profile, - plot_segment_analysis, + plot_crossdate_results, plot_sample_reference_overlay, + plot_segment_analysis, save_diagnostic_plots, )