Skip to content

Commit 05a32a7

Browse files
feat: split structural neighbors vs profile-based substitutes (#222)
Two distinct nearest-neighbor searches, previously conflated under 'substitutes': - Structural neighbors (/api/neighbors) — Tanimoto over Morgan fingerprints. The look-alikes. - Substitutes (/api/substitutes) — NEW: cosine over predicted 6 taste + 164 aroma head scores. A molecule that tastes and smells like the target is the real drop-in swap regardless of structure. Vanillin -> ethyl vanillin #1 (0.94), then isovanillin, coniferaldehyde. Widen the reference set from 3,845 taste-labelled molecules to the full ~8.8k universe so any molecule can surface (ethyl vanillin wasn't a candidate before). The 170-head profile matrix is precomputed by build_profile_index.py -> profile_index.npz (gitignored) for instant startup. Surfaces synced: workbench two sections; MCP find_structural_neighbors + find_substitutes; skill; README. 22 unit tests pass. Closes #213. Signed-off-by: Austin L. <86896075+rvnminers-A-and-N@users.noreply.github.com>
1 parent b8046c2 commit 05a32a7

9 files changed

Lines changed: 276 additions & 66 deletions

File tree

.gitignore

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,3 +67,4 @@ Thumbs.db
6767
!training/samples/*.csv
6868
!training/samples/*.md
6969
training/resolve_cache.json
70+
training/profile_index.npz

README.md

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -80,8 +80,11 @@ is tagged by how it was derived**, so nothing reads as more certain than its sou
8080
acrylamide, ethyl carbamate, furan, and more) and an OAV dosing-balance analysis
8181
that flags the component about to overpower a blend (quantitative when threshold
8282
tables are loaded).
83-
- **Substitution search** — nearest-neighbor lookup over the labeled set for
84-
reformulation and cost-down ("find me a molecule that behaves like this one").
83+
- **Substitutes & structural neighbors** — two nearest-neighbor searches over the whole
84+
molecule universe for reformulation and cost-down: **substitutes** rank by *taste + aroma
85+
profile* match (a molecule that tastes and smells like the target — e.g. ethyl vanillin for
86+
vanillin — regardless of structure), while **structural neighbors** rank by Tanimoto/Morgan
87+
structure similarity (the look-alikes).
8588
- **Flavor Studio** — one hub to pick any mix of everyday **flavors** (banana,
8689
saffron, pumpkin, bubble gum…) *and* **notes** (citrus, floral…) → ranked food-safe
8790
molecules + drop-in swaps. A flavor *is* a set of notes, so they live in one picker.

mcp-server/server.py

Lines changed: 21 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -96,8 +96,15 @@ def _read_full(molecule: str) -> dict:
9696
return d
9797

9898

99-
def _find_substitutes(molecule: str, k: int) -> dict:
99+
def _find_structural_neighbors(molecule: str, k: int) -> dict:
100100
r = _post("/api/neighbors", {"smiles": molecule, "k": k})
101+
if r.status_code != 200:
102+
return {"error": f"structural-neighbor search failed (HTTP {r.status_code})"}
103+
return _strip(r.json())
104+
105+
106+
def _find_substitutes(molecule: str, k: int) -> dict:
107+
r = _post("/api/substitutes", {"smiles": molecule, "k": k})
101108
if r.status_code != 200:
102109
return {"error": f"substitute search failed (HTTP {r.status_code})"}
103110
return _strip(r.json())
@@ -253,12 +260,23 @@ def read_full(molecule: str) -> dict:
253260

254261
@mcp.tool()
255262
def find_substitutes(molecule: str, k: int = 8) -> dict:
256-
"""Find the k structurally-nearest molecules (drop-in swaps / reformulation candidates)
257-
to the given molecule, each with its similarity and known tastes.
263+
"""Find the k best SUBSTITUTES — molecules whose predicted taste+aroma PROFILE is closest to
264+
the given molecule (cosine over the head scores). These are the drop-in swaps: a molecule that
265+
tastes and smells like the target, regardless of structure (e.g. ethyl vanillin for vanillin).
266+
Each with its profile_match, known tastes, and aromas.
258267
"""
259268
return _find_substitutes(molecule, k)
260269

261270

271+
@mcp.tool()
272+
def find_structural_neighbors(molecule: str, k: int = 8) -> dict:
273+
"""Find the k STRUCTURAL neighbors — molecules most similar in structure (Tanimoto / Morgan
274+
fingerprint) to the given molecule. Structural look-alikes (contrast find_substitutes, which
275+
matches by taste+aroma profile). Each with its similarity and known tastes.
276+
"""
277+
return _find_structural_neighbors(molecule, k)
278+
279+
262280
@mcp.tool()
263281
def list_stereoisomers(molecule: str) -> dict:
264282
"""List every stereoisomer (R/S centers and E/Z bonds) of a molecule, with any

skills/flavormancer/scripts/flavormancer.py

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -118,6 +118,12 @@ def cmd_flavor(a):
118118

119119

120120
def cmd_substitutes(a):
121+
# profile-based swaps: closest taste+aroma head-score match (e.g. ethyl vanillin for vanillin)
122+
return _req("/api/substitutes", "POST", {"smiles": a.molecule, "k": a.k})
123+
124+
125+
def cmd_structural_neighbors(a):
126+
# structural look-alikes: Tanimoto / Morgan nearest neighbors
121127
return _req("/api/neighbors", "POST", {"smiles": a.molecule, "k": a.k})
122128

123129

@@ -188,9 +194,12 @@ def mol(name):
188194
mol("read").set_defaults(fn=cmd_read)
189195
mol("read-full").set_defaults(fn=cmd_read_full)
190196
mol("stereoisomers").set_defaults(fn=cmd_stereoisomers)
191-
s = mol("substitutes")
197+
s = mol("substitutes") # taste+aroma profile match (the drop-in swaps)
192198
s.add_argument("-k", type=int, default=8)
193199
s.set_defaults(fn=cmd_substitutes)
200+
s = mol("structural-neighbors") # Tanimoto structural look-alikes
201+
s.add_argument("-k", type=int, default=8)
202+
s.set_defaults(fn=cmd_structural_neighbors)
194203

195204
s = sub.add_parser("formulate")
196205
s.add_argument("ingredients", nargs="+")

tests/test_substitute.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ def test_substitute_rejects_bad_smiles():
1616

1717

1818
def test_substitute_graceful_without_data(monkeypatch):
19-
monkeypatch.setattr(predict, "_SUB_INDEX", ([], [], [], []))
19+
monkeypatch.setattr(predict, "_SUB_INDEX", ([], [], [], [], None, []))
2020
out = predict.substitute("CCO")
2121
assert out["neighbors"] == []
2222
assert "note" in out
@@ -26,7 +26,7 @@ def test_substitute_ranks_by_similarity(monkeypatch):
2626
mols = ["CCO", "CCCO", "c1ccccc1"] # ethanol, propanol, benzene
2727
fps = [predict._MORGAN.GetFingerprint(Chem.MolFromSmiles(s)) for s in mols]
2828
canon = [_canon(s) for s in mols]
29-
monkeypatch.setattr(predict, "_SUB_INDEX", (fps, canon, [[], [], []], [[], [], []]))
29+
monkeypatch.setattr(predict, "_SUB_INDEX", (fps, canon, [[], [], []], [[], [], []], None, []))
3030
out = predict.substitute("CCO", k=2)
3131
neighbors = out["neighbors"]
3232
# the query itself is excluded

training/app.py

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -300,6 +300,24 @@ def api_neighbors(q: Query):
300300
return res
301301

302302

303+
@app.post("/api/substitutes")
304+
def api_substitutes(q: Query):
305+
"""Profile-based substitutes: molecules whose predicted taste+aroma head scores line up
306+
closest with the query — the taste/smell-alikes (vs /api/neighbors' structural look-alikes)."""
307+
smi = _resolve(q.smiles)
308+
if not smi:
309+
return {"substitutes": []}
310+
res = P.substitutes(smi, k=q.k)
311+
for n in res.get("substitutes", []): # same enrichment as neighbors: structure + names + aroma + GRAS
312+
n["svg"] = _svg(n["smiles"], 132, 96)
313+
nm = _names(n["smiles"])
314+
n["name"], n["iupac"] = nm[0], nm[1]
315+
n["aroma"] = _aroma_tags_cheap(n["smiles"], n.pop("aromas", []))
316+
_m = Chem.MolFromSmiles(n["smiles"])
317+
n["gras"] = bool(_m is not None and Chem.MolToInchiKey(_m).split("-")[0] in P._GRAS)
318+
return res
319+
320+
303321
@app.post("/api/names")
304322
def api_names(q: Query):
305323
"""Common (PubChem Title) + IUPAC names for the queried molecule."""

training/build_profile_index.py

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
"""build_profile_index.py — precompute the neighbor / substitute reference index.
2+
3+
The profile-based substitutes rank molecules by cosine similarity over their predicted taste +
4+
aroma head SCORES. Computing that 170-head inference over the whole ~8.8k-molecule universe at app
5+
startup is slow (~3 min), so we precompute it here once and cache to profile_index.npz. predict.py
6+
loads it instantly (rebuilding the cheap Morgan fingerprints from SMILES on load).
7+
8+
Output profile_index.npz:
9+
smiles (N,) canonical SMILES, deduped by connectivity skeleton
10+
taste_documented (N,) comma-separated documented tastes ("" if none)
11+
aromas (N,) comma-separated confident aroma heads (score >= 0.5)
12+
profiles (N,D) float32 head-score matrix: taste heads then aroma heads
13+
dims (D,) column labels ("taste:sweet", "aroma:citrus", ...)
14+
15+
Usage: python build_profile_index.py
16+
"""
17+
import numpy as np
18+
import pandas as pd
19+
import predict as P
20+
from rdkit import Chem
21+
22+
SRC = "master_enrichment.parquet"
23+
OUT = "profile_index.npz"
24+
25+
26+
def main():
27+
m = pd.read_parquet(SRC)
28+
taste_heads, aroma_heads = P._profile_heads()
29+
smis, td, feats, seen = [], [], [], set()
30+
for _, r in m.iterrows():
31+
mol = Chem.MolFromSmiles(str(r["smiles"]))
32+
if mol is None:
33+
continue
34+
skel = Chem.MolToInchiKey(mol).split("-")[0]
35+
if skel in seen:
36+
continue
37+
seen.add(skel)
38+
smis.append(Chem.MolToSmiles(mol))
39+
td.append(str(r.get("taste_documented") or ""))
40+
feats.append(P._feat(mol)[0])
41+
X = np.vstack(feats)
42+
print(f"{len(smis)} unique structures — running {len(taste_heads)}+{len(aroma_heads)} heads...", flush=True)
43+
cols, aromas = [], [[] for _ in smis]
44+
for t in taste_heads:
45+
cols.append(P._CLASSIFIERS[t].predict_proba(X)[:, 1])
46+
for a in aroma_heads:
47+
col = P._AROMA_MODELS[a].predict_proba(X)[:, 1]
48+
cols.append(col)
49+
for i in range(len(smis)):
50+
if col[i] >= 0.5:
51+
aromas[i].append(a)
52+
profiles = np.column_stack(cols).astype("float32")
53+
dims = [f"taste:{t}" for t in taste_heads] + [f"aroma:{a}" for a in aroma_heads]
54+
np.savez_compressed(
55+
OUT,
56+
smiles=np.array(smis, dtype=object),
57+
taste_documented=np.array(td, dtype=object),
58+
aromas=np.array([",".join(a) for a in aromas], dtype=object),
59+
profiles=profiles,
60+
dims=np.array(dims, dtype=object),
61+
)
62+
print(f"wrote {OUT}: {profiles.shape[0]} molecules x {profiles.shape[1]} heads")
63+
64+
65+
if __name__ == "__main__":
66+
main()

training/predict.py

Lines changed: 111 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -239,6 +239,11 @@ def _load_rf(path):
239239
# molecule is in our labeled set, we report the verified fact instead of a guess.
240240
_KNOWN = {} # inchikey -> {taste: 1}
241241
_MASTER = Path("taste_master.parquet")
242+
# The neighbor / substitute reference set: the FULL molecule universe (every structure we know,
243+
# ~8.8k) so structural neighbors and profile substitutes can surface ANY molecule — e.g. ethyl
244+
# vanillin as the top vanillin substitute — not just the taste-labelled subset. Falls back to
245+
# taste_master when the enrichment table hasn't been built yet.
246+
_UNIVERSE = Path("master_enrichment.parquet")
242247
if _MASTER.exists():
243248
import pandas as pd
244249
_m = pd.read_parquet(_MASTER)
@@ -1061,57 +1066,99 @@ def _taste_profile(out):
10611066
# tool (swap an expensive or supply-constrained ingredient for a close analogue,
10621067
# with its known tastes shown). This is the clean Track-A core; the product
10631068
# (Track B, #22) mirrors it as a pgvector ANN query over the same fingerprints.
1064-
_SUB_INDEX = None # lazily built: (fps, smiles, known_tastes, predicted_aromas)
1069+
# lazily built: (fps, smiles, known_tastes, predicted_aromas, profiles, profile_dims)
1070+
# profiles: an (N x D) float32 matrix of predicted head SCORES — taste heads then aroma heads —
1071+
# the "flavor profile" vector used for profile-based substitutes (vs the fingerprint fps used for
1072+
# structural neighbors). profile_dims labels the columns.
1073+
_SUB_INDEX = None
10651074
_SUB_LOCK = _threading.Lock() # guards the one-time index build against concurrent callers
10661075

10671076

1077+
def _profile_heads():
1078+
"""The ordered head list backing a flavor-profile vector: taste heads then aroma heads.
1079+
Same order is used at index-build and query time so the vectors line up."""
1080+
return sorted(_CLASSIFIERS), list(_AROMA_MODELS)
1081+
1082+
10681083
def _build_sub_index():
10691084
global _SUB_INDEX
1085+
import numpy as np
1086+
# Fast path: load the precomputed profile index (build_profile_index.py). The 170-head
1087+
# inference over ~8.8k molecules is slow (~3 min); the cache makes startup instant. We only
1088+
# rebuild the cheap Morgan fingerprints from SMILES on load.
1089+
cache = Path("profile_index.npz")
1090+
if cache.exists():
1091+
z = np.load(cache, allow_pickle=True)
1092+
smis = [str(s) for s in z["smiles"]]
1093+
tastes = [[t for t in str(s).split(",") if t.strip()] for s in z["taste_documented"]]
1094+
aromas = [[a for a in str(s).split(",") if a.strip()] for s in z["aromas"]]
1095+
fps = [_MORGAN.GetFingerprint(Chem.MolFromSmiles(s)) for s in smis]
1096+
_SUB_INDEX = (fps, smis, tastes, aromas, z["profiles"], list(z["dims"]))
1097+
return
10701098
fps, smis, tastes, feats = [], [], [], []
1071-
if _MASTER.exists():
1072-
import numpy as np
1099+
profiles, profile_dims = None, []
1100+
src = _UNIVERSE if _UNIVERSE.exists() else _MASTER
1101+
if src.exists():
10731102
import pandas as pd
1074-
m = pd.read_parquet(_MASTER)
1103+
m = pd.read_parquet(src)
10751104
basic = [t for t in ("sweet", "bitter", "umami", "sour", "salty") if t in m.columns]
1105+
has_documented = "taste_documented" in m.columns # enrichment: comma-separated string
1106+
seen_skel = set()
10761107
for _, r in m.iterrows():
10771108
mol = Chem.MolFromSmiles(str(r["smiles"]))
10781109
if mol is None:
10791110
continue
1111+
skel = Chem.MolToInchiKey(mol).split("-")[0]
1112+
if skel in seen_skel: # dedupe by connectivity so the universe doesn't repeat a molecule
1113+
continue
1114+
seen_skel.add(skel)
10801115
fps.append(_MORGAN.GetFingerprint(mol))
10811116
smis.append(Chem.MolToSmiles(mol))
1082-
tastes.append([t for t in basic if r.get(t) == 1])
1117+
if has_documented:
1118+
tastes.append([t for t in str(r.get("taste_documented") or "").split(",") if t.strip()])
1119+
else:
1120+
tastes.append([t for t in basic if r.get(t) == 1])
10831121
feats.append(_feat(mol)[0])
1084-
# predicted aroma descriptors per molecule (batched over the 16 heads) — so the palette
1085-
# match can score aroma as well as taste, once, at first use
1122+
# One-time batched inference over the whole labeled set. We keep BOTH:
1123+
# - the thresholded confident-aroma NAME list per molecule (for display / palette match)
1124+
# - the full head-SCORE matrix (taste heads + aroma heads) for profile-based substitutes.
1125+
# Kept at n_jobs=1: a single big predict_proba is already vectorized C; joblib fan-out here
1126+
# thrashed under concurrency. Reusing this index is what makes the endpoints fast.
10861127
aromas = [[] for _ in smis]
1087-
if _AROMA_MODELS and feats:
1128+
if feats:
10881129
X = np.vstack(feats)
1089-
# One-time batch over the whole labeled set (~8k rows x 24 heads). Kept at n_jobs=1:
1090-
# a single big predict_proba is already vectorized C, and joblib fan-out here just
1091-
# thrashed under concurrency. Reusing these aromas is what makes the endpoint fast.
1092-
for name, clf in _AROMA_MODELS.items():
1093-
col = clf.predict_proba(X)[:, 1]
1130+
taste_heads, aroma_heads = _profile_heads()
1131+
cols = [_CLASSIFIERS[t].predict_proba(X)[:, 1] for t in taste_heads]
1132+
for name in aroma_heads:
1133+
col = _AROMA_MODELS[name].predict_proba(X)[:, 1]
1134+
cols.append(col)
10941135
for i in range(len(smis)):
10951136
if col[i] >= 0.5:
10961137
aromas[i].append(name)
1138+
profile_dims = [f"taste:{t}" for t in taste_heads] + [f"aroma:{a}" for a in aroma_heads]
1139+
profiles = np.column_stack(cols).astype("float32") if cols else None
10971140
else:
10981141
aromas = []
1099-
_SUB_INDEX = (fps, smis, tastes, aromas)
1142+
_SUB_INDEX = (fps, smis, tastes, aromas, profiles, profile_dims)
11001143

11011144

1102-
def substitute(smiles: str, k: int = 8, min_similarity: float = 0.0) -> dict:
1103-
"""Nearest-neighbor substitution: the k labeled molecules most structurally
1104-
similar to the query (Tanimoto over Morgan fingerprints), each with its known
1105-
tastes. The reformulation / cost-down tool — swap an ingredient for a close
1106-
analogue. Returns {'neighbors': [...]} ranked by similarity (self excluded)."""
1107-
mol = Chem.MolFromSmiles(smiles)
1108-
if mol is None:
1109-
return {"error": f"unparseable SMILES: {smiles}"}
1145+
def _ensure_sub_index():
11101146
if _SUB_INDEX is None: # double-checked lock: a request during the startup build waits for
11111147
with _SUB_LOCK: # that one build instead of kicking off a second (which would thrash cores)
11121148
if _SUB_INDEX is None:
11131149
_build_sub_index()
1114-
fps, smis, tastes, _aromas = _SUB_INDEX
1150+
1151+
1152+
def structural_neighbors(smiles: str, k: int = 8, min_similarity: float = 0.0) -> dict:
1153+
"""STRUCTURAL neighbors: the k labeled molecules most structurally similar to the query
1154+
(Tanimoto over Morgan fingerprints), each with its known tastes. Structural look-alikes —
1155+
contrast with profile-based `substitutes` (taste/aroma-alikes). Returns {'neighbors': [...]}
1156+
ranked by similarity (self excluded)."""
1157+
mol = Chem.MolFromSmiles(smiles)
1158+
if mol is None:
1159+
return {"error": f"unparseable SMILES: {smiles}"}
1160+
_ensure_sub_index()
1161+
fps, smis, tastes, _aromas, _profiles, _dims = _SUB_INDEX
11151162
if not fps:
11161163
return {"neighbors": [], "note": "no reference set loaded (taste_master.parquet absent)"}
11171164
q = _MORGAN.GetFingerprint(mol)
@@ -1138,7 +1185,46 @@ def substitute(smiles: str, k: int = 8, min_similarity: float = 0.0) -> dict:
11381185
if len(neighbors) >= k:
11391186
break
11401187
return {"query": self_smi, "neighbors": neighbors,
1141-
"basis": "Tanimoto / Morgan r2 2048-bit over labeled molecules"}
1188+
"basis": "structural — Tanimoto / Morgan r2 2048-bit over labeled molecules"}
1189+
1190+
1191+
# back-compat alias: the endpoint / callers historically called this `substitute`
1192+
substitute = structural_neighbors
1193+
1194+
1195+
def substitutes(smiles: str, k: int = 8) -> dict:
1196+
"""PROFILE-based substitutes: the k molecules whose predicted FLAVOR profile (taste + aroma
1197+
head scores) is closest to the query's — the drop-in reformulation list. A molecule that
1198+
*tastes and smells* like the target is a likely substitute regardless of its structure, so
1199+
this ranks by cosine similarity over the head-score vectors (not fingerprint distance).
1200+
Returns {'substitutes': [...]} ranked by profile match (self excluded)."""
1201+
import numpy as np
1202+
mol = Chem.MolFromSmiles(smiles)
1203+
if mol is None:
1204+
return {"error": f"unparseable SMILES: {smiles}"}
1205+
_ensure_sub_index()
1206+
_fps, smis, tastes, aromas, profiles, _dims = _SUB_INDEX
1207+
if profiles is None or not len(smis):
1208+
return {"substitutes": [], "note": "no reference set / models loaded"}
1209+
x = _feat(mol)
1210+
taste_heads, aroma_heads = _profile_heads()
1211+
qv = np.array([_CLASSIFIERS[t].predict_proba(x)[0, 1] for t in taste_heads]
1212+
+ [_AROMA_MODELS[a].predict_proba(x)[0, 1] for a in aroma_heads], dtype="float32")
1213+
qn = qv / (float(np.linalg.norm(qv)) + 1e-9)
1214+
pn = profiles / (np.linalg.norm(profiles, axis=1, keepdims=True) + 1e-9)
1215+
sims = pn @ qn
1216+
self_skel = Chem.MolToInchiKey(mol).split("-")[0]
1217+
subs = []
1218+
for i in np.argsort(-sims):
1219+
ni = Chem.MolFromSmiles(smis[i])
1220+
if ni is None or Chem.MolToInchiKey(ni).split("-")[0] == self_skel:
1221+
continue
1222+
subs.append({"smiles": smis[i], "profile_match": round(float(sims[i]), 3),
1223+
"known_tastes": tastes[i], "aromas": aromas[i] if i < len(aromas) else []})
1224+
if len(subs) >= k:
1225+
break
1226+
return {"query": Chem.MolToSmiles(mol), "substitutes": subs,
1227+
"basis": "profile — cosine over predicted taste + aroma head scores"}
11421228

11431229

11441230
def palette_match(tastes, aromas=None, k=5):
@@ -1150,7 +1236,7 @@ def palette_match(tastes, aromas=None, k=5):
11501236
predicted descriptors."""
11511237
if _SUB_INDEX is None:
11521238
_build_sub_index()
1153-
_, smis, tlist, alist = _SUB_INDEX
1239+
_, smis, tlist, alist, _profiles, _dims = _SUB_INDEX
11541240
t_target, a_target = set(tastes or []), set(aromas or [])
11551241
if not (t_target or a_target) or not smis:
11561242
return {"target": {"tastes": sorted(t_target), "aromas": sorted(a_target)}, "matches": []}

0 commit comments

Comments
 (0)