cdef enum:

THIELE_PI_EMPTY = 0 # an empty p orbital: C+, neutral B

THIELE_PI_ONE = 1 # one electron, from a ring double bond

THIELE_PI_PAIR = 2 # a lone pair: pyrrole N, furan O, thiophene S, C-, B-

cdef enum:

THIELE_INVALID = 3

cdef enum:

THIELE_CYCLE_MAX = 7

cdef uint64_t THIELE_CYCLE_BUDGET = 2000000

cdef enum:

THIELE_WHY_NONE = 0

THIELE_WHY_RADICAL = 1

THIELE_WHY_H_UNKNOWN = 2

cdef struct thiele_t:

void *block

uint8_t *ring_ok # [nrings] 1 while the ring is still a candidate

uint8_t *he_arom # [nhalf] 1 when this half-edge is in the candidate edge set

uint8_t *he_share # [nhalf] candidate rings holding this half-edge, saturating at 2

uint8_t *cand # [n] 1 when the atom is in a candidate ring

uint8_t *cls # [n] AROM_MUST / AROM_MAY / AROM_MUST_NOT / THIELE_INVALID

int32_t *comp # [n] component label over the candidate edge set, -1 outside

uint32_t *clist # [n] one component's atoms

uint32_t *stack # [n] breadth-first queue

cdef int arom_thiele_alloc(thiele_t *t, uint32_t n, uint32_t nhalf, uint32_t nrings) except -1:

# one struct, one malloc, one check, one free -- RULES.md 5.2

cdef size_t n_u8 = align8(<size_t> n * sizeof(uint8_t))

cdef size_t n_u32 = align8(<size_t> n * sizeof(uint32_t))

cdef size_t n_i32 = align8(<size_t> n * sizeof(int32_t))

# `or 1`: a molecule with no rings and a molecule with no bonds both reach here, and a zero-byte

# PyMem_Malloc may return NULL, which this function would report as MemoryError

cdef size_t half_u8 = align8(<size_t> (nhalf if nhalf else 1) * sizeof(uint8_t))

cdef size_t ring_u8 = align8(<size_t> (nrings if nrings else 1) * sizeof(uint8_t))

cdef size_t total = ring_u8 + 2 * half_u8 + 2 * n_u8 + n_i32 + 2 * n_u32

cdef char *block = <char *> PyMem_Malloc(total if total else 1)

if block is NULL:

raise MemoryError('aromatisation scratch allocation failed')

memset(block, 0, total)

t.block = <void *> block

cdef size_t off = 0

t.ring_ok = <uint8_t *> (block + off); off += ring_u8

t.he_arom = <uint8_t *> (block + off); off += half_u8

t.he_share = <uint8_t *> (block + off); off += half_u8

t.cand = <uint8_t *> (block + off); off += n_u8

t.cls = <uint8_t *> (block + off); off += n_u8

t.comp = <int32_t *> (block + off); off += n_i32

t.clist = <uint32_t *> (block + off); off += n_u32

t.stack = <uint32_t *> (block + off); off += n_u32

return 0

cdef class _ThieleRun:

self.rings = NULL

cdef list arom_thiele_short_cycles(Structure structure, uint64_t budget):

"""Every simple cycle of 3 to `THIELE_CYCLE_MAX` atoms, as index tuples in ring order, or `None`.

cdef uint32_t n = structure.header.atom_count

cdef uint32_t *ptr = csr_ptr(structure)

cdef halfedge_t *edges = csr_edges(structure)

cdef uint32_t path[THIELE_CYCLE_MAX]

cdef uint32_t slot[THIELE_CYCLE_MAX]

cdef uint8_t *onpath = <uint8_t *> PyMem_Malloc(n if n else 1)

if onpath is NULL:

raise MemoryError('aromatisation cycle scratch allocation failed')

memset(onpath, 0, n)

cdef list out = []

cdef uint64_t steps = 0

cdef uint32_t s, v, w, k, j

cdef int depth

cdef list cycle

try:

for s in range(n):

path[0] = s

slot[0] = ptr[s]

onpath[s] = 1

depth = 0

while depth >= 0:

v = path[depth]

if slot[depth] == ptr[v + 1]:

onpath[v] = 0

depth -= 1

continue

k = slot[depth]

slot[depth] += 1

if not edges[k].flags & HE_IN_RING:

continue

steps += 1

if steps > budget:

return None

w = edges[k].to

if w == s:

if depth >= 2 and path[1] < path[depth]:

cycle = []

for j in range(<uint32_t> depth + 1):

cycle.append(path[j])

out.append(tuple(cycle))

continue

if w < s or onpath[w] or depth + 2 > THIELE_CYCLE_MAX:

continue

depth += 1

path[depth] = w

slot[depth] = ptr[w]

onpath[w] = 1

finally:

PyMem_Free(onpath)

return out

cdef list arom_thiele_relevant(list cycles):

"""The relevant ones among `cycles`: those no set of strictly shorter cycles sums to over GF(2).

cdef dict bit = {}

cdef dict pivots = {}

cdef list relevant = []

cdef list keyed = []

cdef list sized = []

cdef list pending

cdef Py_ssize_t i = 0, j, k, size, count = len(cycles)

cdef object vec, key, top

cdef tuple cycle

for cycle in cycles:

keyed.append((len(cycle), sorted(cycle), cycle))

keyed.sort()

for key in keyed:

sized.append(key[2])

while i < count:

size = len(<tuple> sized[i])

j = i

pending = []

# reduced against the shorter cycles only: two cycles of one size are both relevant when

# neither shorter set spans them, even if they span each other

while j < count and len(<tuple> sized[j]) == size:

cycle = <tuple> sized[j]

vec = 0

for k in range(size):

key = (cycle[k], cycle[k + 1 if k + 1 < size else 0])

if key[0] > key[1]:

key = (key[1], key[0])

if key not in bit:

bit[key] = len(bit)

vec ^= 1 << bit[key]

while vec:

top = vec.bit_length() - 1

if top not in pivots:

break

vec ^= pivots[top]

if vec:

relevant.append(cycle)

pending.append(vec)

j += 1

for vec in pending:

while vec:

top = vec.bit_length() - 1

if top not in pivots:

pivots[top] = vec

break

vec ^= pivots[top]

i = j

return relevant

cdef uint32_t *arom_thiele_pack(list cycles) except NULL:

"""`cycles` in the `structure_rings` layout: count, `count + 1` offsets, then the atoms."""

cdef uint32_t nrings = <uint32_t> len(cycles)

cdef size_t total = 2 + nrings

cdef tuple cycle

for cycle in cycles:

total += len(cycle)

cdef uint32_t *out = <uint32_t *> PyMem_Malloc(total * sizeof(uint32_t))

if out is NULL:

raise MemoryError('aromatisation ring allocation failed')

cdef uint32_t i, k, off = 0, base = 2 + nrings

out[0] = nrings

for i in range(nrings):

cycle = <tuple> cycles[i]

out[1 + i] = off

for k in range(<uint32_t> len(cycle)):

out[base + off + k] = <uint32_t> cycle[k]

off += <uint32_t> len(cycle)

out[1 + nrings] = off

return out

cdef class ThieleResult:

f'refused={self.refused!r})')

cdef ThieleResult arom_thiele_result(MoleculeContainer mol, bint changed, list log, list refused):

cdef ThieleResult r = ThieleResult.__new__(ThieleResult)

r.changed = changed

r.log = log

r.refused = refused

# Unconditional, for the reason `arom_result` states.

mol.log.absorb('thiele', log, rule='thiele')

return r

cdef inline uint8_t arom_thiele_pi(atom_t *a, bint carries_double) noexcept nogil:

if carries_double:

return THIELE_PI_ONE

if a.charge > 0:

return THIELE_PI_EMPTY

if a.charge < 0:

return THIELE_PI_PAIR

if a.element == AROM_Z_B:

return THIELE_PI_EMPTY

return THIELE_PI_PAIR

cdef inline void arom_thiele_bonds(Structure structure, uint32_t i, uint32_t *nbrs,

uint32_t *doubles, uint32_t *triples,

uint32_t *aromatics) noexcept nogil:

cdef uint32_t *ptr = csr_ptr(structure)

cdef halfedge_t *edges = csr_edges(structure)

cdef uint32_t k

cdef uint8_t order

nbrs[0] = 0

doubles[0] = 0

triples[0] = 0

aromatics[0] = 0

for k in range(ptr[i], ptr[i + 1]):

order = edges[k].order

if order == 8:

continue # a dative bond is not a neighbour for this census

nbrs[0] += 1

if order == 2:

doubles[0] += 1

elif order == 3:

triples[0] += 1

elif order == 4:

aromatics[0] += 1

cdef void arom_thiele_share(Structure structure, thiele_t *t, uint32_t *rings,

uint32_t nrings) noexcept nogil:

cdef uint32_t base = 2 + nrings

cdef uint32_t i, k, u, v, size, slot

cdef uint32_t *ptr = csr_ptr(structure)

cdef halfedge_t *edges = csr_edges(structure)

memset(t.he_share, 0, ptr[structure.header.atom_count] * sizeof(uint8_t))

for i in range(nrings):

size = rings[2 + i] - rings[1 + i]

for k in range(size):

u = rings[base + rings[1 + i] + k]

v = rings[base + rings[1 + i] + (k + 1 if k + 1 < size else 0)]

for slot in range(ptr[u], ptr[u + 1]):

if edges[slot].to == v and t.he_share[slot] < 2:

t.he_share[slot] += 1

for slot in range(ptr[v], ptr[v + 1]):

if edges[slot].to == u and t.he_share[slot] < 2:

t.he_share[slot] += 1

cdef bint arom_thiele_ring_ok(Structure structure, uint32_t *ring, uint32_t size,

const uint8_t *he_share, uint8_t *cls_out, uint8_t *why,

uint32_t *culprit) noexcept nogil:

cdef uint32_t i, idx, k, nxt

cdef uint32_t *ptr = csr_ptr(structure)

cdef halfedge_t *edges = csr_edges(structure)

cdef uint32_t pi = 0

cdef bint carrier = False # some ring atom carries a double bond, wherever it points

cdef bint donor[7] # `size` is 5..7 by the gate below; no double bond here

# initialised at the declaration although `arom_thiele_bonds` fills all four: it fills them

# through pointers, which Cython cannot see, and the module compiles with warn.undeclared

cdef uint32_t nbrs = 0, doubles = 0, triples = 0, aromatics = 0

cdef atom_t *atoms = structure.atoms()

cdef atom_t *a

cdef uint8_t invalid

why[0] = THIELE_WHY_NONE

if size < 5 or size > 7:

return False

for i in range(size):

idx = ring[i]

a = atoms + idx

if a.element != AROM_Z_B and a.element != AROM_Z_C and a.element != AROM_Z_N \

and a.element != AROM_Z_O and a.element != AROM_Z_P and a.element != AROM_Z_S \

and a.element != AROM_Z_AS and a.element != AROM_Z_SE and a.element != AROM_Z_TE:

return False

if at_radical(a):

# a radical brings one electron and not two, so `arom_thiele_pi`'s charge rule does not

# describe it; refused rather than guessed

why[0] = THIELE_WHY_RADICAL

culprit[0] = idx

return False

if at_implicit_h_unknown(a):

# the pyrrole / pyridine choice IS the hydrogen count, so an unstated one cannot be

# resolved into a stated aromatic form. Note 4: no token on a prediction.

why[0] = THIELE_WHY_H_UNKNOWN

culprit[0] = idx

return False

arom_thiele_bonds(structure, idx, &nbrs, &doubles, &triples, &aromatics)

if aromatics:

return False # already aromatic here; not a Kekule ring

if nbrs > 3:

return False # over-coordinated for a ring atom, or a spiro atom

donor[i] = doubles == 0

if doubles == 1 and not triples:

continue

invalid = 0

cls_out[idx] = arom_classify_atom(a.element, a.charge, False, nbrs, False,

<int> at_implicit_h(a), &invalid)

if invalid or cls_out[idx] == AROM_MUST:

return False

# no donor is flanked by two donors; see the docstring's paragraph on the phosphite bridge

for i in range(size):

if donor[i] and donor[size - 1 if i == 0 else i - 1] and donor[0 if i + 1 == size else i + 1]:

return False

# the ring holds pi of its own, or borrows a sextet; see the docstring's last two paragraphs

for i in range(size):

idx = ring[i]

nxt = ring[i + 1 if i + 1 < size else 0]

pi += arom_thiele_pi(atoms + idx, not donor[i])

if not donor[i]:

carrier = True

for k in range(ptr[idx], ptr[idx + 1]):

if edges[k].to == nxt:

if edges[k].order == 2 and he_share[k] < 2:

return True # a double bond of the ring's own: nothing is borrowed

break

return carrier and pi % 4 == 2

cdef void arom_thiele_mark(Structure structure, thiele_t *t, uint32_t *rings,

uint32_t nrings) noexcept nogil:

cdef uint32_t base = 2 + nrings

cdef uint32_t i, k, u, v, size, slot

cdef uint32_t *ptr = csr_ptr(structure)

cdef halfedge_t *edges = csr_edges(structure)

memset(t.cand, 0, structure.header.atom_count * sizeof(uint8_t))

memset(t.he_arom, 0, ptr[structure.header.atom_count] * sizeof(uint8_t))

for i in range(nrings):

if not t.ring_ok[i]:

continue

size = rings[2 + i] - rings[1 + i]

for k in range(size):

u = rings[base + rings[1 + i] + k]

v = rings[base + rings[1 + i] + (k + 1 if k + 1 < size else 0)]

t.cand[u] = 1

for slot in range(ptr[u], ptr[u + 1]):

if edges[slot].to == v:

t.he_arom[slot] = 1

for slot in range(ptr[v], ptr[v + 1]):

if edges[slot].to == u:

t.he_arom[slot] = 1

cdef bint arom_thiele_exo(Structure structure, thiele_t *t, uint32_t i) noexcept nogil:

"""True when atom `i` carries a double or triple bond outside the candidate edge set."""

cdef uint32_t *ptr = csr_ptr(structure)

cdef halfedge_t *edges = csr_edges(structure)

cdef uint32_t k

cdef uint8_t order

for k in range(ptr[i], ptr[i + 1]):

order = edges[k].order

if (order == 2 or order == 3) and not t.he_arom[k]:

return True

return False

cdef uint32_t arom_thiele_prune(Structure structure, thiele_t *t, uint32_t *rings,

uint32_t nrings) noexcept nogil:

cdef uint32_t base = 2 + nrings

cdef uint32_t i, k, size, dropped = 0

cdef bint again = True

while again:

again = False

for i in range(nrings):

if not t.ring_ok[i]:

continue

size = rings[2 + i] - rings[1 + i]

for k in range(size):

if arom_thiele_exo(structure, t, rings[base + rings[1 + i] + k]):

t.ring_ok[i] = 0

dropped += 1

again = True

break

if again:

arom_thiele_mark(structure, t, rings, nrings)

return dropped

cdef uint32_t arom_thiele_component(Structure structure, thiele_t *t, uint32_t seed,

uint32_t *nedges) noexcept nogil:

"""One connected component of the candidate edge set, breadth first. Atoms into `t.clist`."""

cdef uint32_t *ptr = csr_ptr(structure)

cdef halfedge_t *edges = csr_edges(structure)

cdef uint32_t head = 0, count = 1, v, w, k

t.clist[0] = seed

t.comp[seed] = <int32_t> seed

nedges[0] = 0

while head < count:

v = t.clist[head]

head += 1

for k in range(ptr[v], ptr[v + 1]):

if not t.he_arom[k]:

continue

w = edges[k].to

# each candidate bond is two half-edges; the v < w test counts it once

if v < w:

nedges[0] += 1

if t.comp[w] < 0:

t.comp[w] = <int32_t> seed

t.clist[count] = w

count += 1

return count

cdef bint arom_thiele_check(Structure structure, thiele_t *t, uint32_t count,

uint32_t nedges, uint32_t *pi_out) noexcept nogil:

cdef uint32_t *ptr = csr_ptr(structure)

cdef halfedge_t *edges = csr_edges(structure)

cdef uint32_t j, i, k, deg, dbl

cdef atom_t *atoms = structure.atoms()

cdef bint cycle = nedges == count

pi_out[0] = 0

for j in range(count):

i = t.clist[j]

if t.cls[i] == THIELE_INVALID:

return False

deg = 0

dbl = 0

for k in range(ptr[i], ptr[i + 1]):

if not t.he_arom[k]:

continue

deg += 1

if edges[k].order == 2:

dbl += 1

if deg != 2:

cycle = False

if t.cls[i] == AROM_MUST:

if dbl != 1:

return False

elif t.cls[i] == AROM_MUST_NOT:

if dbl:

# `kekule()` would never have put a ring double bond on a lone-pair donor, so these

# orders are not a Kekule form of this edge set

return False

elif dbl > 1:

return False # MAY takes one or none, never two

pi_out[0] += arom_thiele_pi(atoms + i, dbl == 1)

if cycle and pi_out[0] % 4 != 2:

return False

return True

def thiele(MoleculeContainer mol not None):

return arom_thiele_result(mol, changed, log, refused)