API Reference#

Quick reference#

Function

Description

compute_scaledRI

Compute Scaled RI (RI₀₋₃ by default) from six isoGDGT abundances

predict_proxy_from_T

Forward: temperature → proxy percentiles (pure Python)

predict_T_from_proxyObs

Inverse: proxy → temperature with full uncertainty (runs Stan)

download_posteriors

Download forward posteriors from Zenodo

download_training_data

Download training CSVs + CMEMS NO₃ field

list_posteriors

Print and return .nc stems in the local cache

build_fwd_data

Build validated Stan data dict for forward calibration

get_posterior

Run forward calibration Stan sampling

save_posterior

Persist forward posterior as compressed NetCDF

load_posterior

Load a forward or invT posterior from the cache

summarize_sampler_diagnostics

Divergences, R-hat, ESS, E-BFMI


Prediction#

Compute Scaled Ring Index#

compute_scaledRI(gdgt0, gdgt1, gdgt2, gdgt3, cren, cren_prime, *, cren_rings=3)[source]#

Compute Scaled Ring Index from six isoGDGT abundances.

Accepts raw LC/MS peak areas or fractional abundances — both give identical results because the formula divides by the total sum of all six GDGTs, so any common scale factor drops out. Default cren_rings=3 produces scaledRI_cren3 (RI₀₋₃), the canonical proxy used in TEXAS calibration posteriors.

Parameters:
  • gdgt0 (float or array-like) – isoGDGT abundances — GDGT-0, GDGT-1, GDGT-2, GDGT-3, crenarchaeol, crenarchaeol regioisomer (cren’). Raw LC/MS peak areas and fractional abundances give the same result (see above).

  • gdgt1 (float or array-like) – isoGDGT abundances — GDGT-0, GDGT-1, GDGT-2, GDGT-3, crenarchaeol, crenarchaeol regioisomer (cren’). Raw LC/MS peak areas and fractional abundances give the same result (see above).

  • gdgt2 (float or array-like) – isoGDGT abundances — GDGT-0, GDGT-1, GDGT-2, GDGT-3, crenarchaeol, crenarchaeol regioisomer (cren’). Raw LC/MS peak areas and fractional abundances give the same result (see above).

  • gdgt3 (float or array-like) – isoGDGT abundances — GDGT-0, GDGT-1, GDGT-2, GDGT-3, crenarchaeol, crenarchaeol regioisomer (cren’). Raw LC/MS peak areas and fractional abundances give the same result (see above).

  • cren (float or array-like) – isoGDGT abundances — GDGT-0, GDGT-1, GDGT-2, GDGT-3, crenarchaeol, crenarchaeol regioisomer (cren’). Raw LC/MS peak areas and fractional abundances give the same result (see above).

  • cren_prime (float or array-like) – isoGDGT abundances — GDGT-0, GDGT-1, GDGT-2, GDGT-3, crenarchaeol, crenarchaeol regioisomer (cren’). Raw LC/MS peak areas and fractional abundances give the same result (see above).

  • cren_rings (int) – Ring count assigned to both crenarchaeol and its regioisomer. 3 → scaledRI_cren3 / RI₀₋₃ (default, recommended). 4 → scaledRI / RI₀₋₄ (Zhang et al. 2016 convention).

Returns:

Scaled Ring Index, dimensionless, nominally in [0, 1].

Return type:

numpy.ndarray or float

Notes

The formula is:

RI      = (1·GDGT1 + 2·GDGT2 + 3·GDGT3 + cren_rings·cren + cren_rings·cren')
          / (GDGT0 + GDGT1 + GDGT2 + GDGT3 + cren + cren')
scaledRI = RI / cren_rings

Examples

>>> compute_scaledRI(0.45, 0.10, 0.08, 0.05, 0.30, 0.02)
array(0.45666...)
>>> import pandas as pd
>>> df = pd.read_csv("my_gdgt_data.csv")
>>> df["scaledRI_cren3"] = compute_scaledRI(
...     df["GDGT-0"], df["GDGT-1"], df["GDGT-2"], df["GDGT-3"],
...     df["cren"],   df["cren_prime"],
... )

Predict proxy from T#

predict_proxy_from_T(temperatures, posterior, *, n_draws=500, percentiles=[5, 50, 95], return_full=False, seed=42, gdgt23ratio=None, no3=None, no3_cutoff=None, suffix=None)[source]#

Forward prediction: temperature → proxy percentiles (Scaled RI, TEX86, or any fitted proxy).

Samples n_draws self-consistent parameter sets from the forward calibration posterior (all parameters drawn from the same posterior index, preserving correlations) and evaluates the calibration curve at each requested temperature. Corresponds to the forward model described in Eq. 1 / Eq. 6–7 of the manuscript.

Parameters:
  • temperatures (array-like) – Temperatures (°C) at which to evaluate the calibration curve.

  • posterior (xr.Dataset or str) – Forward calibration posterior — either a loaded xr.Dataset or a saved-file name string (looked up in the posterior cache).

  • n_draws (int) – Number of posterior draws to sample. Default 500.

  • percentiles (list of float) – Percentiles to return, e.g. [5, 50, 95].

  • return_full (bool) – If True, also return the full (n_draws × len(temperatures)) ensemble array and run metadata under keys "ensemble" and "metadata".

  • seed (int) – Random seed for reproducible draw sampling.

  • gdgt23ratio (array-like, optional) – GDGT-2/GDGT-3 ratio values (one per temperature point). Required only when the posterior was fitted with the multivariate model (β_{G₂/₃} correction).

  • no3 (array-like, optional) – Nitrate concentration values (one per temperature point). Required only when the posterior was fitted with NO₃ correction.

  • no3_cutoff (float, optional) – Nitrate threshold (μmol/L) below which the NO₃ correction applies. Defaults to the value stored in the posterior attributes.

  • suffix (str, optional) – Force a specific parameter suffix (e.g. "crtp"). Auto-detected by priority order when omitted.

Returns:

"x_vals" — temperature array (°C) "pN" — one key per requested percentile, e.g. "p5", "p50", "p95" "ensemble" — full array, shape (n_draws, len(temperatures)), if return_full=True "metadata" — run metadata dict, if return_full=True

Return type:

dict with keys

Predict T from proxy observations#

predict_T_from_proxyObs(proxyObs, prior_mu_t, prior_sigma_t, fwd_posterior=None, *, proxy_name=None, temptype=None, site_name=None, predictors=None, no3=None, gdgt23ratio=None, site_lat=None, site_lon=None, no3_dataset=None, no3_dataset_var='no3_sf2tc_avg', config=None, chains=4, iter_warmup=500, iter_sampling=1000, seed=42, constraint_type='unconstrained', min_temp=None, threads_per_chain=None, save_results=False, save_draws=False, filename_tag=None, cache_dir=None)[source]#

Inverse reconstruction: scaled RI → temperature percentiles.

Runs the TEXAS-Bay inverse Stan model to infer paleotemperature from observed scaled Ring Index values. Marginalises over M draws from the forward calibration posterior to propagate calibration uncertainty into the temperature reconstruction. Corresponds to Section 8 (Applications to Paleothermometry) of the manuscript.

Parameters:
  • proxyObs (array-like, shape (N,)) – Observed proxy values from downcore or coretop samples (e.g. scaledRI, TEX86).

  • prior_mu_t (float or array-like, shape (N,)) – Prior mean temperature (°C). Scalar applies the same prior to all N observations; array sets a site-specific prior per sample.

  • prior_sigma_t (float) – Prior temperature uncertainty (°C). Use a diffuse value (e.g. 10) when little prior information is available.

  • fwd_posterior (str or xr.Dataset, optional) –

    The forward calibration posterior. Accepts either:

    • str — name of the saved posterior (without .nc extension) in the posterior cache directory. The file is loaded automatically.

    • xr.Dataset — a pre-loaded posterior Dataset. No file I/O or Zenodo download is attempted; pass this when the cache is unavailable (e.g. Google Colab with a Drive-mounted .nc):

      ds = xr.open_dataset("my_drive/posterior.nc")
      result = predict_T_from_proxyObs(..., fwd_posterior=ds)
      

  • temptype (str, optional) – Temperature type: "SST" or "thermoT". Used for metadata and output file naming.

  • site_name (str, optional) – Label attached to result metadata and output filenames.

  • predictors (dict, optional) – Non-thermal predictor arrays for the N observations, e.g. {"gdgt23ratio": array, "no3": array}. Must be provided when the forward posterior was fitted with the multivariate model. Overridden by no3 / gdgt23ratio shorthands when both are given.

  • no3 (float or array-like, optional) –

    Nitrate concentration (µmol/L) for the N observations.

    • Array (length N): per-observation values — use modern WOA23 values extracted at each sample’s location (ocean_prop_ds column "no3_sf2tc_avg").

    • Scalar: broadcast to all N observations. Pass a value above no3_cutoff (e.g. no3=10.0 when no3_cutoff=1.0) to effectively disable the NO₃ correction — all observations fall outside the correction window.

    Overrides any "no3" key in predictors. Ignored when site_lat / site_lon / no3_dataset are also provided (the lookup result takes priority).

  • gdgt23ratio (float or array-like, optional) – GDGT-2/GDGT-3 ratio for the N observations. Scalar or array, same broadcast rules as no3. Overrides any "gdgt23ratio" key in predictors.

  • site_lat (float or array-like, optional) – Decimal latitude(s) of the study site(s). Scalar for a single drill core; array of length N to assign a distinct location to each observation. Requires site_lon and no3_dataset.

  • site_lon (float or array-like, optional) – Decimal longitude(s) of the study site(s). Same shape rules as site_lat.

  • no3_dataset (xr.Dataset, optional) – WOA23-derived dataset with a (lat, lon) grid, typically the ocean_prop_ds generated in the preprocessing notebook (SI_code1). Must contain no3_dataset_var. When provided together with site_lat / site_lon, the NO₃ value at those coordinates is looked up via bilinear interpolation and used as the predictor. The result is a scalar (one drill site) or array (per-obs sites), and is broadcast to all N observations when scalar.

  • no3_dataset_var (str) – Variable name to extract from no3_dataset. Default "no3_sf2tc_avg".

  • config (InvTConfig, optional) – Controls number of forward-posterior draws (M), seed, etc. Defaults to InvTConfig() (M=100).

  • chains (int) – Number of MCMC chains. Default 4.

  • iter_warmup (int) – Warmup iterations per chain. Default 500.

  • iter_sampling (int) – Sampling iterations per chain. Default 1000.

  • seed (int) – Random seed. Default 42.

  • constraint_type (str) –

    Temperature constraint applied in the Stan model:

    • "unconstrained" (default): no lower bound; P5 can be unrealistically cold near the calibration curve’s lower asymptote.

    • "hard_constraint": hard lower bound via <lower=min_temp>; prevents sub-freezing samples but the Jacobian biases P50 warm for polar sites.

    • "truncated_prior" (recommended when min_temp is set): proper truncated Normal prior via inverse-CDF reparameterization — P50 is data-driven and P5 is bounded at min_temp without warm bias.

    • "reparameterized", "soft": experimental variants.

  • min_temp (float, optional) – Lower temperature bound (°C). Required for "hard_constraint" and "truncated_prior". Typically −1.8 (seawater freezing point). When provided without an explicit constraint_type, automatically selects "truncated_prior".

  • threads_per_chain (int, optional) – Enable within-chain parallelism via Stan’s reduce_sum.

  • save_results (bool) – If True, save the quantile posterior .nc and results .npz to the invT cache directory.

  • save_draws (bool) – If True, also save the raw posterior draws (pre-quantile) as a separate {base}_draws.nc file in the invT cache directory. The file contains t_est with dims (chain, draw, obs_idx) and is suitable for kernel-density plots or custom quantile calculation. Default False.

  • filename_tag (str or list of str, optional) – Extra tag(s) appended to the output filename.

  • cache_dir (Path or str, optional) – Directory where .nc and .npz files are written when save_results or save_draws is True. Defaults to the standard invT cache (~/.texas/cache/TEXAS_invT_posterior_cache/ for pip installs, or data/cache/TEXAS_invT_posterior_cache/ in the repo).

  • proxy_name (str | None)

Returns:

"proxyObs" — input proxy array "proxy_name" — proxy type label (e.g. "scaledRI", "TEX86") "p5" — 5th percentile temperature (°C), shape (N,) "p50" — median temperature (°C), shape (N,) "p95" — 95th percentile temperature (°C), shape (N,) "metadata" — run metadata dict (model name, attrs, etc.)

Return type:

dict with keys


Download and cache#

Download posteriors#

download_posteriors(names=None, cache_dir=None, force=False)[source]#

Download forward calibration posteriors from Zenodo.

Parameters:
  • names (list of str, optional) – Subset of POSTERIOR_REGISTRY keys to download. Downloads all five posteriors when omitted (~158 MB total; the two EIV multivariate posteriors are ~78 MB each — pass names= to download only the univariate ones if you don’t need the EIV model).

  • cache_dir (Path or str, optional) – Destination directory. Defaults to the standard posterior cache.

  • force (bool) – Re-download files that already exist locally.

Returns:

Local paths of the downloaded .nc files.

Return type:

list of Path

Examples

Download only the univariate SST posterior (~0.3 MB):

>>> download_posteriors(["gen_logi_fixed_hier_crtp_univ_priorApprox_SST_scaledRI_cren3"])

Download all#

download_all(cache_dir=None, data_dir=None, force=False)[source]#

Download everything from Zenodo: forward posteriors + training data.

Files are downloaded individually; already-cached files are skipped unless force=True. Total download is ~158 MB (dominated by the two EIV multivariate posteriors at ~78 MB each).

Parameters:
  • cache_dir (Path or str, optional) – Destination for .nc posteriors. Defaults to the standard posterior cache directory.

  • data_dir (Path or str, optional) – Destination for training data files. Defaults to data/spreadsheets/.

  • force (bool) – Re-download files that already exist locally.

Return type:

None

Download training data#

download_training_data(dest_dir=None, force=False)[source]#

Download GDGT training data files from Zenodo.

Downloads the coretop/culture/mesocosm training CSVs and the CMEMS NO₃ uncertainty field used in the EIV calibration. These are needed only to re-run the SI preprocessing and calibration notebooks from scratch; they are NOT required for inverse temperature reconstructions — use download_posteriors() for that.

Parameters:
  • dest_dir (Path or str, optional) – Destination directory. Defaults to data/spreadsheets/ in the repo (or ~/.texas/data/spreadsheets/ when pip-installed).

  • force (bool) – Re-download files that already exist locally.

Returns:

Local paths of the downloaded files.

Return type:

list of Path

List posteriors#

list_posteriors(model_type='both', cache_dir=None)[source]#

List available posterior files in the cache directory.

Prints a summary and returns a dict of stem names that can be passed directly to predict_T_from_proxyObs(fwd_posterior=...).

Parameters:
  • model_type ("forward", "invT", or "both") – Which cache to inspect. Default "both".

  • cache_dir (Path or str, optional) – Override the default cache root. When given, both forward and invT subdirectories are looked for under this path.

Returns:

{"forward": [...], "invT": [...]} — lists of stem names (no .nc).

Return type:

dict

Set cache directory#

set_cache_dir(path)[source]#

Override TEXAS cache directories at runtime.

Call this before any posterior I/O. For a persistent override, set the TEXAS_CACHE_DIR environment variable instead.

Parameters:

path (str | Path) – Root directory for all TEXAS caches. Two subdirectories will be used inside it: TEXAS_posterior_cache/ and TEXAS_invT_posterior_cache/.

Return type:

None


Data builders#

Build forward data#

build_fwd_data(*, t_cul=None, proxy_cul=None, t_meso=None, proxy_meso=None, t_crtp=None, proxy_crtp=None, gdgt23ratio_crtp=None, sd_gdgt23ratio_crtp=None, no3_crtp=None, sd_no3_crtp=None, no3_cutoff=None, proxy_residuals_crtp=None, sd_proxyObs=None, R2_thermal=None, culmeso_posterior=None, prior_mean_t0=None, prior_sd_t0=None, prior_mean_k=None, prior_sd_k=None, prior_mean_b=None, prior_sd_b=None, prior_mean_v=None, prior_sd_v=None)[source]#

Build the Stan data dictionary for forward calibration models.

Handles all forward model variants:
  • culmeso / Q1_culmeso / v1_culmeso : pass t_cul, proxy_cul, t_meso, proxy_meso

  • culmesocore : add t_crtp, proxy_crtp

  • hier_crtp_multiv : add gdgt23ratio_crtp, no3_crtp

  • hier_crtp_multiv_priorApprox : add culmeso_posterior (extracts hyperpriors)

  • hier_crtp_univ_priorApprox : add culmeso_posterior (no predictors needed)

Parameters:
  • t_cul – Culture temperature and proxy arrays.

  • proxy_cul – Culture temperature and proxy arrays.

  • t_meso – Mesocosm temperature and proxy arrays.

  • proxy_meso – Mesocosm temperature and proxy arrays.

  • t_crtp – Coretop temperature and proxy arrays.

  • proxy_crtp – Coretop temperature and proxy arrays.

  • gdgt23ratio_crtp – GDGT-2/GDGT-3 ratio for coretop samples. Sets use_gdgt23ratio=1 if non-zero/non-NaN.

  • sd_gdgt23ratio_crtp – Per-site measurement SE of gdgt23ratio (same units, linear). Required for the _eiv model; always included in the data dict (defaults to zeros when not provided, which disables the G₂/₃ EIV measurement model).

  • no3_crtp – Nitrate concentration for coretop samples. Sets use_no3=1 if non-zero/non-NaN.

  • sd_no3_crtp – Per-site measurement SE of NO₃ (μmol/L, linear space). Required for the _eiv model. Always included (defaults to zeros; sites with sd=0 receive only the lognormal prior and skip the normal measurement model).

  • no3_cutoff (float | None) – NO3 threshold for the nonthermal correction. Priority: (1) this arg, (2) culmeso_posterior attrs, (3) auto-calculated via Spearman method.

  • proxy_residuals_crtp (ndarray | None) – Pre-computed proxy residuals for NO3 threshold calculation. If omitted, residuals are computed internally by fitting a generalized logistic curve. Warning: Stan models use generalized logistic — residuals from other functional forms may shift the threshold.

  • culmeso_posterior (Dataset | None) – xr.Dataset from a completed culmeso forward run. Auto-extracts prior_mean_*/prior_sd_* hyperpriors and no3_cutoff (if saved in attrs).

  • prior_mean_*/prior_sd_* – Manual hyperprior values. Override auto-extracted values from culmeso_posterior for individual params.

  • R2_thermal (float | None)

Returns:

Stan-ready data dict with proxyObs_* keys, N_* counts, use_* flags,

and hyperpriors — ready for get_posterior().

Return type:

dict

Build invT input data#

build_invT_inputData(proxyObs=None, prior_mu_t=None, prior_sigma_t=None, *, scaledRI=None, fwd_posterior_name=None, predictors=None, config=None, fwd_posterior=None)[source]#

Build the data dictionary for Stan’s inverse model and sampler configuration.

WORKFLOW: ───────── 1. Load forward calibration posterior from .nc file (or accept a pre-loaded Dataset) 2. Randomly sample M parameter sets from that posterior 3. Extract calibration curve parameters (t0, k, b, v, sigma) 4. Package optional environmental predictors (GDGT-2/3, NO3) if used 5. Return data dict (for Stan) + sampler_kwargs (for CmdStanPy)

Parameters:
  • proxyObs (ndarray | List[float]) – Observed proxy values to predict temperature from (length N). Any proxy is accepted: scaledRI, TEX86, ringIndex, etc.

  • prior_mu_t (ndarray | float) – Prior mean temperature (scalar or array of length N)

  • prior_sigma_t (float) – Prior temperature uncertainty (e.g., 10°C)

  • fwd_posterior_name (str | None) – Name of saved forward calibration (without .nc extension). Not required when fwd_posterior is supplied directly.

  • predictors (Dict[str, ndarray] | None) – Optional environmental covariates {‘gdgt23ratio’: array, ‘no3’: array}

  • config (InvTConfig | None) – Configuration object controlling M, seed, etc.

  • fwd_posterior (Dataset | None) – Pre-loaded forward posterior xr.Dataset. When provided, fwd_posterior_name is ignored and no file I/O is performed. Useful when running from Google Colab or any pip-install context where the posterior cache is not available.

  • scaledRI (ndarray | List[float])

Returns:

Dictionary for Stan’s data block sampler_kwargs: Dictionary for CmdStanPy sampling configuration

Return type:

data

WOA23 NO₃ lookup#

lookup_no3_from_woa(lat, lon, woa_dataset, variable='no3_sf2tc_avg', method='linear')[source]#

Look up modern NO₃ at one or more lat/lon coordinates from a WOA23-derived xarray Dataset.

The dataset is typically the preprocessed ocean_prop_ds generated in SI_code1, which contains thermocline-depth-integrated WOA23 climatology on a regular (lat, lon) grid. The returned value(s) are time-invariant (climatological mean) and intended as a modern-ocean proxy for the NO₃ correction in paleo reconstructions.

Parameters:
  • lat (float or array-like) – Latitude(s) in decimal degrees (−90 to 90). Pass a scalar for a single drill site; pass an array of length N to match N observations.

  • lon (float or array-like) – Longitude(s) in decimal degrees. Both −180–180 and 0–360 conventions are accepted — the function normalises to match the dataset’s convention automatically.

  • woa_dataset (xr.Dataset) – WOA23-derived dataset with a (lat, lon) grid containing variable. Dimensions must be named "lat" and "lon".

  • variable (str) – Name of the NO₃ variable to extract. Default "no3_sf2tc_avg" (thermocline depth-integrated annual average from SI_code1).

  • method ({"linear", "nearest"}) – Interpolation method. "linear" (default) performs bilinear interpolation and is preferred for smooth fields. "nearest" snaps to the closest grid cell and is useful when the dataset is sparse or has NaN-masked shelves.

Returns:

NO₃ value(s) in µmol/L. Shape matches the scalar/array input: a 0-d array for scalar inputs, 1-d array of length N for array inputs. NaN is returned for locations outside the dataset’s valid range (e.g. continental shelves masked in WOA23).

Return type:

np.ndarray

Raises:
  • KeyError – If variable is not found in woa_dataset.

  • ValueError – If woa_dataset does not have "lat" and "lon" dimensions.

Examples

Single drill site:

>>> no3_val = lookup_no3_from_woa(15.3, -23.7, ocean_prop_ds)
>>> # returns scalar-equivalent float; broadcasts to all N obs automatically
>>> result = predict_T_from_proxyObs(..., no3=no3_val)

Multi-site stack (per-obs lookup):

>>> no3_arr = lookup_no3_from_woa(core_df["lat"].values,
...                                core_df["lon"].values,
...                                ocean_prop_ds)
>>> result = predict_T_from_proxyObs(..., no3=no3_arr)

Forward calibration#

Get posterior#

get_posterior(data, stan_file, temptype, proxy_name, *, iter_warmup=None, iter_sampling=None, threads_per_chain=None, chains=None, parallel_chains=None, adapt_delta=None, max_treedepth=None, **kwargs)[source]#

Run forward calibration Stan sampling and return the posterior.

Wraps StanSampler with automatic predictor detection, CPU configuration, and metadata attachment. The returned dataset can be passed directly to predict_proxy_from_T or saved with save_posterior.

Parameters:
  • data (dict) – Stan data dict built by build_fwd_data(). Predictor flags (use_gdgt23ratio, use_no3) are auto-detected from the arrays present; you do not need to set them manually.

  • stan_file (str) – Stan model name (without .stan), e.g. "gen_logi_fixed_hier_crtp_multiv_priorApprox_eiv".

  • temptype (str) – Temperature variable type, e.g. "SST" or "thermoT". Stored in the posterior metadata.

  • proxy_name (str) – Proxy type, e.g. "scaledRI_cren3". Required — stored in the .nc attrs and validated downstream when the posterior is used for inverse reconstruction.

  • iter_warmup (int, optional) – HMC warmup iterations per chain (default: CmdStan default, 1000).

  • iter_sampling (int, optional) – Post-warmup sampling iterations per chain (default: 1000).

  • chains (int, optional) – Number of independent chains (default: 4).

  • parallel_chains (int, optional) – Chains to run simultaneously (auto-detected from CPU count).

  • threads_per_chain (int, optional) – Threads per chain for reduce_sum models (auto-enabled for models whose filename contains reduce_sum).

  • adapt_delta (float, optional) – Target acceptance rate (default: 0.8). Increase toward 0.99 to reduce divergences at the cost of more leapfrog steps.

  • max_treedepth (int, optional) – Maximum tree depth for HMC (default: 10).

  • **kwargs – Additional keyword arguments forwarded to CmdStanModel.sample.

Returns:

  • posterior (xr.Dataset) – Forward calibration posterior with parameter draws and metadata attrs (model name, temptype, proxy_name, priors, diagnostics).

  • diagnostics (str) – Human-readable sampler diagnostic summary (divergences, R-hat, ESS, E-BFMI).

Raises:

ValueError – If active predictors are present but a univariate model is requested, or if an EIV model is requested without R2_thermal.

Return type:

Tuple[Dataset, str]

Examples

>>> data = build_fwd_data(t_crtp=..., proxy_crtp=..., ...)
>>> posterior, diag = get_posterior(
...     data,
...     stan_file="gen_logi_fixed_hier_crtp_univ_priorApprox",
...     temptype="SST",
...     proxy_name="scaledRI_cren3",
... )
>>> save_posterior(posterior)

Save posterior#

save_posterior(posterior, cache_dir=None, overwrite=True, filename_suffix='')[source]#

Save a forward-model posterior to disk as compressed NetCDF.

The filename is auto-generated from the posterior’s metadata attrs: {model}_{temptype}[_gdgt23ratio][_no3_{cutoff}][_{proxy_name}]{suffix}.nc

Parameters:
  • posterior (xr.Dataset) – Forward calibration posterior returned by get_posterior(). Must have stan_model_name, temptype, and proxy_name attrs set (proxy_name is required — a warning is raised if missing).

  • cache_dir (str or Path, optional) – Directory to write the file. Defaults to the standard forward posterior cache (data/cache/TEXAS_posterior_cache/ for source installs, ~/.texas/cache/TEXAS_posterior_cache/ for pip installs).

  • overwrite (bool) – If False, raise FileExistsError when the output path already exists. Default True.

  • filename_suffix (str, optional) – Extra tag appended before .nc, e.g. "032326" for a date-stamped run. Leading/trailing underscores are stripped.

Returns:

Absolute path of the saved .nc file.

Return type:

Path

Load posterior#

load_posterior(model_name, model_type='forward', cache_dir=None)[source]#

Load a posterior from disk: {model_name}.nc in the appropriate cache directory.

Parameters:
  • model_name (str) – Name of the model file (without .nc extension)

  • model_type (Literal['forward', 'invT']) – Type of posterior (“forward” or “invT”)

  • cache_dir (str | Path | None) – Custom cache directory (overrides default locations)

Returns:

xarray.Dataset containing the posterior

Raises:

FileNotFoundError – If the posterior file doesn’t exist

Return type:

Dataset


Ensemble#

Generate ensemble (auto)#

generate_ensemble_auto(post_ds, x_vals, model_type='auto', gdgt23ratio=None, no3=None, no3_cutoff=None, return_full_ensemble=False, suffix=None, **kwargs)[source]#

Sample draws from a forward posterior and compute calibration-curve percentiles.

Inspects the posterior’s stan_model_name attr and data_vars to determine the model function, parameter names, and optional-predictor flags automatically, then delegates to generate_ensemble.

Parameters:
  • post_ds (xr.Dataset) – Forward calibration posterior returned by get_posterior() or loaded with load_posterior().

  • x_vals (np.ndarray) – Temperature values (°C) at which to evaluate the calibration curve.

  • model_type ({"auto", "forward", "inverse"}) – Force forward or inverse dispatch; "auto" (default) infers from the posterior. InvT posteriors are not supported — use predict_T_from_proxyObs() instead.

  • gdgt23ratio (np.ndarray, optional) – GDGT-2/3 ratio values; required when the posterior was fitted with a multivariate (GDGT-2/3) model.

  • no3 (np.ndarray, optional) – NO₃ concentrations (µmol/L); required when the posterior uses the NO₃ correction.

  • no3_cutoff (float, optional) – Override the NO₃ cutoff from the posterior attrs.

  • return_full_ensemble (bool) – If True, return the full M × N draw matrix in addition to percentiles. Default False.

  • suffix (str, optional) – Force a specific parameter suffix (e.g. "crtp"); overrides auto-detection.

  • **kwargs – Forwarded to generate_ensemble.

Returns:

Keys "p1""p99" (and optionally "ensemble") — each a numpy array of length len(x_vals).

Return type:

dict

Raises:

NotImplementedError – If called with an invT posterior.

Detect model and params#

detect_model_and_params(posterior_ds, suffix=None)[source]#

Auto-detect which logistic model and parameters to use. Uses a shared suffix priority via choose_suffix().

Parameters:
  • posterior_ds (Dataset)

  • suffix (str)


Diagnostics#

Sampler diagnostics#

summarize_sampler_diagnostics(fit)[source]#

Extract divergent__, treedepth__, E-BFMI, R_hat, and ESS_bulk from a CmdStanPy fit and return them as stan_diag_* attrs.

Return type:

dict

Summary table#

create_summary_table(datasets)[source]#

Build a DataFrame summarizing the stan_diag_* attrs from each xarray.Dataset.

Parameters:

datasets (list)

Return type:

DataFrame


Plotting#

Plot prior distributions#

plot_prior_distributions(priors_list=None, posterior_datasets=None, posterior_labels_list=None, show_suptitle=True, kde_bw=0.3, focus_on_posterior=True, include_groups=('t0', 'k', 'b', 'v', 'a', 'beta_G23', 'beta_NO3'), suffix_include=None, zoomin_suffix=None, zoomin_dataset_idx=None, use_linestyle_by_param=False, show_histogram=True, show_annotation=False, set_linewidth=1.5, set_fig_width_factor=3, set_fig_height_factor=3.5, set_leg_max_ncol=3, color_list=None, param_source_map=None, annotation_style='ci95', show_subplot_legend=True, show_figure_legend=True, show_prior_expression=True)[source]#

Plot priors + any number of posterior distributions in a grid, split by parameter group (t0, k, b, etc.).

Parameters:
  • param_source_map (Dict[str, int] | None) –

    Optional dict mapping a param group name to the index of the dataset in posterior_datasets that should be used as the sole source for that group. All other datasets are skipped for that group.

    Use this when different parameters come from different posteriors — e.g. logistic params (t0, k, b…) from a culmeso run and beta coefficients from a multivariate crtp run:

    plot_prior_distributions(
        posterior_datasets=[culmeso_ds, crtp_multiv_ds],
        param_source_map={"beta_G23": 1, "beta_NO3": 1},
    )
    

    When a group is not in param_source_map, all datasets are searched as usual.

  • priors_list (List[str] | Dict[str, str] | None)

  • posterior_datasets (List[xr.Dataset] | None)

  • posterior_labels_list (List[str] | None)

  • show_suptitle (bool)

  • kde_bw (float)

  • focus_on_posterior (bool)

  • include_groups (Sequence[str])

  • suffix_include (List[str] | None)

  • zoomin_suffix (str | List[str] | None)

  • zoomin_dataset_idx (int | None)

  • use_linestyle_by_param (bool)

  • show_histogram (bool)

  • show_annotation (bool)

  • set_linewidth (float)

  • set_fig_width_factor (float)

  • set_fig_height_factor (float)

  • set_leg_max_ncol (int)

  • color_list (Sequence[str] | None)

  • annotation_style (Literal['ci95', 'ci68', 'sigma'])

  • show_subplot_legend (bool)

  • show_figure_legend (bool)

  • show_prior_expression (bool)