Added UnitRings

src/apps/structure-info/model.ts

@@ -18,6 +18,7 @@ import { mmCIF_Database } from 'mol-io/reader/cif/schema/mmcif';
 import { openCif, downloadCif } from './helpers';
 import { BitFlags } from 'mol-util';
 import { SecondaryStructureType } from 'mol-model/structure/model/types';
+import { UnitRings } from 'mol-model/structure/structure/unit/rings';
 
 
 async function downloadFromPdb(pdb: string) {
@@ -100,6 +101,21 @@ export function printSequence(model: Model) {
     console.log();
 }
 
+export function printRings(structure: Structure) {
+    console.log('\nRings\n=============');
+    for (const unit of structure.units) {
+        if (!Unit.isAtomic(unit)) continue;
+        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]);
+        }
+        if (all.length > 5) fps.push('...')
+        console.log(`Unit ${unit.id}, ${all.length} ring(s), ${byFingerprint.size} different fingerprint(s).\n  ${fps.join(', ')}`);
+    }
+    console.log();
+}
+
 export function printUnits(structure: Structure) {
     console.log('\nUnits\n=============');
     const l = Element.Location();
@@ -143,6 +159,7 @@ async function run(mmcif: mmCIF_Database) {
     printSequence(models[0]);
     //printIHMModels(models[0]);
     printUnits(structure);
+    printRings(structure);
     //printBonds(structure);
     //printSecStructure(models[0]);
 }

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

@@ -12,6 +12,7 @@ import { idFactory } from 'mol-util/id-factory';
 import { IntraUnitBonds, computeIntraUnitBonds } from './unit/bonds'
 import { CoarseElements, CoarseSphereConformation, CoarseGaussianConformation } from '../model/properties/coarse';
 import { ValueRef } from 'mol-util';
+import { UnitRings } from './unit/rings';
 
 // A building block of a structure that corresponds to an atomic or a coarse grained representation
 // 'conveniently grouped together'.
@@ -97,6 +98,12 @@ namespace Unit {
             return this.props.bonds.ref;
         }
 
+        get rings() {
+            if (this.props.rings.ref) return this.props.rings.ref;
+            this.props.rings.ref = UnitRings.create(this);
+            return this.props.rings.ref;
+        }
+
         constructor(id: number, invariantId: number, model: Model, elements: SortedArray, conformation: SymmetryOperator.ArrayMapping, props: AtomicProperties) {
             this.id = id;
             this.invariantId = invariantId;
@@ -113,10 +120,11 @@ namespace Unit {
     interface AtomicProperties {
         lookup3d: ValueRef<Lookup3D | undefined>,
         bonds: ValueRef<IntraUnitBonds | undefined>,
+        rings: ValueRef<UnitRings | undefined>
     }
 
     function AtomicProperties() {
-        return { lookup3d: ValueRef.create(void 0), bonds: ValueRef.create(void 0) };
+        return { lookup3d: ValueRef.create(void 0), bonds: ValueRef.create(void 0), rings: ValueRef.create(void 0) };
     }
 
     class Coarse<K extends Kind.Gaussians | Kind.Spheres, C extends CoarseSphereConformation | CoarseGaussianConformation> implements Base {

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

@@ -0,0 +1,107 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import computeRings from './rings/compute'
+import Unit from '../unit';
+
+interface UnitRings {
+    /** Each ring is specified as an array of indices in Unit.elements. */
+    readonly all: ReadonlyArray<ReadonlyArray<number>>,
+    readonly byFingerprint: Map<string, ReadonlyArray<number>>
+}
+
+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);
+    }
+
+    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++;
+        }
+
+        return { all: rings, byFingerprint };
+    }
+}
+
+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('');
+}

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

@@ -0,0 +1,177 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import Unit from '../../unit';
+import { IntraUnitBonds } from '../bonds/intra-data';
+import { Segmentation } from 'mol-data/int';
+
+export default function computeRings(unit: Unit.Atomic) {
+    const size = largestResidue(unit);
+    const state = State(unit, size);
+
+    const residuesIt = Segmentation.transientSegments(unit.model.atomicHierarchy.residueSegments, unit.elements);
+    while (residuesIt.hasNext) {
+        const seg = residuesIt.move();
+        processResidue(state, seg.start, seg.end);
+    }
+
+    return state.rings;
+}
+
+const enum Constants {
+    MaxDepth = 4
+}
+
+interface State {
+    startVertex: number,
+    endVertex: number,
+    count: number,
+    visited: Int32Array,
+    queue: Int32Array,
+    color: Int32Array,
+    pred: Int32Array,
+
+    left: Int32Array,
+    right: Int32Array,
+
+    currentColor: number,
+
+    rings: number[][],
+    bonds: IntraUnitBonds,
+    unit: Unit.Atomic
+}
+
+function State(unit: Unit.Atomic, capacity: number): State {
+    return {
+        startVertex: 0,
+        endVertex: 0,
+        count: 0,
+        visited: new Int32Array(capacity),
+        queue: new Int32Array(capacity),
+        pred: new Int32Array(capacity),
+        left: new Int32Array(Constants.MaxDepth),
+        right: new Int32Array(Constants.MaxDepth),
+        color: new Int32Array(capacity),
+        currentColor: 0,
+        rings: [],
+        unit,
+        bonds: unit.bonds
+    };
+}
+
+function resetState(state: State) {
+    state.count = state.endVertex - state.startVertex;
+    const { visited, pred, color } = state;
+    for (let i = 0; i < state.count; i++) {
+        visited[i] = -1;
+        pred[i] = -1;
+        color[i] = 0;
+    }
+    state.currentColor = 0;
+}
+
+function largestResidue(unit: Unit.Atomic) {
+    const residuesIt = Segmentation.transientSegments(unit.model.atomicHierarchy.residueSegments, unit.elements);
+    let size = 0;
+    while (residuesIt.hasNext) {
+        const seg = residuesIt.move();
+        size = Math.max(size, seg.end - seg.start);
+    }
+    return size;
+}
+
+function processResidue(state: State, start: number, end: number) {
+    const { visited } = state;
+    state.startVertex = start;
+    state.endVertex = end;
+
+    // no two atom rings
+    if (state.endVertex - state.startVertex < 3) return;
+
+    resetState(state);
+
+    for (let i = 0; i < state.count; i++) {
+        if (visited[i] >= 0) continue;
+        findRings(state, i);
+    }
+}
+
+function addRing(state: State, a: number, b: number) {
+    // only "monotonous" rings
+    if (b < a) return;
+
+    const { pred, color, left, right } = state;
+    const nc = ++state.currentColor;
+
+    let current = a;
+
+    for (let t = 0; t < Constants.MaxDepth; t++) {
+        color[current] = nc;
+        current = pred[current];
+        if (current < 0) break;
+    }
+
+    let leftOffset = 0, rightOffset = 0;
+
+    let found = false, target = 0;
+    current = b;
+    for (let t = 0; t < Constants.MaxDepth; t++) {
+        if (color[current] === nc) {
+            target = current;
+            found = true;
+            break;
+        }
+        right[rightOffset++] = current;
+        current = pred[current];
+        if (current < 0) break;
+    }
+    if (!found) return;
+
+    current = a;
+    for (let t = 0; t < Constants.MaxDepth; t++) {
+        left[leftOffset++] = current;
+        if (target === current) break;
+        current = pred[current];
+        if (current < 0) break;
+    }
+
+    const ring = new Int32Array(leftOffset + rightOffset);
+    let ringOffset = 0;
+    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 number[]);
+}
+
+function findRings(state: State, from: number) {
+    const { bonds, startVertex, endVertex, visited, queue, pred } = state;
+    const { b: neighbor, flags: bondFlags, offset } = bonds;
+    visited[from] = 1;
+    queue[0] = from;
+    let head = 0, size = 1;
+
+    while (head < size) {
+        const top = queue[head++];
+        const a = startVertex + top;
+        const start = offset[a], end = offset[a + 1];
+
+        for (let i = start; i < end; i++) {
+            const b = neighbor[i];
+            if (b < startVertex || b >= endVertex || !IntraUnitBonds.isCovalent(bondFlags[i])) continue;
+
+            const other = b - startVertex;
+
+            if (visited[other] > 0) {
+                if (pred[other] !== top && pred[top] !== other) addRing(state, top, other);
+                continue;
+            }
+
+            visited[other] = 1;
+            queue[size++] = other;
+            pred[other] = top;
+        }
+    }
+}