"""Band-ratio weathering indicators (carbonyl index) and composition tables.
Two important caveats, both demonstrated in the study:
1. `spectral_moments` normalization centers intensities around zero, so raw
integrated band areas can be negative and a naive band *ratio* is unstable.
`carbonyl_index` therefore per-spectrum min-max renormalizes to a positive
[0, 1] scale before integrating, and floors the denominator.
2. The carbonyl index conflates weathering with intrinsic polymer chemistry
(pristine PET/PMMA have strong native C=O bands), so it is only a valid
oxidation marker *within a single polyolefin class* (PE/PP). On this
heterogeneous, multi-instrument dataset it does not cleanly separate
environmental from pristine spectra — which is itself a finding: hand-crafted
physical indices are unreliable across instruments, motivating data-driven
OOD detection instead.
"""
from __future__ import annotations
import numpy as np
import pandas as pd
from scipy.stats import spearmanr
_trapz = getattr(np, "trapezoid", None) or np.trapz
[docs]
def band_area(
X: np.ndarray,
wavenumbers: np.ndarray,
lo: float,
hi: float,
mode: str = "area",
subtract_local_baseline: bool = True,
) -> np.ndarray:
"""Integrated area (or peak height) of a wavenumber window.
With subtract_local_baseline, a straight line between the window's
endpoints is removed first — the usual convention for band indices.
"""
idx = np.flatnonzero((wavenumbers >= lo) & (wavenumbers <= hi))
if idx.size < 3:
raise ValueError(f"Band [{lo}, {hi}] covers fewer than 3 grid points.")
wn = wavenumbers[idx]
band = X[:, idx].astype(float)
if subtract_local_baseline:
left, right = band[:, [0]], band[:, [-1]]
frac = ((wn - wn[0]) / (wn[-1] - wn[0]))[None, :]
band = band - (left + (right - left) * frac)
if mode == "height":
return band.max(axis=1)
return _trapz(band, wn, axis=1)
[docs]
def carbonyl_index(
X: np.ndarray,
wavenumbers: np.ndarray,
carbonyl: tuple[float, float] = (1650.0, 1800.0),
ref: tuple[float, float] = (1420.0, 1470.0),
mode: str = "area",
renormalize: bool = True,
) -> np.ndarray:
"""C=O band relative to the CH2-scissoring reference band.
The default bands match ``TrustConfig.carbonyl_band`` / ``reference_band``
so every call site uses one band definition unless overridden explicitly.
With renormalize=True (default), each spectrum is min-max scaled to [0, 1]
first so the band areas are non-negative and the ratio is stable on
moments-normalized data. See module docstring for the within-class caveat.
"""
Xc = X.astype(float)
if renormalize:
lo = Xc.min(axis=1, keepdims=True)
hi = Xc.max(axis=1, keepdims=True)
Xc = (Xc - lo) / np.clip(hi - lo, 1e-9, None)
num = band_area(Xc, wavenumbers, *carbonyl, mode=mode, subtract_local_baseline=False)
den = band_area(Xc, wavenumbers, *ref, mode=mode, subtract_local_baseline=False)
return num / np.clip(den, 1e-6, None)
[docs]
def weathering_report(
ci: np.ndarray,
confidence: np.ndarray,
ood_scores: dict[str, np.ndarray],
) -> pd.DataFrame:
"""Spearman correlations of carbonyl index vs confidence and OOD scores."""
rows = []
valid = np.isfinite(ci)
rho, p = spearmanr(ci[valid], confidence[valid])
rows.append({"pair": "carbonyl_index_vs_confidence", "spearman_rho": rho,
"p_value": p, "n": int(valid.sum())})
for name, scores in ood_scores.items():
rho, p = spearmanr(ci[valid], scores[valid])
rows.append({"pair": f"carbonyl_index_vs_{name}", "spearman_rho": rho,
"p_value": p, "n": int(valid.sum())})
return pd.DataFrame(rows)
[docs]
def unknown_composition(
labels: np.ndarray,
class_names: list[str],
abstain_label: str = "ABSTAIN",
) -> pd.DataFrame:
"""Class fractions among the unknowns, with and without abstentions."""
labels = np.asarray(labels, dtype=object)
n = len(labels)
kept = labels != abstain_label
rows = []
for name in list(class_names) + [abstain_label]:
count = int(np.sum(labels == name))
rows.append({
"class": name,
"count": count,
"fraction_of_all": count / n,
"fraction_of_classified": (count / kept.sum()) if (name != abstain_label and kept.any()) else np.nan,
})
return pd.DataFrame(rows)