1616
1717import matplotlib .pyplot as plt
1818import numpy as np
19- from rdkit import Chem
20- from rdkit .Chem import Descriptors
21- from tqdm .auto import tqdm
22-
23-
24- def _exact_mass (smi : str ) -> float :
25- mol = Chem .MolFromSmiles (smi )
26- return Descriptors .ExactMolWt (mol ) if mol else float ("nan" )
27-
28-
29- def _heavy_atoms (smi : str ) -> int :
30- mol = Chem .MolFromSmiles (smi )
31- return mol .GetNumHeavyAtoms () if mol else 0
3219
3320
3421def main ():
@@ -39,7 +26,6 @@ def main():
3926
4027 out = Path (args .output )
4128 out .parent .mkdir (parents = True , exist_ok = True )
42- fig2_path = out .with_name (out .stem + "_pop.png" )
4329
4430 print ("Loading CSV ..." )
4531 rows = []
@@ -53,7 +39,6 @@ def main():
5339 simba_pred = np .array ([float (r ["simba_pred_mces" ]) for r in rows ])
5440 simba_gt = np .array ([float (r ["simba_gt_mces" ]) for r in rows ])
5541 oracle_pred = np .array ([float (r ["oracle_pred_mces" ]) for r in rows ])
56- oracle_gt = np .array ([float (r ["oracle_gt_mces" ]) for r in rows ])
5742 simba_tani = np .array ([float (r ["simba_tanimoto" ]) for r in rows ])
5843 covered = np .array ([r ["covered" ] == "True" for r in rows ])
5944
@@ -71,22 +56,6 @@ def main():
7156 cosine_oracle = 1.0 - oracle_pred / 40.0
7257 cosine_label = "Cosine sim (embedding, 1−pred/40) ← re-run diagnose_retrieval to get spectral"
7358
74- # ── RDKit properties (cached) ─────────────────────────────────────────────
75- print ("Computing molecular properties ..." )
76- mass_cache : dict [str , float ] = {}
77- atom_cache : dict [str , int ] = {}
78- unique_smi = {r ["test_smi" ] for r in rows } | {r ["simba_smi" ] for r in rows }
79- for s in tqdm (unique_smi , desc = "RDKit" ):
80- mass_cache [s ] = _exact_mass (s )
81- atom_cache [s ] = _heavy_atoms (s )
82-
83- mass_test = np .array ([mass_cache [r ["test_smi" ]] for r in rows ])
84- mass_simba = np .array ([mass_cache [r ["simba_smi" ]] for r in rows ])
85- atoms_test = np .array ([atom_cache [r ["test_smi" ]] for r in rows ])
86- atoms_simba = np .array ([atom_cache [r ["simba_smi" ]] for r in rows ])
87- mass_diff = np .abs (mass_test - mass_simba )
88- atoms_diff = np .abs (atoms_test - atoms_simba ).astype (float )
89-
9059 om = covered # covered mask
9160
9261 # ── Print summary stats ───────────────────────────────────────────────────
@@ -97,30 +66,13 @@ def main():
9766 print (
9867 f" { 'spectral' if has_spectral else 'embedding' } cosine oracle pick: mean={ cosine_oracle [om ].mean ():.3f} median={ np .median (cosine_oracle [om ]):.3f} "
9968 )
100- print (
101- f" mass_diff pair: mean={ mass_diff [om ].mean ():.1f} Da median={ np .median (mass_diff [om ]):.1f} Da"
102- )
103- print (
104- f" atoms_diff pair: mean={ atoms_diff [om ].mean ():.1f} median={ np .median (atoms_diff [om ]):.1f} "
105- )
10669 print (f" Tanimoto SIMBA pair: mean={ simba_tani [om ].mean ():.3f} " )
10770
108- HIGH = om & (simba_err > 25 )
109- LOW = om & (simba_err < 5 )
110- print (
111- f"\n High error (>25): n={ HIGH .sum ():,} mass_diff mean={ mass_diff [HIGH ].mean ():.1f} Da "
112- f"atoms_diff mean={ atoms_diff [HIGH ].mean ():.1f} tani mean={ simba_tani [HIGH ].mean ():.3f} "
113- )
114- print (
115- f" Low error (<5): n={ LOW .sum ():,} mass_diff mean={ mass_diff [LOW ].mean ():.1f} Da "
116- f"atoms_diff mean={ atoms_diff [LOW ].mean ():.1f} tani mean={ simba_tani [LOW ].mean ():.3f} "
117- )
118-
119- # ── Figure 1: Error anatomy (2×3) ─────────────────────────────────────────
120- fig , axes = plt .subplots (2 , 3 , figsize = (18 , 11 ))
71+ # ── Figure 1: panels 1 + 6 ────────────────────────────────────────────────
72+ fig , axes = plt .subplots (1 , 2 , figsize = (14 , 5 ))
12173
12274 # 1. Cosine sim distribution: SIMBA pick vs oracle pick
123- ax = axes [0 , 0 ]
75+ ax = axes [0 ]
12476 bins1 = np .linspace (0 , 1 , 51 )
12577 ax .hist (
12678 cosine_oracle [om ],
@@ -144,50 +96,8 @@ def main():
14496 ax .legend (fontsize = 9 )
14597 ax .grid (True , alpha = 0.2 )
14698
147- # 2. simba_err vs mass_diff (hexbin)
148- ax = axes [0 , 1 ]
149- hb = ax .hexbin (mass_diff [om ], simba_err [om ], gridsize = 30 , cmap = "Reds" , mincnt = 1 )
150- plt .colorbar (hb , ax = ax , label = "count" )
151- ax .axhline (0 , color = "k" , lw = 1 , ls = "--" )
152- ax .set_xlabel ("|mass_test − mass_SIMBA_pick| (Da)" )
153- ax .set_ylabel ("simba_err (GT MCES − pred MCES)" )
154- ax .set_title ("2 · Calibration error vs mass difference" , fontweight = "bold" )
155- ax .grid (True , alpha = 0.2 )
156-
157- # 3. simba_err vs atoms_diff
158- ax = axes [0 , 2 ]
159- hb = ax .hexbin (atoms_diff [om ], simba_err [om ], gridsize = 30 , cmap = "Purples" , mincnt = 1 )
160- plt .colorbar (hb , ax = ax , label = "count" )
161- ax .axhline (0 , color = "k" , lw = 1 , ls = "--" )
162- ax .set_xlabel ("|atoms_test − atoms_SIMBA_pick|" )
163- ax .set_ylabel ("simba_err" )
164- ax .set_title ("3 · Calibration error vs atom count diff" , fontweight = "bold" )
165- ax .grid (True , alpha = 0.2 )
166-
167- # 4. simba_err vs atoms_test
168- ax = axes [1 , 0 ]
169- hb = ax .hexbin (
170- atoms_test [om ].astype (float ), simba_err [om ], gridsize = 30 , cmap = "Blues" , mincnt = 1
171- )
172- plt .colorbar (hb , ax = ax , label = "count" )
173- ax .axhline (0 , color = "k" , lw = 1 , ls = "--" )
174- ax .set_xlabel ("Heavy atom count of test molecule" )
175- ax .set_ylabel ("simba_err" )
176- ax .set_title ("4 · Calibration error vs test mol size" , fontweight = "bold" )
177- ax .grid (True , alpha = 0.2 )
178-
179- # 5. simba_err vs Tanimoto(test, SIMBA pick)
180- ax = axes [1 , 1 ]
181- hb = ax .hexbin (simba_tani [om ], simba_err [om ], gridsize = 30 , cmap = "Oranges" , mincnt = 1 )
182- plt .colorbar (hb , ax = ax , label = "count" )
183- ax .axhline (0 , color = "k" , lw = 1 , ls = "--" )
184- ax .set_xlabel ("Tanimoto similarity (test, SIMBA pick)" )
185- ax .set_ylabel ("simba_err" )
186- ax .set_title ("5 · Calibration error vs structural similarity" , fontweight = "bold" )
187- ax .grid (True , alpha = 0.2 )
188-
18999 # 6. GT MCES distribution of SIMBA picks with pred mean
190- ax = axes [1 , 2 ]
100+ ax = axes [1 ]
191101 bins6 = np .arange (0 , 42.5 , 2.5 )
192102 ax .hist (simba_gt [om ], bins = bins6 , color = "#5B8DB8" , edgecolor = "none" , alpha = 0.85 )
193103 ax .axvline (
@@ -218,113 +128,7 @@ def main():
218128 )
219129 plt .tight_layout ()
220130 plt .savefig (out , dpi = 140 , bbox_inches = "tight" )
221- print (f"\n Figure 1 saved → { out } " )
222-
223- # ── Figure 2: Population analysis (2×2) ──────────────────────────────────
224- fig2 , axes2 = plt .subplots (2 , 2 , figsize = (14 , 11 ))
225-
226- # A. mass_diff distribution: high vs low error
227- ax = axes2 [0 , 0 ]
228- bins_m = np .linspace (0 , np .percentile (mass_diff [om ], 99 ), 60 )
229- ax .hist (
230- mass_diff [LOW ],
231- bins = bins_m ,
232- color = "#4E9A7A" ,
233- alpha = 0.75 ,
234- density = True ,
235- edgecolor = "none" ,
236- label = f"err<5 n={ LOW .sum ():,} μ={ mass_diff [LOW ].mean ():.1f} Da" ,
237- )
238- ax .hist (
239- mass_diff [HIGH ],
240- bins = bins_m ,
241- color = "#E07B54" ,
242- alpha = 0.75 ,
243- density = True ,
244- edgecolor = "none" ,
245- label = f"err>25 n={ HIGH .sum ():,} μ={ mass_diff [HIGH ].mean ():.1f} Da" ,
246- )
247- ax .set_xlabel ("|mass_test − mass_SIMBA_pick| (Da)" )
248- ax .set_ylabel ("density" )
249- ax .set_title ("A · Mass diff: low-error vs high-error pairs" , fontweight = "bold" )
250- ax .legend (fontsize = 9 )
251- ax .grid (True , alpha = 0.2 )
252-
253- # B. 2D hexbin (mass_test, mass_simba) colored by mean simba_err
254- ax = axes2 [0 , 1 ]
255- lim = np .percentile (np .concatenate ([mass_test [om ], mass_simba [om ]]), 99 )
256- hb = ax .hexbin (
257- mass_test [om ],
258- mass_simba [om ],
259- C = simba_err [om ],
260- gridsize = 35 ,
261- cmap = "RdYlGn_r" ,
262- reduce_C_function = np .mean ,
263- mincnt = 3 ,
264- )
265- plt .colorbar (hb , ax = ax , label = "mean simba_err" )
266- ax .plot ([0 , lim ], [0 , lim ], "k--" , lw = 1.2 , alpha = 0.7 , label = "test = SIMBA pick" )
267- ax .set_xlim (0 , lim )
268- ax .set_ylim (0 , lim )
269- ax .set_xlabel ("Exact mass — test molecule (Da)" )
270- ax .set_ylabel ("Exact mass — SIMBA pick (Da)" )
271- ax .set_title ("B · Mass pairs colored by mean calibration error" , fontweight = "bold" )
272- ax .legend (fontsize = 8 )
273- ax .grid (True , alpha = 0.2 )
274-
275- # C. Tanimoto distribution: high vs low error
276- ax = axes2 [1 , 0 ]
277- bins_t = np .linspace (0 , 1 , 41 )
278- ax .hist (
279- simba_tani [LOW ],
280- bins = bins_t ,
281- color = "#4E9A7A" ,
282- alpha = 0.75 ,
283- density = True ,
284- edgecolor = "none" ,
285- label = f"err<5 μ={ simba_tani [LOW ].mean ():.3f} " ,
286- )
287- ax .hist (
288- simba_tani [HIGH ],
289- bins = bins_t ,
290- color = "#E07B54" ,
291- alpha = 0.75 ,
292- density = True ,
293- edgecolor = "none" ,
294- label = f"err>25 μ={ simba_tani [HIGH ].mean ():.3f} " ,
295- )
296- ax .set_xlabel ("Tanimoto similarity (test, SIMBA pick)" )
297- ax .set_ylabel ("density" )
298- ax .set_title ("C · Tanimoto: low-error vs high-error pairs" , fontweight = "bold" )
299- ax .legend (fontsize = 9 )
300- ax .grid (True , alpha = 0.2 )
301-
302- # D. SIMBA pick GT vs oracle GT colored by simba_err
303- ax = axes2 [1 , 1 ]
304- hb = ax .hexbin (
305- oracle_gt [om ],
306- simba_gt [om ],
307- C = simba_err [om ],
308- gridsize = 30 ,
309- cmap = "RdYlGn_r" ,
310- reduce_C_function = np .mean ,
311- mincnt = 3 ,
312- )
313- plt .colorbar (hb , ax = ax , label = "mean simba_err" )
314- ax .plot ([0 , 40 ], [0 , 40 ], "k--" , lw = 1.2 , alpha = 0.7 , label = "SIMBA = oracle" )
315- ax .set_xlabel ("Oracle GT MCES" )
316- ax .set_ylabel ("SIMBA pick GT MCES" )
317- ax .set_title ("D · SIMBA GT vs oracle GT (colored by error)" , fontweight = "bold" )
318- ax .legend (fontsize = 8 )
319- ax .grid (True , alpha = 0.2 )
320-
321- fig2 .suptitle (
322- "Population Analysis: High-error (err>25) vs Low-error (err<5) SIMBA Retrieval Pairs" ,
323- fontsize = 12 ,
324- )
325- plt .tight_layout ()
326- plt .savefig (fig2_path , dpi = 140 , bbox_inches = "tight" )
327- print (f"Figure 2 saved → { fig2_path } " )
131+ print (f"\n Saved → { out } " )
328132
329133
330134if __name__ == "__main__" :
0 commit comments