nx_chem_fingerprint.nx source
↩ module page · 344 lines · 14054 B
1// nx_chem_fingerprint.nx -- C2.6 milestone: ECFP-style Morgan
2// circular-substructure fingerprint + Tanimoto similarity.
3//
4// Algorithm: extended-connectivity hashing (Rogers & Hahn 2010).
5// 1. Each atom gets an initial invariant from chemistry-relevant
6// features: (heavy-neighbor count, valence-sum, atomic number,
7// formal charge, implicit-H count, in-ring flag, aromaticity).
8// 2. For each refinement iteration r in 1..radius:
9// new_inv[i] = poly_hash(old_inv[i], [(bond_order, neighbor_old_inv)] sorted)
10// The sorted-list-of-(bond, neighbor) construction makes the
11// hash atom-order-invariant; the polynomial multiplier 31 is
12// pinned by source identity ("ECFP4_v1_mul_31").
13// 3. After radius iterations, accumulate ALL invariants seen
14// across all atoms at all radii into the fingerprint bit-set.
15// 4. Fold into a bit vector of length N_BITS by setting
16// bit[h & (N_BITS-1)] for each accumulated invariant h.
17//
18// EXCEED axis hit (landscape doc):
19// E5 -- hash-pinned fingerprint definitions. Algorithm identity is
20// the source bytes of this file. RDKit's Morgan/ECFP doesn't
21// match the ECFP-paper Figure 8 (RDKit Issue #2018); the nx
22// substrate ships content-addressed semantics: ECFP4_v1
23// computed in 2026 matches ECFP4_v1 computed in 2030 because
24// the function bytes are bit-identical.
25//
26// Honest gaps (deferred):
27// - Stereo-aware ECFP (folds Atom.stereo + Bond.stereo into invariants):
28// C2.6.1 future iteration; current MVP is stereo-blind.
29// - FCFP (functional-class equivalence rather than atomic-equivalence):
30// C2.7 future primitive.
31// - Variable bit length (current 1024 only; 512 / 2048 / 4096 hash-pinned
32// variants are simple parameter changes when needed).
33//
34// nx_safety_envelope:
35// intended_use: chemoinformatic similarity + substructure
36// hashing; foundation for nx_chem_similarity_search
37// and nx_chem_chemometrics descriptor matrix
38// sil_target: SIL1
39// evidence: [C2.6 KAT in nx_chem_fingerprint_test.nx;
40// algorithm hash-pinned by source bytes;
41// Tanimoto identity/symmetry verified by KAT]
42// verdict: BENCH-PENDING
43
44import "nx_chem.nx"
45import "nx_chem_molecule.nx"
46const NX_MAGIC_10000: i64 = 10000
47
48// ---- algorithm identity constants (hash-pinned by source) ----
49const NX_CHEM_ECFP_V1_MUL: nx_int = 31 // polynomial-hash multiplier
50const NX_CHEM_ECFP_V1_INIT_SEED: nx_int = 5381 // djb2-style seed
51const NX_CHEM_ECFP_DEFAULT_BITS: nx_int = 1024
52const NX_CHEM_ECFP_DEFAULT_RADIUS: nx_int = 2 // ECFP4 standard
53
54// ---- bit-vector storage uses u64 slots ----
55const NX_CHEM_FP_BITS_PER_SLOT: nx_int = 64
56
57// =================================================================
58// Initial atom invariant for ECFP. Packs chemistry-relevant features
59// into a single i64 deterministically.
60//
61// Layout (LSB first):
62// bits 0-7 : z (atomic number, 0..118)
63// bits 8-11 : degree (heavy-neighbor count, 0..15)
64// bits 12-15 : charge + 8 (encoding -8..+7)
65// bits 16-19 : h_count + 2 (encoding -2..+13; -1 sentinel as 1)
66// bit 20 : aromaticity flag (any model)
67// bit 21 : in_ring flag
68// bits 22-29 : valence-sum (0..255)
69//
70// Differences from morgan canonical-rank invariant:
71// - includes valence-sum (sum of bond orders incident)
72// - includes in_ring (Morgan canonical doesn't need this)
73// - omits isotope (less load-bearing for substructure similarity;
74// stereo-aware variant in C2.6.1 will optionally include it)
75// =================================================================
76func nx_chem_fp_initial_invariant(m: *MolGraph, atom_idx: nx_int) -> nx_int {
77 let a: *Atom = ((m.atoms as nx_int) + (atom_idx * NX_ATOM_BYTES)) as *Atom
78 var deg: nx_int = 0
79 var vsum: nx_int = 0
80 var bi: nx_int = 0
81 while bi < m.n_bonds {
82 let b: *Bond = ((m.bonds as nx_int) + (bi * NX_BOND_BYTES)) as *Bond
83 if b.a == atom_idx {
84 deg = deg + 1
85 if b.order == NX_BOND_SINGLE { vsum = vsum + 1 }
86 if b.order == NX_BOND_DOUBLE { vsum = vsum + 2 }
87 if b.order == NX_BOND_TRIPLE { vsum = vsum + 3 }
88 if b.order == NX_BOND_QUAD { vsum = vsum + 4 }
89 if b.order == NX_BOND_AROMATIC { vsum = vsum + 1 }
90 }
91 if b.b == atom_idx {
92 deg = deg + 1
93 if b.order == NX_BOND_SINGLE { vsum = vsum + 1 }
94 if b.order == NX_BOND_DOUBLE { vsum = vsum + 2 }
95 if b.order == NX_BOND_TRIPLE { vsum = vsum + 3 }
96 if b.order == NX_BOND_QUAD { vsum = vsum + 4 }
97 if b.order == NX_BOND_AROMATIC { vsum = vsum + 1 }
98 }
99 bi = bi + 1
100 }
101 if deg > 15 { deg = 15 }
102 if vsum > 255 { vsum = 255 }
103 var v: nx_int = a.z & 0xff
104 v = v | ((deg & 0xf) << 8)
105 var c: nx_int = a.charge + 8
106 if c < 0 { c = 0 }
107 if c > 15 { c = 15 }
108 v = v | ((c & 0xf) << 12)
109 var h: nx_int = a.h_count + 2
110 if h < 0 { h = 0 }
111 if h > 15 { h = 15 }
112 v = v | ((h & 0xf) << 16)
113 var arom: nx_int = 0
114 if a.aromaticity != 0 { arom = 1 }
115 v = v | ((arom & 0x1) << 20)
116 v = v | ((a.in_ring & 0x1) << 21)
117 v = v | ((vsum & 0xff) << 22)
118 return v
119}
120
121// =================================================================
122// Insertion sort for an nx_int array ascending (defensive `while
123// done == 0` pattern per session-codified gotcha).
124// =================================================================
125func nx_chem_fp_sort_asc(arr: *nx_int, n: nx_int) -> nx_int {
126 var i: nx_int = 1
127 while i < n {
128 let key: nx_int = arr[i]
129 var j: nx_int = i - 1
130 var done: nx_int = 0
131 while done == 0 {
132 if j < 0 { done = 1 }
133 else {
134 if arr[j] <= key { done = 1 }
135 else {
136 arr[j + 1] = arr[j]
137 j = j - 1
138 }
139 }
140 }
141 arr[j + 1] = key
142 i = i + 1
143 }
144 return 0
145}
146
147// =================================================================
148// Polynomial hash of (own_invariant, sorted [(bond_order, neighbor_inv)] list).
149//
150// h = seed
151// h = h * 31 + own_inv
152// for each (b, n) in sorted_pairs:
153// h = h * 31 + b
154// h = h * 31 + n
155// =================================================================
156func nx_chem_fp_hash_neighborhood(
157 own_inv: nx_int,
158 pair_bonds: *nx_int, pair_nbrs: *nx_int, n_pairs: nx_int
159) -> nx_int {
160 var h: nx_int = NX_CHEM_ECFP_V1_INIT_SEED
161 h = (h * NX_CHEM_ECFP_V1_MUL) + own_inv
162 var i: nx_int = 0
163 while i < n_pairs {
164 h = (h * NX_CHEM_ECFP_V1_MUL) + pair_bonds[i]
165 h = (h * NX_CHEM_ECFP_V1_MUL) + pair_nbrs[i]
166 i = i + 1
167 }
168 return h
169}
170
171// =================================================================
172// One refinement step: for each atom, compute new invariant from
173// old_inv + sorted-by-pair list of (bond_order, neighbor_old_inv).
174// Returns count of distinct invariants in new_inv (not strictly used
175// here but useful for diagnostic + future convergence-based variants).
176// =================================================================
177func nx_chem_fp_refine_step(m: *MolGraph, old_inv: *nx_int, new_inv: *nx_int) -> nx_int {
178 let pair_bonds: *nx_int = (sys_mmap(128) as *nx_int)
179 let pair_nbrs: *nx_int = (sys_mmap(128) as *nx_int)
180 let pair_keys: *nx_int = (sys_mmap(128) as *nx_int)
181 var i: nx_int = 0
182 while i < m.n_atoms {
183 var n_pairs: nx_int = 0
184 var bi: nx_int = 0
185 while bi < m.n_bonds {
186 let b: *Bond = ((m.bonds as nx_int) + (bi * NX_BOND_BYTES)) as *Bond
187 var other: nx_int = -1
188 if b.a == i { other = b.b }
189 if b.b == i { other = b.a }
190 if other >= 0 {
191 if n_pairs < 16 {
192 // pack (bond_order << 32) | neighbor_inv into pair_keys for sorting,
193 // and remember separate arrays for the actual fingerprint hash.
194 let key_packed: nx_int = (b.order << 32) | (old_inv[other] & 0xffffffff)
195 pair_keys[n_pairs] = key_packed
196 n_pairs = n_pairs + 1
197 }
198 }
199 bi = bi + 1
200 }
201 let _s: nx_int = nx_chem_fp_sort_asc(pair_keys, n_pairs)
202 // Unpack into bond + nbr arrays for hash function
203 var k: nx_int = 0
204 while k < n_pairs {
205 pair_bonds[k] = (pair_keys[k] >> 32) & 0xff
206 pair_nbrs[k] = pair_keys[k] & 0xffffffff
207 k = k + 1
208 }
209 new_inv[i] = nx_chem_fp_hash_neighborhood(old_inv[i], pair_bonds, pair_nbrs, n_pairs)
210 i = i + 1
211 }
212 return 0
213}
214
215// =================================================================
216// Set a bit in the bit-vector storage (slots of u64).
217// =================================================================
218func nx_chem_fp_set_bit(slots: *nx_int, bit_idx: nx_int) -> nx_int {
219 let slot: nx_int = bit_idx / NX_CHEM_FP_BITS_PER_SLOT
220 let offset: nx_int = bit_idx - (slot * NX_CHEM_FP_BITS_PER_SLOT)
221 slots[slot] = slots[slot] | (1 << offset)
222 return 0
223}
224
225// =================================================================
226// Test a bit in the bit-vector. Returns 0 or 1.
227// =================================================================
228func nx_chem_fp_get_bit(slots: *nx_int, bit_idx: nx_int) -> nx_int {
229 let slot: nx_int = bit_idx / NX_CHEM_FP_BITS_PER_SLOT
230 let offset: nx_int = bit_idx - (slot * NX_CHEM_FP_BITS_PER_SLOT)
231 return (slots[slot] >> offset) & 1
232}
233
234// =================================================================
235// Population count (number of 1-bits) of a u64 value. Naive shift+test
236// loop; the substrate uses i64 throughout so this is portable. For
237// large fingerprints with many slots, a table-driven popcount could
238// be substituted, but for 1024-bit fingerprints over 16 slots, the
239// loop is fine.
240// =================================================================
241func nx_chem_fp_popcount_u64(v: nx_int) -> nx_int {
242 var count: nx_int = 0
243 var x: nx_int = v
244 var i: nx_int = 0
245 while i < 64 {
246 if (x & 1) == 1 { count = count + 1 }
247 x = x >> 1
248 i = i + 1
249 }
250 return count
251}
252
253// =================================================================
254// Main entry: compute ECFP fingerprint for molecule m into a bit
255// vector of n_bits length (must be a power of 2 and a multiple of 64).
256//
257// out_slots is a caller-provided array of (n_bits / 64) i64 slots.
258// =================================================================
259func nx_chem_compute_ecfp_fingerprint(
260 m: *MolGraph, radius: nx_int, n_bits: nx_int, out_slots: *nx_int
261) -> nx_int {
262 let n_slots: nx_int = n_bits / NX_CHEM_FP_BITS_PER_SLOT
263 var zi: nx_int = 0
264 while zi < n_slots {
265 out_slots[zi] = 0
266 zi = zi + 1
267 }
268 if m.n_atoms == 0 { return 0 }
269 let inv_a: *nx_int = (sys_mmap((m.n_atoms * 8) as i64)) as *nx_int
270 let inv_b: *nx_int = (sys_mmap((m.n_atoms * 8) as i64)) as *nx_int
271 // Initial invariants
272 var i: nx_int = 0
273 while i < m.n_atoms {
274 inv_a[i] = nx_chem_fp_initial_invariant(m, i)
275 // Each radius-0 invariant contributes to the fingerprint
276 let bit_idx_0: nx_int = inv_a[i] & (n_bits - 1)
277 let _s0: nx_int = nx_chem_fp_set_bit(out_slots, bit_idx_0)
278 i = i + 1
279 }
280 // Iterative refinement; each radius contributes more invariants
281 var r: nx_int = 0
282 while r < radius {
283 let _ref: nx_int = nx_chem_fp_refine_step(m, inv_a, inv_b)
284 var j: nx_int = 0
285 while j < m.n_atoms {
286 let bit_idx: nx_int = inv_b[j] & (n_bits - 1)
287 let _s: nx_int = nx_chem_fp_set_bit(out_slots, bit_idx)
288 inv_a[j] = inv_b[j]
289 j = j + 1
290 }
291 r = r + 1
292 }
293 return 0
294}
295
296// =================================================================
297// Tanimoto similarity (Jaccard for binary fingerprints).
298// Returns Q4 fixed-point (0..10000) representing |A ∩ B| / |A ∪ B|.
299//
300// Both fingerprints must be the same size (n_bits / 64 slots).
301// =================================================================
302func nx_chem_fp_tanimoto_q4(fp_a: *nx_int, fp_b: *nx_int, n_bits: nx_int) -> nx_int {
303 let n_slots: nx_int = n_bits / NX_CHEM_FP_BITS_PER_SLOT
304 var intersect_count: nx_int = 0
305 var union_count: nx_int = 0
306 var i: nx_int = 0
307 while i < n_slots {
308 let a: nx_int = fp_a[i]
309 let b: nx_int = fp_b[i]
310 intersect_count = intersect_count + nx_chem_fp_popcount_u64(a & b)
311 union_count = union_count + nx_chem_fp_popcount_u64(a | b)
312 i = i + 1
313 }
314 if union_count == 0 { return 0 } // both empty: define as 0 (matches RDKit edge case)
315 // tanimoto in Q4 (0..10000)
316 return (intersect_count * NX_MAGIC_10000) / union_count
317}
318
319// =================================================================
320// Convenience: total bit count of a fingerprint.
321// =================================================================
322func nx_chem_fp_popcount(fp: *nx_int, n_bits: nx_int) -> nx_int {
323 let n_slots: nx_int = n_bits / NX_CHEM_FP_BITS_PER_SLOT
324 var total: nx_int = 0
325 var i: nx_int = 0
326 while i < n_slots {
327 total = total + nx_chem_fp_popcount_u64(fp[i])
328 i = i + 1
329 }
330 return total
331}
332
333// =================================================================
334// Convenience: bit-equality check. Returns 1 if all slots equal, 0 otherwise.
335// =================================================================
336func nx_chem_fp_equal(fp_a: *nx_int, fp_b: *nx_int, n_bits: nx_int) -> nx_int {
337 let n_slots: nx_int = n_bits / NX_CHEM_FP_BITS_PER_SLOT
338 var i: nx_int = 0
339 while i < n_slots {
340 if fp_a[i] != fp_b[i] { return 0 }
341 i = i + 1
342 }
343 return 1
344}