Train and ship a public-domain aroma descriptor model

docs/AROMA.md

@@ -1,9 +1,33 @@
-# Aroma — evaluated and deferred
+# Aroma — a public-domain model, and the intensity ceiling
 
 Flavor = taste **+** aroma. The taste side ships and is clean (sweet/bitter/umami
-~0.95 AUROC + intensity + sour/salty). The aroma side is **deliberately deferred** —
-not for lack of effort, but because the public, commercially-clean data is
-insufficient. This documents the evaluation, so the decision is legible.
+~0.95 AUROC + intensity + sour/salty). A commercial-clean **aroma descriptor model now
+ships too** (see the update below); what still needs licensed/customer data is scored
+**intensity**. This documents the evaluation and the decision, so both are legible.
+
+## UPDATE (2026-07): a public-domain aroma model now ships
+
+The earlier audit missed a clean source hiding in plain sight: **PubChem's "Odor"
+annotation is public domain (HSDB / Haz-Map, US-government) for ~2,200 compounds**,
+enumerable via the annotations API. Pipeline:
+- `build_odor_notes.py --pubchem-all` — pulls that corpus (odor text + detection
+  thresholds), source-filtered to public domain, dropping GoodScents/Leffingwell/FEMA.
+- `build_aroma_dataset.py` — keyword-normalizes the free text into a controlled
+  multi-label descriptor vocabulary (presence/absence).
+- `train_aroma.py` — one RandomForest per descriptor on Morgan fingerprints (the taste
+  stack). **13 heads clear CV-AUROC ≥ 0.70**: ammoniacal 0.96, medicinal 0.96, almond
+  0.91, citrus 0.90, garlic 0.88, fishy 0.85, fruity/minty/camphor 0.84, ethereal 0.82,
+  floral 0.80, pungent 0.80. `predict_aroma()` surfaces them for **any** molecule.
+
+**Honest ceiling:** this is **presence/absence** (which notes apply), not **intensity**
+(how strong). HSDB free text carries no scored ratings, and the corpus skews industrial
+(pungent/ammoniacal/ethereal are the big classes; vanilla/caramel are sparse). A **scored
+intensity** map still needs expert-labeled panel data — licensed **PMP 2001** or the
+customer's own — the "comes with your data" upgrade. Everything below explains why that
+richer data is walled off. (The earlier `keller_2016` regressor is superseded — it scored
+negative-R² on naive-subject 0-100 ratings; the HSDB presence/absence framing works.)
+
+---
 
 ## What an aroma model needs
 

tests/test_predict.py

@@ -18,10 +18,16 @@ def test_predict_rejects_bad_smiles():
     assert "error" in out
 
 
-def test_predict_aroma_is_honest_placeholder():
+def test_predict_aroma_honest_about_model_presence():
+    # With no trained heads present (as in CI), aroma is honestly unavailable rather than
+    # fabricated; when aroma_models/ IS built, it returns a descriptor list.
     out = predict.predict_aroma("CCO")
-    assert out["available"] is False
-    assert "AROMA.md" in out["note"]
+    if predict._AROMA_MODELS:
+        assert out["available"] is True
+        assert isinstance(out.get("descriptors"), list)
+    else:
+        assert out["available"] is False
+        assert "note" in out
 
 
 def test_predict_includes_tox_screen():

training/app.py

@@ -125,9 +125,10 @@ def api_neighbors(q: Query):
     if not smi:
         return {"neighbors": []}
     res = P.substitute(smi, k=q.k)
-    for n in res.get("neighbors", []):  # enrich each candidate with a structure + a name
+    for n in res.get("neighbors", []):  # enrich each candidate with a structure + names
         n["svg"] = _svg(n["smiles"], 132, 96)
-        n["name"] = _names(n["smiles"])[0]
+        nm = _names(n["smiles"])
+        n["name"], n["iupac"] = nm[0], nm[1]
     return res
 
 
@@ -295,22 +296,29 @@ def col(name):
 
 @app.post("/api/aroma")
 def api_aroma(q: Query):
-    """Real, cited documented odor (public-domain HSDB/CAMEO) when available — NOT a trained
-    model and NOT invented scores. A trained descriptor model comes next (see docs/AROMA.md)."""
+    """Aroma read: (1) real cited DOCUMENTED odor + threshold (public-domain HSDB/Haz-Map) when
+    the molecule is in the corpus, and (2) PREDICTED descriptors from RandomForest heads trained
+    on that corpus — which work for ANY molecule, including ones with no documented entry. The
+    predicted heads are presence/absence (not intensity); each carries its CV-AUROC."""
     smi = _resolve(q.smiles)
     mol = Chem.MolFromSmiles(smi) if smi else None
     if mol is None:
         return {"available": False}
+    out = {"available": False}
     rec = _ODOR_TABLE.get(Chem.MolToInchiKey(mol).split("-")[0])
-    if not rec:
-        return {"available": False}
-    out = {"available": True}
-    if rec.get("odor"):
-        notes = [s.strip() for s in rec["odor"].split("\n") if s.strip()]
-        concise = [n for n in notes if len(n) <= 90] or notes  # lead with punchy descriptors
-        out["documented"] = {"notes": concise[:4], "source": rec.get("odor_source")}
-    if rec.get("threshold_ppm") is not None:
-        out["threshold"] = {"ppm": rec["threshold_ppm"], "source": rec.get("threshold_source")}
+    if rec:
+        if rec.get("odor"):
+            notes = [s.strip() for s in rec["odor"].split("\n") if s.strip()]
+            concise = [n for n in notes if len(n) <= 90] or notes  # lead with punchy descriptors
+            out["documented"] = {"notes": concise[:4], "source": rec.get("odor_source")}
+            out["available"] = True
+        if rec.get("threshold_ppm") is not None:
+            out["threshold"] = {"ppm": rec["threshold_ppm"], "source": rec.get("threshold_source")}
+            out["available"] = True
+    pa = P.predict_aroma(smi)
+    if pa.get("available") and pa.get("descriptors"):
+        out["predicted"] = {"descriptors": pa["descriptors"], "note": pa.get("note")}
+        out["available"] = True
     return out
 
 

training/build_aroma_dataset.py

@@ -1,60 +1,82 @@
 """
-build_aroma_dataset.py — assemble the clean aroma (odor-descriptor) training table.
+build_aroma_dataset.py — turn documented odor free-text into a multi-label descriptor dataset.
 
-Source: keller_2016 (Keller & Vosshall 2016, BMC Neuroscience, CC-BY-4.0) — the only
-commercially-clean odor-descriptor dataset available (~480 molecules, 20 descriptors).
-Each molecule gets a mean 0-100 panel rating per descriptor.
+Reads odor_notes.parquet (public-domain HSDB/Haz-Map odor descriptions + SMILES, built by
+build_odor_notes.py --pubchem-all) and normalizes each molecule's free-text odor into a
+CONTROLLED descriptor vocabulary by word-boundary keyword matching. Output aroma_train.parquet
+= (inchikey, smiles, <one 0/1 column per descriptor>) — a multi-label presence/absence dataset
+for train_aroma.py.
 
-Output: aroma_master.parquet  (smiles + 20 descriptor columns, 0-100)
+Honest scope: PRESENCE/ABSENCE learned from free text, not intensity — HSDB text has no scored
+descriptors. It's weak labeling (a missing keyword is treated as a negative), so noisier than
+expert panel data (GS-LF), but real and commercial-clean (public domain). (The earlier
+keller_2016 CC-BY approach was dropped — negative-R² on naive-subject ratings; this HSDB set is
+larger and yes/no rather than noisy 0-100 scores.)
+
+Usage: python build_aroma_dataset.py            # odor_notes.parquet -> aroma_train.parquet
 """
+import re
+import sys
 from pathlib import Path
 
 import pandas as pd
 from rdkit import Chem
 
-KELLER = Path("aroma/keller_2016")
-DESCRIPTORS = ["ACID", "AMMONIA/URINOUS", "BAKERY", "BURNT", "CHEMICAL", "COLD",
-               "DECAYED", "EDIBLE", "FISH", "FLOWER", "FRUIT", "GARLIC", "GRASS",
-               "MUSKY", "SOUR", "SPICES", "SWEATY", "SWEET", "WARM", "WOOD"]
-
-
-def canon(smiles):
-    if not isinstance(smiles, str):
-        return None
-    m = Chem.MolFromSmiles(smiles)
-    return Chem.MolToSmiles(m) if m else None
-
-
-def build():
-    mol = pd.read_csv(KELLER / "molecules.csv")
-    sti = pd.read_csv(KELLER / "stimuli.csv")
-    beh = pd.read_csv(KELLER / "behavior.csv")
-
-    beh = beh[beh["MeasurementValue"].isin(DESCRIPTORS)].copy()
-    beh["Value"] = pd.to_numeric(beh["Value"], errors="coerce")
-    beh = beh.dropna(subset=["Value"])
-
-    # mean rating per (Stimulus, descriptor) across the panel
-    sm = beh.groupby(["Stimulus", "MeasurementValue"])["Value"].mean().unstack()
-
-    # Stimulus -> CID (single-molecule) -> canonical SMILES
-    sti = sti[["Stimulus", "CIDs"]].copy()
-    sti["CID"] = pd.to_numeric(sti["CIDs"], errors="coerce")
-    sm = sm.join(sti.set_index("Stimulus")[["CID"]]).dropna(subset=["CID"])
-    sm["CID"] = sm["CID"].astype(int)
-    mol_map = mol.dropna(subset=["CID"]).drop_duplicates("CID").set_index("CID")["CanonicalSMILES"]
-    sm["smiles"] = sm["CID"].map(mol_map).map(canon)
-    sm = sm.dropna(subset=["smiles"])
+# controlled odor-descriptor vocabulary: descriptor -> keyword cues matched with word boundaries
+VOCAB = {
+    "fruity": ["fruity", "fruit"], "sweet": ["sweet"], "floral": ["floral", "flower", "flowery"],
+    "citrus": ["citrus", "lemon", "orange", "lime", "grapefruit"], "green": ["green"],
+    "minty": ["mint", "minty", "menthol", "peppermint", "spearmint"], "herbal": ["herb", "herbal"],
+    "woody": ["wood", "woody"], "spicy": ["spicy", "spice"], "rose": ["rose", "rosy"],
+    "almond": ["almond"], "vanilla": ["vanilla"], "caramel": ["caramel", "caramellic"],
+    "nutty": ["nut", "nutty"], "buttery": ["butter", "buttery"], "fatty": ["fatty"],
+    "rancid": ["rancid"], "sulfurous": ["sulfur", "sulphur", "sulfurous", "sulfury"],
+    "garlic": ["garlic"], "onion": ["onion"], "fishy": ["fish", "fishy"],
+    "earthy": ["earth", "earthy", "musty", "moldy"], "camphor": ["camphor", "camphoraceous"],
+    "pine": ["pine", "piney"], "balsamic": ["balsam", "balsamic"],
+    "medicinal": ["medicinal", "phenol", "phenolic"], "smoky": ["smoke", "smoky", "smoked"],
+    "pungent": ["pungent", "sharp", "acrid"], "ethereal": ["ether", "ethereal", "solvent"],
+    "winey": ["wine", "winey", "vinous"], "honey": ["honey"], "coconut": ["coconut"],
+    "banana": ["banana"], "apple": ["apple"], "cherry": ["cherry"], "clove": ["clove"],
+    "cinnamon": ["cinnamon", "cinnamic"], "anise": ["anise", "licorice", "aniseed"],
+    "coffee": ["coffee"], "cocoa": ["cocoa", "chocolate"], "meaty": ["meat", "meaty"],
+    "cheesy": ["cheese", "cheesy"], "creamy": ["cream", "creamy"], "waxy": ["wax", "waxy"],
+    "fresh": ["fresh"], "grassy": ["grass", "grassy", "hay"],
+    "burnt": ["burnt", "roasted", "toasted"], "tarry": ["tar", "tarry", "creosote"],
+    "ammoniacal": ["ammonia", "ammoniacal"], "fecal": ["fecal", "feces", "faecal", "manure"],
+}
+PATS = {d: re.compile(r"\b(" + "|".join(k) + r")\b", re.I) for d, k in VOCAB.items()}
 
-    # molecule level: mean across the molecule's stimuli (concentrations)
-    agg = sm.groupby("smiles")[DESCRIPTORS].mean().reset_index()
-    agg.to_parquet("aroma_master.parquet")
 
-    print(f"aroma_master.parquet: {len(agg)} molecules x {len(DESCRIPTORS)} descriptors")
-    for d in DESCRIPTORS:
-        print(f"  {d:16s} mean={agg[d].mean():5.1f}  max={agg[d].max():5.1f}  "
-              f">=15: {int((agg[d] >= 15).sum())}")
+def tag(text):
+    """Set of descriptors whose keywords appear in the odor text."""
+    t = (text or "").replace("\n", " ")
+    return {d for d, p in PATS.items() if p.search(t)}
 
 
 if __name__ == "__main__":
-    build()
+    src = "odor_notes.parquet"
+    if not Path(src).exists():
+        print(f"{src} not found — run build_odor_notes.py --pubchem-all first")
+        sys.exit(1)
+    df = pd.read_parquet(src)
+    df = df[df["odor"].notna() & df["smiles"].notna()].copy()
+    rows = []
+    for _, r in df.iterrows():
+        if Chem.MolFromSmiles(str(r["smiles"])) is None:
+            continue
+        tags = tag(r["odor"])
+        if not tags:  # no descriptor keyword -> uninformative ("characteristic odor"); skip to
+            continue  # limit false-negative noise (absence-as-negative only among tagged mols)
+        row = {"inchikey": r["inchikey"], "smiles": r["smiles"]}
+        for d in VOCAB:
+            row[d] = int(d in tags)
+        rows.append(row)
+    out = pd.DataFrame(rows)
+    out.to_parquet("aroma_train.parquet")
+    counts = {d: int(out[d].sum()) for d in VOCAB}
+    print(f"aroma_train.parquet: {len(out)} molecules with >=1 descriptor "
+          f"(of {len(df)} with odor text), {len(VOCAB)} descriptors")
+    print("descriptor positives (sorted):")
+    for d, n in sorted(counts.items(), key=lambda x: -x[1]):
+        print(f"  {n:5d}  {d}")

training/build_odor_notes.py

@@ -169,15 +169,17 @@ def annotation_records(heading):
         page += 1
 
 
-def cids_to_inchikeys(cids, chunk=100):
-    """{cid: InChIKey} via batched PubChem property calls (~1 request per 100 CIDs)."""
+def cids_to_structs(cids, chunk=100):
+    """{cid: (InChIKey, CanonicalSMILES)} via batched PubChem property calls (~1 per 100 CIDs).
+    SMILES is needed downstream to fingerprint these molecules for the aroma descriptor model."""
     out = {}
     for i in range(0, len(cids), chunk):
         ch = cids[i:i + chunk]
-        d = _get(f"{_BASE}/pug/compound/cid/{','.join(map(str, ch))}/property/InChIKey/JSON")
+        d = _get(f"{_BASE}/pug/compound/cid/{','.join(map(str, ch))}"
+                 f"/property/InChIKey,SMILES/JSON")
         for p in (d or {}).get("PropertyTable", {}).get("Properties", []):
-            if p.get("InChIKey"):
-                out[p["CID"]] = p["InChIKey"]
+            if p.get("InChIKey"):  # PubChem renamed the SMILES fields; accept either
+                out[p["CID"]] = (p["InChIKey"], p.get("SMILES") or p.get("ConnectivitySMILES"))
     return out
 
 
@@ -193,16 +195,17 @@ def build_from_pubchem_annotations():
         thr_by_cid.setdefault(cid, []).extend((src, t) for t in strs)
     all_cids = sorted(set(odor_by_cid) | set(thr_by_cid))
     print(f"annotations: {len(odor_by_cid)} odor + {len(thr_by_cid)} threshold CIDs "
-          f"(public-domain); resolving {len(all_cids)} InChIKeys...", flush=True)
-    cid2ik = cids_to_inchikeys(all_cids)
+          f"(public-domain); resolving {len(all_cids)} structures...", flush=True)
+    cid2s = cids_to_structs(all_cids)
     rows = []
     for cid in all_cids:
-        ik = cid2ik.get(cid)
-        if not ik:
+        st = cid2s.get(cid)
+        if not st:
             continue
+        ik, smi = st
         ostrs, osrcs = odor_by_cid.get(cid, ([], set()))
         thr_pairs = thr_by_cid.get(cid, [])
-        rows.append({"inchikey": ik,
+        rows.append({"inchikey": ik, "smiles": smi,
                      "odor": "\n".join(sorted(set(ostrs), key=len)) or None,
                      "odor_source": "; ".join(sorted(osrcs)) or None,
                      "odor_threshold_ppm": _parse_threshold_ppm(thr_pairs),
@@ -211,7 +214,7 @@ def build_from_pubchem_annotations():
     return rows
 
 
-COLS = ["inchikey", "odor", "odor_source", "odor_threshold_ppm",
+COLS = ["inchikey", "smiles", "odor", "odor_source", "odor_threshold_ppm",
         "odor_threshold_note", "odor_threshold_source"]
 
 if __name__ == "__main__":

training/predict.py

@@ -193,6 +193,19 @@ def _ik1(*smiles):
     for p in _TOX_DIR.glob("*_rf.joblib"):
         _TOX_MODELS[p.stem.replace("_rf", "")] = joblib.load(p)
 
+# Odor-descriptor heads (public-domain HSDB corpus, presence/absence). Loaded if trained.
+# Real, commercial-clean aroma predictions — NOT intensity (see docs/AROMA.md).
+_AROMA_MODELS = {}
+_AROMA_META = {}
+_AROMA_DIR = Path("aroma_models")
+if _AROMA_DIR.exists():
+    for p in _AROMA_DIR.glob("*_clf.joblib"):
+        _AROMA_MODELS[p.stem.replace("_clf", "")] = joblib.load(p)
+    _mf = _AROMA_DIR / "manifest.json"
+    if _mf.exists():
+        import json as _json
+        _AROMA_META = _json.loads(_mf.read_text()).get("descriptors", {})
+
 # Known-label lookup: ground truth for molecules we actually have data on. This
 # is how the salty/sour data works as a FLAG without a model — if a queried
 # molecule is in our labeled set, we report the verified fact instead of a guess.
@@ -594,17 +607,29 @@ def analyze_balance(ingredients):
     }
 
 
-def predict_aroma(smiles, top_k=8):
-    """Aroma is deferred. No commercially-clean *public* odor data yields a working
-    model (see docs/AROMA.md), so rather than fabricate smells we return an honest
-    'not available'. A real head gets trained on licensed (PMP 2001) or customer
-    odor data — OpenPOM's MIT architecture for large sets, RandomForest for small —
-    and wired in here then. predict() never calls this unless include_aroma=True."""
+def predict_aroma(smiles, top_k=8, threshold=0.5):
+    """Predicted odor descriptors from RandomForest heads trained on the PUBLIC-DOMAIN HSDB
+    odor corpus (see docs/AROMA.md). Returns the descriptors the model scores above threshold,
+    each with its probability and the head's CV-AUROC. This is PRESENCE/ABSENCE (the free-text
+    corpus carries no intensity), not a scored intensity map — honest about that ceiling; a
+    stronger intensity model needs licensed (PMP 2001) or customer panel data. Returns
+    available:False until the heads are trained into aroma_models/ (train_aroma.py)."""
     m = Chem.MolFromSmiles(smiles)
     if m is None:
         return {"error": f"unparseable SMILES: {smiles}"}
-    return {"available": False,
-            "note": "aroma deferred — needs licensed/customer odor data; see docs/AROMA.md"}
+    if not _AROMA_MODELS:
+        return {"available": False,
+                "note": "aroma model not trained here — build with train_aroma.py"}
+    fp = _fp(m)
+    preds = []
+    for name, clf in _AROMA_MODELS.items():
+        p = float(clf.predict_proba(fp)[0][1])
+        if p >= threshold:
+            preds.append({"odor": name, "score": round(p, 3),
+                          "auroc": _AROMA_META.get(name, {}).get("auroc")})
+    preds.sort(key=lambda d: -d["score"])
+    return {"available": True, "predicted": True, "descriptors": preds[:top_k],
+            "note": "presence/absence model on public-domain HSDB odor text; not intensity"}
 
 
 # Plain-language meaning of each Tox21 assay, for caution context.