Browse Source

Merge branch 'master' into cartoon-repr

Alexander Rose 6 years ago
parent
commit
34cbd3fb56

+ 1 - 0
package.json

@@ -15,6 +15,7 @@
     "build": "cpx \"src/**/*.{vert,frag,glsl,scss,woff,woff2,ttf,otf,eot,svg,html}\" build/node_modules/ && tsc",
     "watch": "tsc -watch",
     "watch-extra": "cpx \"src/**/*.{vert,frag,glsl,scss,woff,woff2,ttf,otf,eot,svg,html}\" build/node_modules/ --watch",
+    "watch-all-win": "start cmd /K npm run watch & start cmd /K npm run watch-extra & start cmd /K npm run watch-viewer & start http-server -p 1338",
     "test": "jest",
     "build-viewer": "webpack build/node_modules/apps/viewer/index.js --mode development -o build/viewer/index.js",
     "watch-viewer": "webpack build/node_modules/apps/viewer/index.js -w --mode development -o build/viewer/index.js",

+ 2 - 3
src/apps/structure-info/model.ts

@@ -9,11 +9,10 @@ import * as argparse from 'argparse'
 require('util.promisify').shim();
 
 import { CifFrame } from 'mol-io/reader/cif'
-import { Model, Structure, StructureElement, Unit, Format, StructureProperties } from 'mol-model/structure'
+import { Model, Structure, StructureElement, Unit, Format, StructureProperties, UnitRing } from 'mol-model/structure'
 // import { Run, Progress } from 'mol-task'
 import { OrderedSet } from 'mol-data/int';
 import { openCif, downloadCif } from './helpers';
-import { UnitRings } from 'mol-model/structure/structure/unit/rings';
 import { Vec3 } from 'mol-math/linear-algebra';
 
 
@@ -136,7 +135,7 @@ export function printRings(structure: Structure) {
         const { all, byFingerprint } = unit.rings;
         const fps: string[] = [];
         for (let i = 0, _i = Math.min(5, all.length); i < _i; i++) {
-            fps[fps.length] = UnitRings.getRingFingerprint(unit, all[i]);
+            fps[fps.length] = UnitRing.fingerprint(unit, all[i]);
         }
         if (all.length > 5) fps.push('...')
         console.log(`Unit ${unit.id}, ${all.length} ring(s), ${byFingerprint.size} different fingerprint(s).\n  ${fps.join(', ')}`);

+ 8 - 0
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);
 }

+ 125 - 15
src/mol-math/graph/int-adjacency-graph.ts

@@ -5,7 +5,8 @@
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
 
-import { arrayPickIndices } from 'mol-data/util';
+import { arrayPickIndices, cantorPairing } from 'mol-data/util';
+import { LinkedIndex, SortedArray } 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,41 @@ export namespace IntAdjacencyGraph {
         }
     }
 
+    export class UniqueEdgeBuilder {
+        private xs: number[] = [];
+        private ys: number[] = [];
+        private included = new Set<number>();
+
+        addEdge(i: number, j: number) {
+            let u = i, v = j;
+            if (i > j) { u = j; v = i; }
+            const id = cantorPairing(u, v);
+            if (this.included.has(id)) return false;
+            this.included.add(id);
+            this.xs[this.xs.length] = u;
+            this.ys[this.ys.length] = v;
+            return true;
+        }
+
+        getGraph(): IntAdjacencyGraph {
+            return fromVertexPairs(this.vertexCount, this.xs, this.ys);
+        }
+
+        // if we cant to add custom props as well
+        getEdgeBuiler() {
+            return new EdgeBuilder(this.vertexCount, this.xs, this.ys);
+        }
+
+        constructor(public vertexCount: number) {
+        }
+    }
+
+    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,22 +235,89 @@ 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 };
+    }
+
+    /**
+     * Check if any vertex in `verticesA` is connected to any vertex in `verticesB`
+     * via at most `maxDistance` edges.
+     *
+     * Returns true if A and B are intersecting.
+     */
+    export function areVertexSetsConnected(graph: IntAdjacencyGraph, verticesA: SortedArray<number>, verticesB: SortedArray<number>, maxDistance: number): boolean {
+        // check if A and B are intersecting, this handles maxDepth = 0
+        if (SortedArray.areIntersecting(verticesA, verticesB)) return true;
+        if (maxDistance < 1) return false;
+
+        const visited = new Set<number>();
+        for (let i = 0, il = verticesA.length; i < il; ++i) {
+            visited.add(verticesA[i]);
+        }
+
+        return areVertexSetsConnectedImpl(graph, verticesA, verticesB, maxDistance, visited);
+    }
 }
 
-/**
- * Check if any vertex in `verticesA` is connected to any vertex in `verticesB`
- * via `depth` hops or intermediate vertices
- */
-export function areConnected(verticesA: ReadonlyArray<number>, verticesB: ReadonlyArray<number>, graph: IntAdjacencyGraph, depth: number): boolean {
-    const { b, offset } = graph
-    const linkedVectices: number[] = []
-    for (let i = 0, il = verticesA.length; i < il; ++i) {
-        const vi = verticesA[i]
-        for (let j = offset[vi], jl = offset[vi + 1]; j < jl; ++j) {
-            const li = b[j]
-            if (verticesB.includes(li)) return true
-            if (!verticesA.includes(li)) linkedVectices.push(li)
+function areVertexSetsConnectedImpl(graph: IntAdjacencyGraph, frontier: ArrayLike<number>, target: SortedArray<number>, distance: number, visited: Set<number>): boolean {
+    const { b: neighbor, offset } = graph;
+    const newFrontier: number[] = [];
+
+    for (let i = 0, il = frontier.length; i < il; ++i) {
+        const src = frontier[i];
+
+        for (let j = offset[src], jl = offset[src + 1]; j < jl; ++j) {
+            const other = neighbor[j];
+            if (visited.has(other)) continue;
+            if (SortedArray.has(target, other)) return true;
+
+            visited.add(other);
+            newFrontier[newFrontier.length] = other;
         }
     }
-    return depth > 0 ? areConnected(linkedVectices, verticesB, graph, depth - 1) : false
+
+    return distance > 1 ? areVertexSetsConnectedImpl(graph, newFrontier, target, distance - 1, visited) : false;
 }

+ 2 - 1
src/mol-model/structure/structure.ts

@@ -11,4 +11,5 @@ import StructureSymmetry from './structure/symmetry'
 import { Link } from './structure/unit/links'
 import StructureProperties from './structure/properties'
 
-export { StructureElement, Link, Structure, Unit, StructureSymmetry, StructureProperties }
+export { StructureElement, Link, Structure, Unit, StructureSymmetry, StructureProperties }
+export * from './structure/unit/rings'

+ 16 - 46
src/mol-model/structure/structure/carbohydrates/compute.ts

@@ -6,7 +6,7 @@
 
 import { Segmentation } from 'mol-data/int';
 import { combinations } from 'mol-data/util/combination';
-import { areConnected } from 'mol-math/graph';
+import { IntAdjacencyGraph } from 'mol-math/graph';
 import { Vec3 } from 'mol-math/linear-algebra';
 import PrincipalAxes from 'mol-math/linear-algebra/matrix/principal-axes';
 import { fillSerial } from 'mol-util/array';
@@ -18,43 +18,12 @@ import Structure from '../structure';
 import Unit from '../unit';
 import { SaccharideNameMap, UnknownSaccharideComponent } from './constants';
 import { CarbohydrateElement, CarbohydrateLink, Carbohydrates, CarbohydrateTerminalLink } from './data';
+import { UnitRings, UnitRing } from '../unit/rings';
 
-function getResidueIndex(elementIndex: number, unit: Unit.Atomic) {
-    return unit.model.atomicHierarchy.residueAtomSegments.index[unit.elements[elementIndex]]
-}
-
-function sugarResidueIdx(unit: Unit.Atomic, ring: ArrayLike<StructureElement.UnitIndex>): ResidueIndex {
-    const { elements } = unit;
-    const residueIndex = unit.model.atomicHierarchy.residueAtomSegments.index;
-    const idx = residueIndex[elements[ring[0]]];
-    for (let rI = 1, _rI = ring.length; rI < _rI; rI++) {
-        if (idx !== residueIndex[elements[ring[rI]]]) return -1 as ResidueIndex;
-    }
-    return idx;
-}
+const C = ElementSymbol('C'), O = ElementSymbol('O');
+const SugarRingFps = [UnitRing.elementFingerprint([C, C, C, C, C, O]), UnitRing.elementFingerprint([C, C, C, C, O])]
 
-function addSugarRings(unit: Unit.Atomic, fp: string, sugarResidues: Map<ResidueIndex, number[]>) {
-    const rings = unit.rings;
-    const byFp = rings.byFingerprint.get(fp);
-    if (!byFp) return;
-    for (const r of byFp) {
-        const idx = sugarResidueIdx(unit, rings.all[r]);
-        if (idx >= 0) {
-            if (sugarResidues.has(idx)) sugarResidues.get(idx)!.push(r);
-            else sugarResidues.set(idx, [r]);
-        }
-    }
-}
-
-function getSugarRingIndices(unit: Unit.Atomic) {
-    const sugarResidues = new Map<ResidueIndex, number[]>();
-    addSugarRings(unit, 'C-C-C-C-C-O', sugarResidues);
-    addSugarRings(unit, 'C-C-C-C-O', sugarResidues);
-    return sugarResidues;
-}
-
-const C = ElementSymbol('C')
-function getDirection(direction: Vec3, unit: Unit.Atomic, indices: ReadonlyArray<StructureElement.UnitIndex>, center: Vec3) {
+function getDirection(direction: Vec3, unit: Unit.Atomic, indices: ArrayLike<StructureElement.UnitIndex>, center: Vec3) {
     let indexC1 = -1, indexC1X = -1, indexC = -1
     const { elements } = unit
     const { position } = unit.conformation
@@ -73,8 +42,8 @@ function getDirection(direction: Vec3, unit: Unit.Atomic, indices: ReadonlyArray
     }
     const index = indexC1 !== -1 ? indexC1
         : indexC1X !== -1 ? indexC1X
-        : indexC !== -1 ? indexC
-        : elements[indices[0]]
+            : indexC !== -1 ? indexC
+                : elements[indices[0]]
     Vec3.normalize(direction, Vec3.sub(direction, center, position(index, direction)))
     return direction
 }
@@ -85,6 +54,7 @@ function getAtomId(unit: Unit.Atomic, index: number) {
     return label_atom_id.value(elements[index])
 }
 
+
 export function computeCarbohydrates(structure: Structure): Carbohydrates {
     const links: CarbohydrateLink[] = []
     const terminalLinks: CarbohydrateTerminalLink[] = []
@@ -120,7 +90,7 @@ export function computeCarbohydrates(structure: Structure): Carbohydrates {
         const chainIt = Segmentation.transientSegments(chainAtomSegments, unit.elements)
         const residueIt = Segmentation.transientSegments(residueAtomSegments, unit.elements)
 
-        let sugarResidueMap: Map<ResidueIndex, number[]> | undefined = void 0;
+        let sugarResidueMap: Map<ResidueIndex, UnitRings.Index[]> | undefined = void 0;
 
         while (chainIt.hasNext) {
             residueIt.setSegment(chainIt.move());
@@ -134,7 +104,7 @@ export function computeCarbohydrates(structure: Structure): Carbohydrates {
                 }
 
                 if (!sugarResidueMap) {
-                    sugarResidueMap = getSugarRingIndices(unit);
+                    sugarResidueMap = UnitRings.byFingerprintAndResidue(unit.rings, SugarRingFps);
                 }
 
                 const sugarRings = sugarResidueMap.get(residueIndex);
@@ -163,19 +133,19 @@ export function computeCarbohydrates(structure: Structure): Carbohydrates {
                     ringElements.push(elementIndex)
                     elementsWithRingMap.set(elementKey(residueIndex, unit.id), elementIndex)
                     elements.push({
-                        geometry: { center, normal, direction, },
+                        geometry: { center, normal, direction },
                         hasRing: true,
                         unit, residueIndex, component: saccharideComp
                     })
                 }
 
                 // add carbohydrate links induced by intra-residue bonds
-                const ringCombinations = combinations(fillSerial(new Array(sugarRings.length)), 2)
+                const ringCombinations = combinations(fillSerial(new Array(sugarRings.length) as number[]), 2)
                 for (let j = 0, jl = ringCombinations.length; j < jl; ++j) {
                     const rc = ringCombinations[j];
                     const r0 = rings.all[sugarRings[rc[0]]], r1 = rings.all[sugarRings[rc[1]]];
-                    if (areConnected(r0, r1, unit.links, 2)) {
-                        // fix both directions as it is unlcear where the C1 atom is
+                    if (IntAdjacencyGraph.areVertexSetsConnected(unit.links, r0, r1, 2)) {
+                        // fix both directions as it is unclear where the C1 atom is
                         fixLinkDirection(ringElements[rc[0]], ringElements[rc[1]])
                         fixLinkDirection(ringElements[rc[1]], ringElements[rc[0]])
                         links.push({
@@ -202,8 +172,8 @@ export function computeCarbohydrates(structure: Structure): Carbohydrates {
                 pairBonds.getBonds(indexA).forEach(bondInfo => {
                     const { unitA, unitB } = pairBonds
                     const indexB = bondInfo.indexB
-                    const elementIndexA = elementsWithRingMap.get(elementKey(getResidueIndex(indexA, unitA), unitA.id))
-                    const elementIndexB = elementsWithRingMap.get(elementKey(getResidueIndex(indexB, unitB), unitB.id))
+                    const elementIndexA = elementsWithRingMap.get(elementKey(unitA.getResidueIndex(indexA), unitA.id))
+                    const elementIndexB = elementsWithRingMap.get(elementKey(unitB.getResidueIndex(indexB), unitB.id))
 
                     if (elementIndexA !== undefined && elementIndexB !== undefined) {
                         const atomIdA = getAtomId(unitA, indexA)

+ 4 - 0
src/mol-model/structure/structure/unit.ts

@@ -117,6 +117,10 @@ namespace Unit {
             return this.props.rings.ref;
         }
 
+        getResidueIndex(elementIndex: StructureElement.UnitIndex) {
+            return this.model.atomicHierarchy.residueAtomSegments.index[this.elements[elementIndex]];
+        }
+
         constructor(id: number, invariantId: number, model: Model, elements: StructureElement.Set, conformation: SymmetryOperator.ArrayMapping, props: AtomicProperties) {
             this.id = id;
             this.invariantId = invariantId;

+ 3 - 2
src/mol-model/structure/structure/unit/links/data.ts

@@ -8,6 +8,7 @@
 import { LinkType } from '../../../model/types'
 import { IntAdjacencyGraph } from 'mol-math/graph';
 import Unit from '../../unit';
+import StructureElement from '../../element';
 
 type IntraUnitLinks = IntAdjacencyGraph<{ readonly order: ArrayLike<number>, readonly flags: ArrayLike<LinkType.Flag> }>
 
@@ -76,14 +77,14 @@ namespace InterUnitBonds {
         }
 
         constructor(public unitA: Unit.Atomic, public unitB: Unit.Atomic,
-            public bondCount: number, public linkedElementIndices: ReadonlyArray<number>,
+            public bondCount: number, public linkedElementIndices: ReadonlyArray<StructureElement.UnitIndex>,
             private linkMap: Map<number, BondInfo[]>) {
         }
     }
 
     export interface BondInfo {
         /** indexInto */
-        readonly indexB: number,
+        readonly indexB: StructureElement.UnitIndex,
         readonly order: number,
         readonly flag: LinkType.Flag
     }

+ 3 - 2
src/mol-model/structure/structure/unit/links/inter-compute.ts

@@ -13,6 +13,7 @@ import { InterUnitBonds } from './data';
 import { UniqueArray } from 'mol-data/generic';
 import { SortedArray } from 'mol-data/int';
 import { Vec3, Mat4 } from 'mol-math/linear-algebra';
+import StructureElement from '../../element';
 
 const MAX_RADIUS = 4;
 
@@ -24,8 +25,8 @@ function addMapEntry<A, B>(map: Map<A, B[]>, a: A, b: B) {
 interface PairState {
     mapAB: Map<number, InterUnitBonds.BondInfo[]>,
     mapBA: Map<number, InterUnitBonds.BondInfo[]>,
-    bondedA: UniqueArray<number, number>,
-    bondedB: UniqueArray<number, number>
+    bondedA: UniqueArray<StructureElement.UnitIndex, StructureElement.UnitIndex>,
+    bondedB: UniqueArray<StructureElement.UnitIndex, StructureElement.UnitIndex>
 }
 
 function addLink(indexA: number, indexB: number, order: number, flag: LinkType.Flag, state: PairState) {

+ 97 - 75
src/mol-model/structure/structure/unit/rings.ts

@@ -4,105 +4,127 @@
  * @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';
+import { SortedArray } from 'mol-data/int';
+import { ResidueIndex } from '../../model';
+import { ElementSymbol } from '../../model/types';
 
-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>>
-}
+type UnitRing = SortedArray<StructureElement.UnitIndex>
 
-namespace UnitRings {
-    export function getRingFingerprint(unit: Unit.Atomic, ring: ArrayLike<number>) {
-        const { elements } = unit;
-        const { type_symbol } = unit.model.atomicHierarchy.atoms;
-
-        const symbols: string[] = [];
-        for (let i = 0, _i = ring.length; i < _i; i++) symbols[symbols.length] = type_symbol.value(elements[ring[i]]) as String as string;
-        return getFingerprint(symbols);
+class UnitRings {
+    /** Each ring is specified as an array of indices in Unit.elements. */
+    readonly all: ReadonlyArray<UnitRing>;
+
+    private _byFingerprint?: ReadonlyMap<UnitRing.Fingerprint, ReadonlyArray<UnitRings.Index>>;
+    private _index?: {
+        readonly elementRingIndices: ReadonlyMap<StructureElement.UnitIndex, UnitRings.Index[]>,
+        readonly ringComponentIndex: ReadonlyArray<UnitRings.ComponentIndex>,
+        readonly ringComponents: ReadonlyArray<ReadonlyArray<UnitRings.Index>>
+    };
+
+    private get index() {
+        if (this._index) return this._index;
+        this._index = createIndex(this.all);
+        return this._index;
     }
 
-    export function create(unit: Unit.Atomic): UnitRings {
-        const rings = computeRings(unit);
-        const byFingerprint = new Map<string, number[]>();
-
-        let idx = 0;
-        for (const r of rings) {
-            const fp = getRingFingerprint(unit, r);
-            if (byFingerprint.has(fp)) byFingerprint.get(fp)!.push(idx);
-            else byFingerprint.set(fp, [idx]);
-            idx++;
-        }
+    get byFingerprint() {
+        if (this._byFingerprint) return this._byFingerprint;
+        this._byFingerprint = createByFingerprint(this.unit, this.all);
+        return this._byFingerprint;
+    }
 
-        return { all: rings, byFingerprint };
+    /** Maps atom index inside a Unit to the smallest ring index (an atom can be part of more than one ring) */
+    get elementRingIndices() {
+        return this.index.elementRingIndices;
     }
-}
 
-export { UnitRings }
+    /** Maps UnitRings.Index to index to ringComponents */
+    get ringComponentIndex() {
+        return this.index.ringComponentIndex;
+    }
 
-function getFingerprint(elements: string[]) {
-    const len = elements.length;
-    const reversed: string[] = new Array(len);
+    get ringComponents() {
+        return this.index.ringComponents;
+    }
 
-    for (let i = 0; i < len; i++) reversed[i] = elements[len - i - 1];
+    constructor(all: ReadonlyArray<UnitRing>, public unit: Unit.Atomic) {
+        this.all = all;
+    }
+}
 
-    const rotNormal = getMinimalRotation(elements);
-    const rotReversed = getMinimalRotation(reversed);
+namespace UnitRing {
+    export type Fingerprint = { readonly '@type': 'unit-ring-fingerprint' } & string
 
-    let isNormalSmaller = false;
+    export function fingerprint(unit: Unit.Atomic, ring: UnitRing): Fingerprint {
+        const { elements } = unit;
+        const { type_symbol } = unit.model.atomicHierarchy.atoms;
 
-    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;
-        }
+        const symbols: ElementSymbol[] = [];
+        for (let i = 0, _i = ring.length; i < _i; i++) symbols[symbols.length] = type_symbol.value(elements[ring[i]]);
+        return elementFingerprint(symbols);
     }
 
-    if (isNormalSmaller) return buildFinderprint(elements, rotNormal);
-    return buildFinderprint(reversed, rotReversed);
+    export function elementFingerprint(elements: ArrayLike<ElementSymbol>) {
+        return getFingerprint(elements as ArrayLike<String> as string[]) as Fingerprint;
+    }
 }
 
-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;
+namespace UnitRings {
+    /** Index into UnitRings.all */
+    export type Index = { readonly '@type': 'unit-ring-index' } & number
+    export type ComponentIndex = { readonly '@type': 'unit-ring-component-index' } & number
 
-    let u = '', v = '', k = 0;
+    export function create(unit: Unit.Atomic): UnitRings {
+        const rings = computeRings(unit);
+        return new UnitRings(rings, unit);
+    }
 
-    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];
+    /** Creates a mapping ResidueIndex -> list or rings that are on that residue and have one of the specified fingerprints. */
+    export function byFingerprintAndResidue(rings: UnitRings, fingerprints: ReadonlyArray<UnitRing.Fingerprint>) {
+        const map = new Map<ResidueIndex, Index[]>();
+        for (const fp of fingerprints) {
+            addSingleResidueRings(rings, fp, map);
         }
+        return map;
+    }
+}
 
-        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;
+function createByFingerprint(unit: Unit.Atomic, rings: ReadonlyArray<UnitRing>) {
+    const byFingerprint = new Map<UnitRing.Fingerprint, UnitRings.Index[]>();
+    let idx = 0 as UnitRings.Index;
+    for (const r of rings) {
+        const fp = UnitRing.fingerprint(unit, r);
+        if (byFingerprint.has(fp)) byFingerprint.get(fp)!.push(idx);
+        else byFingerprint.set(fp, [idx]);
+        idx++;
     }
+    return byFingerprint;
+}
 
-    return k;
+function ringResidueIdx(unit: Unit.Atomic, ring: ArrayLike<StructureElement.UnitIndex>): ResidueIndex {
+    const { elements } = unit;
+    const residueIndex = unit.model.atomicHierarchy.residueAtomSegments.index;
+    const idx = residueIndex[elements[ring[0]]];
+    for (let rI = 1, _rI = ring.length; rI < _rI; rI++) {
+        if (idx !== residueIndex[elements[ring[rI]]]) return -1 as ResidueIndex;
+    }
+    return idx;
 }
 
-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('-');
+function addSingleResidueRings(rings: UnitRings, fp: UnitRing.Fingerprint, map: Map<ResidueIndex, UnitRings.Index[]>) {
+    const byFp = rings.byFingerprint.get(fp);
+    if (!byFp) return;
+    for (const r of byFp) {
+        const idx = ringResidueIdx(rings.unit, rings.all[r]);
+        if (idx >= 0) {
+            if (map.has(idx)) map.get(idx)!.push(r);
+            else map.set(idx, [r]);
+        }
     }
-    ret.push(elements[(i + offset) % len]);
-    return ret.join('');
-}
+}
+
+
+export { UnitRing, UnitRings }

+ 123 - 6
src/mol-model/structure/structure/unit/rings/compute.ts

@@ -4,13 +4,15 @@
  * @author David Sehnal <david.sehnal@gmail.com>
  */
 
-import Unit from '../../unit';
-import { IntraUnitLinks } from '../links/data';
-import { Segmentation } from 'mol-data/int';
+import { Segmentation, SortedArray } from 'mol-data/int';
+import { IntAdjacencyGraph } from 'mol-math/graph';
 import { LinkType } from '../../../model/types';
 import { StructureElement } from '../../../structure';
+import Unit from '../../unit';
+import { IntraUnitLinks } from '../links/data';
+import { sortArray } from 'mol-data/util';
 
-export default function computeRings(unit: Unit.Atomic) {
+export function computeRings(unit: Unit.Atomic) {
     const size = largestResidue(unit);
     const state = State(unit, size);
 
@@ -41,7 +43,7 @@ interface State {
 
     currentColor: number,
 
-    rings: StructureElement.UnitIndex[][],
+    rings: SortedArray<StructureElement.UnitIndex>[],
     bonds: IntraUnitLinks,
     unit: Unit.Atomic
 }
@@ -145,7 +147,8 @@ function addRing(state: State, a: number, b: number) {
     for (let t = 0; t < leftOffset; t++) ring[ringOffset++] = state.startVertex + left[t];
     for (let t = rightOffset - 1; t >= 0; t--) ring[ringOffset++] = state.startVertex + right[t];
 
-    state.rings.push(ring as any as StructureElement.UnitIndex[]);
+    sortArray(ring);
+    state.rings.push(SortedArray.ofSortedArray(ring));
 }
 
 function findRings(state: State, from: number) {
@@ -176,4 +179,118 @@ 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: ArrayLike<SortedArray<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 graph = new IntAdjacencyGraph.UniqueEdgeBuilder(rings.length);
+    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;
+                graph.addEdge(rI, rJ);
+            }
+        }
+    }
+
+    const components = IntAdjacencyGraph.connectedComponents(graph.getGraph());
+
+    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 };
 }