diff --git a/package.json b/package.json index 72dd74f55720e4c2e513614661ceaa78b274da4e..ded444b432fbf6c7a1438edce5b5c016843689db 100644 --- a/package.json +++ b/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", diff --git a/src/apps/structure-info/model.ts b/src/apps/structure-info/model.ts index c04627418a078014a7f89f061b8f0d5bbde55fef..fbc466845b97f16694d61f36f8afffd2462e7219 100644 --- a/src/apps/structure-info/model.ts +++ b/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(', ')}`); diff --git a/src/mol-data/util/hash-functions.ts b/src/mol-data/util/hash-functions.ts index 8584baa5946e9a06a4ef1ca8e07b8e43607028b9..d822d8d9c4c1f60d77106d1d2b00526d146fc75b 100644 --- a/src/mol-data/util/hash-functions.ts +++ b/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); } \ No newline at end of file diff --git a/src/mol-math/graph/int-adjacency-graph.ts b/src/mol-math/graph/int-adjacency-graph.ts index 4f74689e4b58fff99786b3d0cec063a87c308fc8..73ee9a0beafe0c853afc61cfaede7707af811ee4 100644 --- a/src/mol-math/graph/int-adjacency-graph.ts +++ b/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; } \ No newline at end of file diff --git a/src/mol-model/structure/structure.ts b/src/mol-model/structure/structure.ts index 3a924335ea86458fea0e1d2cb89d2be92310bdee..a6a25607900a4ac9912f878f6ac5ddc9abf9541b 100644 --- a/src/mol-model/structure/structure.ts +++ b/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 } \ No newline at end of file +export { StructureElement, Link, Structure, Unit, StructureSymmetry, StructureProperties } +export * from './structure/unit/rings' \ No newline at end of file diff --git a/src/mol-model/structure/structure/carbohydrates/compute.ts b/src/mol-model/structure/structure/carbohydrates/compute.ts index 308780d5ae6f341f9a166cd2cb1175ac957a2c4c..e0bcf3e61c7f8c29dd4b349c9750c58ddd323993 100644 --- a/src/mol-model/structure/structure/carbohydrates/compute.ts +++ b/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) diff --git a/src/mol-model/structure/structure/unit.ts b/src/mol-model/structure/structure/unit.ts index 48f429a29dc65e7f1d883669d7ea6d112635a57a..8ec99349308a246bf74320c2b517233d85999e8b 100644 --- a/src/mol-model/structure/structure/unit.ts +++ b/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; diff --git a/src/mol-model/structure/structure/unit/links/data.ts b/src/mol-model/structure/structure/unit/links/data.ts index ed046a026da9011b17f43602898b8c2f3e1fdf69..74ff5bd2d5aade0d8f93d0fed3bedf026ea75d63 100644 --- a/src/mol-model/structure/structure/unit/links/data.ts +++ b/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 } diff --git a/src/mol-model/structure/structure/unit/links/inter-compute.ts b/src/mol-model/structure/structure/unit/links/inter-compute.ts index 03184c454ff8674d8471d65f76c442721c166f30..cd8e37397544128778502efe23fa2b664c1009a1 100644 --- a/src/mol-model/structure/structure/unit/links/inter-compute.ts +++ b/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) { diff --git a/src/mol-model/structure/structure/unit/rings.ts b/src/mol-model/structure/structure/unit/rings.ts index 4b92e5b914f77b682a6c3b5365ae6de87c8affba..64f9f857cb0f736216ea46ed887bbd84ba41cbc9 100644 --- a/src/mol-model/structure/structure/unit/rings.ts +++ b/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(''); -} \ No newline at end of file +} + + +export { UnitRing, UnitRings } \ No newline at end of file diff --git a/src/mol-model/structure/structure/unit/rings/compute.ts b/src/mol-model/structure/structure/unit/rings/compute.ts index 0225846f647a147a8911211fec2d0efef065dca6..4517fb879f826cd4934682f6ab247d5075084f11 100644 --- a/src/mol-model/structure/structure/unit/rings/compute.ts +++ b/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 }; } \ No newline at end of file