Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 66 additions & 0 deletions docs/AROMA.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
# Aroma — evaluated and deferred

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.

## What an aroma model needs

A model mapping a molecule's **structure → odor descriptors** (sweet, floral,
woody…). The state of the art is the **Principal Odor Map** (Lee et al., *Science*
2023; Google/Osmo), reimplemented open-source as **OpenPOM** (`BioMachineLearning/
openpom`, MIT) — a message-passing GNN trained on ~5,000 **expert-labeled** molecules
(the GS-LF dataset).

## The licensing wall (audited)

OpenPOM's *code* is MIT, but its *training data* is restricted. A full audit of the
Pyrfume catalogue (see [`DATA-SOURCES.md`](DATA-SOURCES.md)) found that **every rich
odor-descriptor dataset is proprietary or NonCommercial** — Leffingwell, GoodScents,
Arctander, Flavornet (©Datu), FlavorDB / FooDB (NC), OlfactionBase ("all rights
reserved"), AromaDB (CSIR), Dravnieks (ASTM ©), sharma_2021 (ACS ©), snitz_2019
(CC-BY-NC). The **only** commercially-clean odor-descriptor set is **`keller_2016`**
(Keller & Vosshall 2016, *BMC Neuroscience*, CC-BY-4.0; ~480 molecules, 20 descriptors).

## The empirical result

We aggregated `keller_2016` ([`training/build_aroma_dataset.py`](../training/build_aroma_dataset.py))
and trained one RandomForest regressor per descriptor on 2048-bit Morgan fingerprints
([`training/train_aroma.py`](../training/train_aroma.py)), scored by **honest 5-fold
cross-validation** (400 trees):

| descriptor | CV-R² | descriptor | CV-R² | descriptor | CV-R² | descriptor | CV-R² |
|---|---|---|---|---|---|---|---|
| acid | −0.24 | cold | −0.19 | fruit | −0.05 | sour | −0.05 |
| ammonia | −0.15 | decayed | −0.20 | garlic | −0.12 | spices | −0.12 |
| bakery | −0.22 | edible | −0.19 | grass | −0.45 | sweaty | +0.03 |
| burnt | −0.13 | fish | −0.04 | musky | −0.23 | sweet | −0.05 |
| chemical | −0.10 | flower | −0.12 | warm | −0.28 | wood | −0.08 |

**All 20 descriptors scored CV-R² ≤ 0** (range −0.45 to +0.03) — every model is
*worse than predicting the mean*. **0/20 usable heads.**

## Root cause

`keller_2016` is **naive-subject** data: random volunteers rating *unfamiliar*
molecules on a 0–100 scale. People can't reliably name what a molecule smells like,
so the labels are noise (the per-descriptor means compress to ~22–30 for nearly every
molecule). The *learnable* odor data uses **expert** labels — which is exactly the
data that's restricted. **That is the structural reason the entire field trains on
GS-LF**, and why a clean public aroma model isn't currently possible.

## Decision

We do **not** ship a negative-R² model — it would output confident, wrong smells, and
a flavor chemist would catch it instantly (worse than nothing). `predict_aroma()`
returns an honest "not available" marker, and **we lead with the taste engine**.

The **OpenPOM scaffold (`training/train_odor.py`) is kept** as the aroma engine for
when clean fuel exists:
1. **License PMP 2001** (~$2,775, Leffingwell & Associates) → re-curate → train.
2. **A customer's own odor data** (the paid pilot) → train on-prem.
3. A future open *expert-labeled* dataset, if one emerges.

Until then: **aroma comes with your data** — which is on-thesis (public data proves
the method; the customer's data unlocks the rest).
1 change: 1 addition & 0 deletions docs/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,5 +8,6 @@ Architecture, capability catalogue, and design decisions.
- `SOURCES.md` — data sources, libraries, research, and license attribution
- `API-CONTRACT.md` — the fixed JSON contract the .NET API exposes and the React UI consumes
- `DATA-SOURCES.md` — data acquisition tracker: what each source unlocks and how to get it
- `AROMA.md` — why the aroma model is deferred (clean-data audit + the empirical evaluation)

These document *what we're building and why*, so decisions don't get re-argued.
1 change: 1 addition & 0 deletions docs/SOURCES.md
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ record of provenance, not legal advice. Get an IP/OSS-license review before ship

**Aroma**
- **Pyrfume** + **Leffingwell / GoodScents (GS-LF)** odor datasets — the usual training data behind OpenPOM, but **RESTRICTED and NOT USED**: Leffingwell's manifest cites use restrictions (*John Leffingwell & Google*); GoodScents/Arctander/Flavornet (© Datu Inc.) are likewise proprietary. We **exclude all of them** (the demo may go to a customer / commercial use). The aroma model will use only commercial-clean **open** odor data (CC-BY sets like `keller_2016`; smaller — see `DATA-SOURCES.md`). The OpenPOM *code* is MIT.
- **keller_2016** — Keller & Vosshall (2016), *BMC Neuroscience*, **CC-BY-4.0**; ~480 molecules with naive-subject odor-descriptor ratings. The only commercially-clean odor-descriptor set — evaluated for the aroma model and found too noisy to learn from (CV-R² ≤ 0 across all 20 descriptors; see `docs/AROMA.md`).

**Safety / regulatory (lookups — data-gated)**
- **FEMA GRAS list** — usual/maximum use levels for the dosing analyzer. (FEMA.)
Expand Down
60 changes: 60 additions & 0 deletions training/build_aroma_dataset.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
"""
build_aroma_dataset.py — assemble the clean aroma (odor-descriptor) training table.

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.

Output: aroma_master.parquet (smiles + 20 descriptor columns, 0-100)
"""
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"])

# 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())}")


if __name__ == "__main__":
build()
66 changes: 66 additions & 0 deletions training/train_aroma.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
"""
train_aroma.py — odor-descriptor regressors (one per descriptor) on Morgan fingerprints.

Reads aroma_master.parquet (keller_2016, CC-BY); trains a RandomForestRegressor per
descriptor predicting its 0-100 panel rating; reports HONEST 5-fold CV R2 (small,
noisy panel-mean data — expect modest, descriptor-dependent scores). Saves the ones
that clear a minimum CV-R2 to aroma_models/.
"""
from pathlib import Path

import joblib
import numpy as np
import pandas as pd
from rdkit import Chem
from rdkit.Chem import DataStructs, rdFingerprintGenerator
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score

FP_BITS, FP_RADIUS = 2048, 2
_MORGAN = rdFingerprintGenerator.GetMorganGenerator(radius=FP_RADIUS, fpSize=FP_BITS)
MIN_R2 = 0.10 # below this a descriptor is too noisy to ship an honest head
OUT = Path("aroma_models")
OUT.mkdir(exist_ok=True)
for s in OUT.glob("*_reg.joblib"):
s.unlink()

DESCRIPTORS = ["ACID", "AMMONIA/URINOUS", "BAKERY", "BURNT", "CHEMICAL", "COLD",
"DECAYED", "EDIBLE", "FISH", "FLOWER", "FRUIT", "GARLIC", "GRASS",
"MUSKY", "SOUR", "SPICES", "SWEATY", "SWEET", "WARM", "WOOD"]


def fp(smiles):
m = Chem.MolFromSmiles(smiles)
if m is None:
return None
bv = _MORGAN.GetFingerprint(m)
arr = np.zeros((FP_BITS,), dtype=np.int8)
DataStructs.ConvertToNumpyArray(bv, arr)
return arr


df = pd.read_parquet("aroma_master.parquet")
feats, keep = [], []
for i, s in enumerate(df["smiles"]):
f = fp(s)
if f is not None:
feats.append(f)
keep.append(i)
X = np.array(feats)
df = df.iloc[keep].reset_index(drop=True)

print(f"training odor-descriptor regressors on {len(df)} molecules:")
kept = 0
for d in DESCRIPTORS:
y = df[d].values
mask = ~np.isnan(y)
Xd, yd = X[mask], y[mask]
r2 = cross_val_score(RandomForestRegressor(n_estimators=400, n_jobs=-1, random_state=42),
Xd, yd, cv=5, scoring="r2").mean()
flag = "kept" if r2 >= MIN_R2 else "drop (too noisy)"
if r2 >= MIN_R2:
reg = RandomForestRegressor(n_estimators=400, n_jobs=-1, random_state=42).fit(Xd, yd)
joblib.dump(reg, OUT / f"{d.replace('/', '_')}_reg.joblib")
kept += 1
print(f" {d:16s} n={int(mask.sum()):3d} CV-R2={r2:+.2f} -> {flag}")
print(f"\nkept {kept}/{len(DESCRIPTORS)} descriptor heads (CV-R2 >= {MIN_R2}) -> aroma_models/")
Loading