IntGraph connected components, connected components for UnitRings

src/mol-data/util/hash-functions.ts

@@ -59,4 +59,12 @@ export function hashString(s: string) {
  */
 export function cantorPairing(a: number, b: number) {
     return (a + b) * (a + b + 1) / 2 + b;
+}
+
+/**
+ * A unique number for each sorted pair of integers
+ * Biggest representable pair is (67108863, 67108863) (limit imposed by Number.MAX_SAFE_INTEGER)
+ */
+export function sortedCantorPairing(a: number, b: number) {
+    return a < b ? cantorPairing(a, b) : cantorPairing(b, a);
 }

src/mol-math/graph/int-adjacency-graph.ts

@@ -6,6 +6,7 @@
  */
 
 import { arrayPickIndices } from 'mol-data/util';
+import { LinkedIndex } from 'mol-data/int';
 
 /**
  * Represent a graph using vertex adjacency list.
@@ -126,6 +127,13 @@ export namespace IntAdjacencyGraph {
             this.curB = ob;
         }
 
+        /** Builds property-less graph */
+        addAllEdges() {
+            for (let i = 0; i < this.edgeCount; i++) {
+                this.addNextEdge();
+            }
+        }
+
         assignProperty<T>(prop: { [i: number]: T }, value: T) {
             prop[this.curA] = value;
             prop[this.curB] = value;
@@ -152,6 +160,12 @@ export namespace IntAdjacencyGraph {
         }
     }
 
+    export function fromVertexPairs(vertexCount: number, xs: number[], ys: number[]) {
+        const graphBuilder = new IntAdjacencyGraph.EdgeBuilder(vertexCount, xs, ys);
+        graphBuilder.addAllEdges();
+        return graphBuilder.createGraph();
+    }
+
     export function induceByVertices<P extends IntAdjacencyGraph.EdgePropsBase>(graph: IntAdjacencyGraph<P>, vertexIndices: ArrayLike<number>): IntAdjacencyGraph<P> {
         const { b, offset, vertexCount, edgeProps } = graph;
         const vertexMap = new Int32Array(vertexCount);
@@ -192,6 +206,52 @@ export namespace IntAdjacencyGraph {
 
         return create(newOffsets, newA, newB, newEdgeCount, newEdgeProps);
     }
+
+    export function connectedComponents(graph: IntAdjacencyGraph): { componentCount: number, componentIndex: Int32Array } {
+        const vCount = graph.vertexCount;
+
+        if (vCount === 0) return { componentCount: 0, componentIndex: new Int32Array(0) };
+        if (graph.edgeCount === 0) {
+            const componentIndex = new Int32Array(vCount);
+            for (let i = 0, _i = vCount; i < _i; i++) {
+                componentIndex[i] = i;
+            }
+            return { componentCount: vCount, componentIndex };
+        }
+
+        const componentIndex = new Int32Array(vCount);
+        for (let i = 0, _i = vCount; i < _i; i++) componentIndex[i] = -1;
+        let currentComponent = 0;
+        componentIndex[0] = currentComponent;
+
+        const { offset, b: neighbor } = graph;
+        const stack = [0];
+        const list = LinkedIndex(vCount);
+        list.remove(0);
+
+        while (stack.length > 0) {
+            const v = stack.pop()!;
+            const cIdx = componentIndex[v];
+
+            for (let eI = offset[v], _eI = offset[v + 1]; eI < _eI; eI++) {
+                const n = neighbor[eI];
+                if (!list.has(n)) continue;
+                list.remove(n);
+                stack.push(n);
+                componentIndex[n] = cIdx;
+            }
+
+            // check if we visited all vertices.
+            // If not, create a new component and continue.
+            if (stack.length === 0 && list.head >= 0) {
+                stack.push(list.head);
+                componentIndex[list.head] = ++currentComponent;
+                list.remove(list.head);
+            }
+        }
+
+        return { componentCount: vCount, componentIndex };
+    }
 }
 
 /**

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

@@ -4,18 +4,33 @@
  * @author David Sehnal <david.sehnal@gmail.com>
  */
 
-import computeRings from './rings/compute'
+import { computeRings, getFingerprint, createIndex } from './rings/compute'
 import Unit from '../unit';
 import StructureElement from '../element';
 
+type UnitRing = ReadonlyArray<StructureElement.UnitIndex>
+
 interface UnitRings {
     /** Each ring is specified as an array of indices in Unit.elements. */
-    readonly all: ReadonlyArray<ReadonlyArray<StructureElement.UnitIndex>>,
-    readonly byFingerprint: ReadonlyMap<string, ReadonlyArray<number>>
+    readonly all: ReadonlyArray<UnitRing>,
+    readonly byFingerprint: ReadonlyMap<string, ReadonlyArray<UnitRings.Index>>,
+
+    readonly index: {
+        /** Maps atom index inside a Unit to the smallest ring index (an atom can be part of more than one ring) */
+        readonly elementRingIndices: ReadonlyMap<StructureElement.UnitIndex, UnitRings.Index[]>,
+
+        /** Maps UnitRings.Index to index to ringComponents */
+        readonly ringComponentIndex: ReadonlyArray<UnitRings.ComponentIndex>,
+        readonly ringComponents: ReadonlyArray<ReadonlyArray<UnitRings.Index>>
+    }
 }
 
 namespace UnitRings {
-    export function getRingFingerprint(unit: Unit.Atomic, ring: ArrayLike<number>) {
+    /** Index into UnitRings.all */
+    export type Index = { readonly '@type': 'unit-ring-index' } & number
+    export type ComponentIndex = { readonly '@type': 'unit-ring-component-index' } & number
+
+    export function getRingFingerprint(unit: Unit.Atomic, ring: UnitRing) {
         const { elements } = unit;
         const { type_symbol } = unit.model.atomicHierarchy.atoms;
 
@@ -26,9 +41,9 @@ namespace UnitRings {
 
     export function create(unit: Unit.Atomic): UnitRings {
         const rings = computeRings(unit);
-        const byFingerprint = new Map<string, number[]>();
+        const byFingerprint = new Map<string, Index[]>();
 
-        let idx = 0;
+        let idx = 0 as Index;
         for (const r of rings) {
             const fp = getRingFingerprint(unit, r);
             if (byFingerprint.has(fp)) byFingerprint.get(fp)!.push(idx);
@@ -36,73 +51,19 @@ namespace UnitRings {
             idx++;
         }
 
-        return { all: rings, byFingerprint };
+        console.log(createIndex(rings));
+
+        let _index: UnitRings['index'] | undefined = void 0;
+        return {
+            all: rings,
+            byFingerprint,
+            get index() {
+                if (_index) return _index;
+                _index = createIndex(rings);
+                return _index;
+            }
+        };
     }
 }
 
-export { UnitRings }
-
-function getFingerprint(elements: string[]) {
-    const len = elements.length;
-    const reversed: string[] = new Array(len);
-
-    for (let i = 0; i < len; i++) reversed[i] = elements[len - i - 1];
-
-    const rotNormal = getMinimalRotation(elements);
-    const rotReversed = getMinimalRotation(reversed);
-
-    let isNormalSmaller = false;
-
-    for (let i = 0; i < len; i++) {
-        const u = elements[(i + rotNormal) % len], v = reversed[(i + rotReversed) % len];
-        if (u !== v) {
-            isNormalSmaller = u < v;
-            break;
-        }
-    }
-
-    if (isNormalSmaller) return buildFinderprint(elements, rotNormal);
-    return buildFinderprint(reversed, rotReversed);
-}
-
-function getMinimalRotation(elements: string[]) {
-    // adapted from http://en.wikipedia.org/wiki/Lexicographically_minimal_string_rotation
-
-    const len = elements.length;
-    const f = new Int32Array(len * 2);
-    for (let i = 0; i < f.length; i++) f[i] = -1;
-
-    let u = '', v = '', k = 0;
-
-    for (let j = 1; j < f.length; j++) {
-        let i = f[j - k - 1];
-        while (i !== -1) {
-            u = elements[j % len]; v = elements[(k + i + 1) % len];
-            if (u === v) break;
-            if (u < v) k = j - i - 1;
-            i = f[i];
-        }
-
-        if (i === -1) {
-            u = elements[j % len]; v = elements[(k + i + 1) % len];
-            if (u !== v) {
-                if (u < v) k = j;
-                f[j - k] = -1;
-            } else f[j - k] = i + 1;
-        } else f[j - k] = i + 1;
-    }
-
-    return k;
-}
-
-function buildFinderprint(elements: string[], offset: number) {
-    const len = elements.length;
-    const ret: string[] = [];
-    let i;
-    for (i = 0; i < len - 1; i++) {
-        ret.push(elements[(i + offset) % len]);
-        ret.push('-');
-    }
-    ret.push(elements[(i + offset) % len]);
-    return ret.join('');
-}
+export { UnitRing, UnitRings }

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

@@ -9,8 +9,10 @@ import { IntraUnitLinks } from '../links/data';
 import { Segmentation } from 'mol-data/int';
 import { LinkType } from '../../../model/types';
 import { StructureElement } from '../../../structure';
+import { sortedCantorPairing } from 'mol-data/util';
+import { IntAdjacencyGraph } from 'mol-math/graph';
 
-export default function computeRings(unit: Unit.Atomic) {
+export function computeRings(unit: Unit.Atomic) {
     const size = largestResidue(unit);
     const state = State(unit, size);
 
@@ -176,4 +178,127 @@ function findRings(state: State, from: number) {
             pred[other] = top;
         }
     }
+}
+
+export function getFingerprint(elements: string[]) {
+    const len = elements.length;
+    const reversed: string[] = new Array(len);
+
+    for (let i = 0; i < len; i++) reversed[i] = elements[len - i - 1];
+
+    const rotNormal = getMinimalRotation(elements);
+    const rotReversed = getMinimalRotation(reversed);
+
+    let isNormalSmaller = false;
+
+    for (let i = 0; i < len; i++) {
+        const u = elements[(i + rotNormal) % len], v = reversed[(i + rotReversed) % len];
+        if (u !== v) {
+            isNormalSmaller = u < v;
+            break;
+        }
+    }
+
+    if (isNormalSmaller) return buildFinderprint(elements, rotNormal);
+    return buildFinderprint(reversed, rotReversed);
+}
+
+function getMinimalRotation(elements: string[]) {
+    // adapted from http://en.wikipedia.org/wiki/Lexicographically_minimal_string_rotation
+
+    const len = elements.length;
+    const f = new Int32Array(len * 2);
+    for (let i = 0; i < f.length; i++) f[i] = -1;
+
+    let u = '', v = '', k = 0;
+
+    for (let j = 1; j < f.length; j++) {
+        let i = f[j - k - 1];
+        while (i !== -1) {
+            u = elements[j % len]; v = elements[(k + i + 1) % len];
+            if (u === v) break;
+            if (u < v) k = j - i - 1;
+            i = f[i];
+        }
+
+        if (i === -1) {
+            u = elements[j % len]; v = elements[(k + i + 1) % len];
+            if (u !== v) {
+                if (u < v) k = j;
+                f[j - k] = -1;
+            } else f[j - k] = i + 1;
+        } else f[j - k] = i + 1;
+    }
+
+    return k;
+}
+
+function buildFinderprint(elements: string[], offset: number) {
+    const len = elements.length;
+    const ret: string[] = [];
+    let i;
+    for (i = 0; i < len - 1; i++) {
+        ret.push(elements[(i + offset) % len]);
+        ret.push('-');
+    }
+    ret.push(elements[(i + offset) % len]);
+    return ret.join('');
+}
+
+type RingIndex = import('../rings').UnitRings.Index
+type RingComponentIndex = import('../rings').UnitRings.ComponentIndex
+
+export function createIndex(rings: StructureElement.UnitIndex[][]) {
+    const elementRingIndices: Map<StructureElement.UnitIndex, RingIndex[]> = new Map();
+
+    // for each ring atom, assign all rings that it is present in
+    for (let rI = 0 as RingIndex, _rI = rings.length; rI < _rI; rI++) {
+        const r = rings[rI];
+        for (let i = 0, _i = r.length; i < _i; i++) {
+            const e = r[i];
+            if (elementRingIndices.has(e)) elementRingIndices.get(e)!.push(rI);
+            else elementRingIndices.set(e, [rI]);
+        }
+    }
+
+    // create a graph where vertices are rings, edge if two rings share at least one atom
+    const addedEdges = new Set<number>();
+    const xs: RingIndex[] = [], ys: RingIndex[] = [];
+
+    for (let rI = 0 as RingIndex, _rI = rings.length; rI < _rI; rI++) {
+        const r = rings[rI];
+
+        for (let i = 0, _i = r.length; i < _i; i++) {
+            const e = r[i];
+
+            const containedRings = elementRingIndices.get(e)!;
+
+            if (containedRings.length === 1) continue;
+
+            for (let j = 0, _j = containedRings.length; j < _j; j++) {
+                const rJ = containedRings[j];
+                if (rI >= rJ) continue;
+
+                const edgeIndex = sortedCantorPairing(rI, rJ);
+                if (addedEdges.has(edgeIndex)) continue;
+
+                addedEdges.add(edgeIndex);
+                xs[xs.length] = rI;
+                ys[ys.length] = rJ;
+            }
+        }
+    }
+
+    const graph = IntAdjacencyGraph.fromVertexPairs(rings.length, xs, ys);
+    const components = IntAdjacencyGraph.connectedComponents(graph);
+
+    const ringComponentIndex = components.componentIndex as any as RingComponentIndex[];
+    const ringComponents: RingIndex[][] = [];
+    for (let i = 0; i < components.componentCount; i++) ringComponents[i] = [];
+
+    for (let rI = 0 as RingIndex, _rI = rings.length; rI < _rI; rI++) {
+        ringComponents[ringComponentIndex[rI]].push(rI);
+    }
+
+    return { elementRingIndices, ringComponentIndex, ringComponents };
 }