ring computation algorithm fix

CHANGELOG.md

@@ -11,6 +11,7 @@ Note that since we don't clearly distinguish between a public and private interf
 - Fix "texture not renderable" & "texture not bound" warnings (#319)
 - Fix visual for bonds between two aromatic rings
 - Fix visual for delocalized bonds (parsed from mmcif and mol2)
+- Fix ring computation algorithm
 
 ## [v3.2.0] - 2022-02-17
 

src/mol-model/structure/structure/unit/rings/compute.ts

@@ -36,9 +36,11 @@ interface State {
     endVertex: number,
     count: number,
     visited: Int32Array,
+    marked: Int32Array,
     queue: Int32Array,
     color: Int32Array,
     pred: Int32Array,
+    depth: Int32Array,
 
     left: Int32Array,
     right: Int32Array,
@@ -60,8 +62,10 @@ function State(unit: Unit.Atomic, capacity: number): State {
         endVertex: 0,
         count: 0,
         visited: new Int32Array(capacity),
+        marked: new Int32Array(capacity),
         queue: new Int32Array(capacity),
         pred: new Int32Array(capacity),
+        depth: new Int32Array(capacity),
         left: new Int32Array(Constants.MaxDepth),
         right: new Int32Array(Constants.MaxDepth),
         color: new Int32Array(capacity),
@@ -78,17 +82,31 @@ function State(unit: Unit.Atomic, capacity: number): State {
 
 function resetState(state: State) {
     state.count = state.endVertex - state.startVertex;
-    const { visited, pred, color } = state;
+    const { visited, pred, color, depth, marked } = state;
+    const { offset: bondOffset } = state.bonds;
     for (let i = 0; i < state.count; i++) {
-        visited[i] = -1;
+        const a = state.startVertex + i;
+        const bStart = bondOffset[a], bEnd = bondOffset[a + 1];
+        // Terminal/standalone atoms (< 2 bonds) can't be part of rings so
+        // mark them as visited
+        visited[i] = bEnd - bStart < 2 ? 1 : -1;
         pred[i] = -1;
+        marked[i] = -1;
         color[i] = 0;
+        depth[i] = 0;
     }
     state.currentColor = 0;
     state.currentAltLoc = '';
     state.hasAltLoc = false;
 }
 
+function resetDepth(state: State) {
+    const { depth } = state;
+    for (let i = 0; i < state.count; i++) {
+        depth[i] = state.count + 1;
+    }
+}
+
 function largestResidue(unit: Unit.Atomic) {
     const residuesIt = Segmentation.transientSegments(unit.model.atomicHierarchy.residueAtomSegments, unit.elements);
     let size = 0;
@@ -117,11 +135,13 @@ function processResidue(state: State, start: number, end: number) {
     }
     arraySetRemove(altLocs, '');
 
+    let mark = 1;
     if (altLocs.length === 0) {
         resetState(state);
         for (let i = 0; i < state.count; i++) {
             if (visited[i] >= 0) continue;
-            findRings(state, i);
+            resetDepth(state);
+            mark = findRings(state, i, mark);
         }
     } else {
         for (let aI = 0; aI < altLocs.length; aI++) {
@@ -134,7 +154,8 @@ function processResidue(state: State, start: number, end: number) {
                 if (altLoc && altLoc !== state.currentAltLoc) {
                     continue;
                 }
-                findRings(state, i);
+                resetDepth(state);
+                mark = findRings(state, i, mark);
             }
         }
     }
@@ -147,7 +168,7 @@ function processResidue(state: State, start: number, end: number) {
 function addRing(state: State, a: number, b: number) {
     // only "monotonous" rings
     if (b < a) {
-        return;
+        return false;
     }
 
     const { pred, color, left, right } = state;
@@ -176,7 +197,7 @@ function addRing(state: State, a: number, b: number) {
         if (current < 0) break;
     }
     if (!found) {
-        return;
+        return false;
     }
 
     current = a;
@@ -190,7 +211,7 @@ function addRing(state: State, a: number, b: number) {
     const len = leftOffset + rightOffset;
     // rings must have at least three elements
     if (len < 3) {
-        return;
+        return false;
     }
 
     const ring = new Int32Array(len);
@@ -200,43 +221,46 @@ function addRing(state: State, a: number, b: number) {
 
     sortArray(ring);
 
-    if (state.hasAltLoc) {
-        // we need to check if the ring was already added because alt locs are present.
-
-        for (let rI = 0, _rI = state.currentRings.length; rI < _rI; rI++) {
-            const r = state.currentRings[rI];
-            if (ring[0] !== r[0]) continue;
-            if (ring.length !== r.length) continue;
-
-            let areSame = true;
-            for (let aI = 0, _aI = ring.length; aI < _aI; aI++) {
-                if (ring[aI] !== r[aI]) {
-                    areSame = false;
-                    break;
-                }
-            }
-            if (areSame) {
-                return;
+    // Check if the ring is unique
+    for (let rI = 0, _rI = state.currentRings.length; rI < _rI; rI++) {
+        const r = state.currentRings[rI];
+        if (ring[0] !== r[0]) continue;
+        if (ring.length !== r.length) continue;
+
+        let areSame = true;
+        for (let aI = 0, _aI = ring.length; aI < _aI; aI++) {
+            if (ring[aI] !== r[aI]) {
+                areSame = false;
+                break;
             }
         }
+        if (areSame) {
+            return true;
+        }
     }
 
     state.currentRings.push(SortedArray.ofSortedArray(ring));
+
+    return true;
 }
 
-function findRings(state: State, from: number) {
-    const { bonds, startVertex, endVertex, visited, queue, pred } = state;
+function findRings(state: State, from: number, mark: number) {
+    const { bonds, startVertex, endVertex, visited, marked, queue, pred, depth } = state;
     const { elements } = state.unit;
     const { b: neighbor, edgeProps: { flags: bondFlags }, offset } = bonds;
-    visited[from] = 1;
+    marked[from] = mark;
+    depth[from] = 0;
     queue[0] = from;
     let head = 0, size = 1;
 
     while (head < size) {
         const top = queue[head++];
+        const d = depth[top];
         const a = startVertex + top;
         const start = offset[a], end = offset[a + 1];
 
+        visited[top] = 1;
+
         for (let i = start; i < end; i++) {
             const b = neighbor[i];
             if (b < startVertex || b >= endVertex || !BondType.isCovalent(bondFlags[i])) continue;
@@ -250,18 +274,26 @@ function findRings(state: State, from: number) {
 
             const other = b - startVertex;
 
-            if (visited[other] > 0) {
+            if (marked[other] === mark) {
                 if (pred[other] !== top && pred[top] !== other) {
-                    addRing(state, top, other);
+                    if (addRing(state, top, other)) {
+                        visited[other] = 1;
+                        return mark + 1;
+                    }
                 }
                 continue;
             }
 
-            visited[other] = 1;
+            const newDepth = Math.min(depth[other], d + 1);
+            if (newDepth > Constants.MaxDepth) continue;
+
+            depth[other] = newDepth;
+            marked[other] = mark;
             queue[size++] = other;
             pred[other] = top;
         }
     }
+    return mark + 1;
 }
 
 export function getFingerprint(elements: string[]) {