"""Paleobiology Database first appearances for the Part 2 innovations, with confidence intervals on the true time of origin. Cleaning of the occurrence downloads (applied identically to every innovation): * drop form taxa, ichnotaxa and other flagged records (PBDB 'flags' not empty); * drop open-nomenclature identifications ('?', 'cf.', 'aff.', 'indet.'); * keep only records resolved to genus or species; * drop records whose age bracket is wider than MAX_SPAN Myr (stage-level resolution required); * drop the accepted names listed in the innovation's pbdb_exclude, each with a stated reason (alt claims kept for the sensitivity run, and database misassignments such as the lungfish 'Pentlandia' in the Silurian). First appearance = the stage (max_ma, min_ma) of the oldest remaining record. Range extension (how much older than its first fossil the innovation may be): Strauss & Sadler (1989) classical confidence interval applied to the K_OLDEST oldest distinct horizons (formation x interval) of the clade, i.e. the gap before the first fossil is assumed to be drawn from the same distribution as the gaps between the earliest horizons. Using only the earliest horizons, rather than the whole range, is the simple way to respect Marshall's (1997) point that recovery potential is not uniform through a range: early members of a clade are rare, so the early gap sizes are the relevant ones. Ages are drawn uniformly inside each stage-level bracket, so the stage uncertainty is propagated. """ import io import os import numpy as np import pandas as pd from config import RAW, DATA, SEED from innovations import INNOVATIONS MAX_SPAN = 15.0 # Myr; widest age bracket accepted for an occurrence K_OLDEST = 10 # number of oldest horizons used for the gap distribution EXT_CAP = 50.0 # Myr; cap on a single range-extension draw (heavy tail when few horizons) NDRAW = 4000 def read_pbdb_csv(path): with open(path, encoding="utf-8") as f: lines = f.read().split("\n") start = [k for k, l in enumerate(lines) if l.startswith('"occurrence_no"') or l.startswith('"orig_no"') or l.startswith('"interval_no"')] if not start: return pd.DataFrame() return pd.read_csv(io.StringIO("\n".join(lines[start[0]:])), low_memory=False) def clean(df, exclude): n0 = len(df) if n0 == 0: return df, dict(n_raw=0) d = df.copy() d["flags"] = d["flags"].fillna("").astype(str) d = d[d["flags"].str.strip() == ""] idn = d["identified_name"].fillna("").astype(str) bad = idn.str.contains(r"\?|cf\.|aff\.|indet\.|informal", regex=True) d = d[~bad] d = d[d["accepted_rank"].isin(["species", "genus", "subgenus"])] d = d[(d["max_ma"] - d["min_ma"]) <= MAX_SPAN] acc = d["accepted_name"].fillna("").astype(str) for name in (exclude or {}): keep = ~acc.str.startswith(name) for col in ("class", "order", "family", "genus"): if col in d.columns: keep &= ~(d[col].fillna("").astype(str) == name) d = d[keep] acc = d["accepted_name"].fillna("").astype(str) return d, dict(n_raw=n0, n_clean=len(d)) def horizons(d): """One row per distinct horizon: formation (or collection when no formation) x interval.""" h = d.copy() h["hz"] = h["formation"].fillna("").astype(str) blank = h["hz"].str.strip() == "" h.loc[blank, "hz"] = "coll" + h.loc[blank, "collection_no"].astype(str) h["hz"] = h["hz"] + "|" + h["early_interval"].astype(str) g = h.groupby("hz").agg(max_ma=("max_ma", "max"), min_ma=("min_ma", "min")).reset_index() return g.sort_values("max_ma", ascending=False).reset_index(drop=True) def extension_draws(hz, rng, n=NDRAW): """Strauss-Sadler range-extension draws (Myr) from the K_OLDEST oldest horizons.""" if len(hz) < 3: return None hi = hz["max_ma"].to_numpy() lo = hz["min_ma"].to_numpy() ages = rng.uniform(lo, hi, size=(n, len(hz))) ages = -np.sort(-ages, axis=1) k = min(K_OLDEST, ages.shape[1]) span = ages[:, 0] - ages[:, k - 1] u = rng.uniform(size=n) ext = span * ((1.0 - u) ** (-1.0 / (k - 1)) - 1.0) return np.minimum(ext, EXT_CAP) def build(verbose=False): rng = np.random.default_rng(SEED) rng_lit = np.random.default_rng(SEED + 77) # separate stream for the all-literature-dates variant rows, ext_store, fad_store, lit_store = [], {}, {}, {} for inn in INNOVATIONS: r = dict(id=inn["id"], name=inn["name"], o2_class=inn["o2_class"], lit_age=inn["lit_age"], lit_old=inn["lit_old"], lit_young=inn["lit_young"], alt_old=inn["alt_old"], verified=inn["verified"], lit_source=inn["lit_source"], requirement=inn["req"]) if inn.get("pbdb"): raw = read_pbdb_csv(os.path.join(RAW, "pbdb", "occs_%s.csv" % inn["id"])) d, info = clean(raw, inn.get("pbdb_exclude")) r.update(pbdb_taxa="+".join(inn["pbdb"]), pbdb_n_raw=info.get("n_raw", 0), pbdb_n_clean=len(d)) if len(d): d = d.sort_values(["max_ma", "min_ma"], ascending=False) first = d.iloc[0] hz = horizons(d) r.update(pbdb_fad_old=float(first["max_ma"]), pbdb_fad_young=float(first["min_ma"]), pbdb_fad_taxon=str(first["accepted_name"]), pbdb_fad_interval=str(first["early_interval"]), pbdb_fad_formation=str(first.get("formation", "")), pbdb_n_horizons=len(hz)) ext = extension_draws(hz, rng) if ext is not None: ext_store[inn["id"]] = ext r.update(ext50=float(np.median(ext)), ext95=float(np.quantile(ext, 0.95))) if verbose: print("====", inn["id"], info, "horizons", len(hz)) print(d[["accepted_name", "identified_name", "early_interval", "max_ma", "min_ma", "formation"]].head(6).to_string(index=False)) rows.append(r) tab = pd.DataFrame(rows) # adopted first-appearance bracket: PBDB when it exists, the literature otherwise has = tab["pbdb_fad_old"].notna() tab["fad_old"] = np.where(has, tab["pbdb_fad_old"], tab["lit_old"]) tab["fad_young"] = np.where(has, tab["pbdb_fad_young"], tab["lit_young"]) tab["fad_basis"] = np.where(has, "PBDB", "literature") # flag innovations whose cleaned PBDB first appearance (bracket midpoint) is more than 5 Myr younger than # the literature date (the stage-resolution and open-nomenclature cuts can remove the oldest accepted # material, e.g. Frasnian tetrapods and Carnian-early Norian pterosaurs) tab["pbdb_minus_lit_Myr"] = np.where(has, 0.5 * (tab["pbdb_fad_old"] + tab["pbdb_fad_young"]) - tab["lit_age"], np.nan) tab["pbdb_younger_than_literature"] = has & (tab["pbdb_minus_lit_Myr"] < -5.0) # pooled extension distribution for innovations without enough horizons of their own pooled = np.concatenate([v for v in ext_store.values()]) if ext_store else np.zeros(1) for i, r in tab.iterrows(): if r["id"] not in ext_store: ext_store[r["id"]] = rng.choice(pooled, size=NDRAW) tab.loc[i, "ext50"] = float(np.median(ext_store[r["id"]])) tab.loc[i, "ext95"] = float(np.quantile(ext_store[r["id"]], 0.95)) tab.loc[i, "ext_basis"] = "pooled" else: tab.loc[i, "ext_basis"] = "own horizons" fad_store[r["id"]] = rng.uniform(r["fad_young"], r["fad_old"], size=NDRAW) lit_store[r["id"]] = rng_lit.uniform(r["lit_young"], r["lit_old"], size=NDRAW) tab["fad_mid"] = 0.5 * (tab["fad_old"] + tab["fad_young"]) tab["origin50"] = [float(np.median(fad_store[i] + ext_store[i])) for i in tab["id"]] tab["origin95_old"] = [float(np.quantile(fad_store[i] + ext_store[i], 0.95)) for i in tab["id"]] tab.to_csv(os.path.join(DATA, "innovations_dated.csv"), index=False) draws = {i: (fad_store[i], ext_store[i], lit_store[i]) for i in tab["id"]} return tab, draws if __name__ == "__main__": t, _ = build(verbose=True) pd.set_option("display.width", 250) print(t[["id", "lit_age", "lit_old", "lit_young", "pbdb_fad_old", "pbdb_fad_young", "pbdb_fad_taxon", "pbdb_n_clean", "pbdb_n_horizons", "ext50", "ext95", "ext_basis", "origin50"]].to_string())