Browse Source

added InterUnitGraph.Builder

Alexander Rose 5 years ago
parent
commit
116d3ec319

+ 51 - 1
src/mol-math/graph/inter-unit-graph.ts

@@ -1,5 +1,7 @@
+import { UniqueArray } from '../../mol-data/generic'
+
 /**
- * Copyright (c) 2017-2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
@@ -123,6 +125,54 @@ namespace InterUnitGraph {
     export function getVertexKey<Unit extends UnitBase, VertexIndex extends number>(index: VertexIndex, unit: Unit) {
         return `${index}|${unit.id}`
     }
+
+    //
+
+    function addMapEntry<A, B>(map: Map<A, B[]>, a: A, b: B) {
+        if (map.has(a)) map.get(a)!.push(b);
+        else map.set(a, [b]);
+    }
+
+
+    export class Builder<Unit extends InterUnitGraph.UnitBase, VertexIndex extends number, EdgeProps extends InterUnitGraph.EdgePropsBase = {}> {
+        private uA: Unit
+        private uB: Unit
+        private mapAB: Map<number, EdgeInfo<VertexIndex, EdgeProps>[]>
+        private mapBA: Map<number, EdgeInfo<VertexIndex, EdgeProps>[]>
+        private linkedA: UniqueArray<VertexIndex, VertexIndex>
+        private linkedB: UniqueArray<VertexIndex, VertexIndex>
+        private linkCount: number
+
+        private map = new Map<number, UnitPairEdges<Unit, VertexIndex, EdgeProps>[]>();
+
+        startUnitPair(unitA: Unit, unitB: Unit) {
+            this.uA = unitA
+            this.uB = unitB
+            this.mapAB = new Map()
+            this.mapBA = new Map()
+            this.linkedA = UniqueArray.create()
+            this.linkedB = UniqueArray.create()
+            this.linkCount = 0
+        }
+
+        finishUnitPair() {
+            if (this.linkCount === 0) return
+            addMapEntry(this.map, this.uA.id, new UnitPairEdges(this.uA, this.uB, this.linkCount, this.linkedA.array, this.mapAB))
+            addMapEntry(this.map, this.uB.id, new UnitPairEdges(this.uB, this.uA, this.linkCount, this.linkedB.array, this.mapBA))
+        }
+
+        add(indexA: VertexIndex, indexB: VertexIndex, props: EdgeProps) {
+            addMapEntry(this.mapAB, indexA, { indexB, props })
+            addMapEntry(this.mapBA, indexB, { indexB: indexA, props })
+            UniqueArray.add(this.linkedA, indexA, indexA)
+            UniqueArray.add(this.linkedB, indexB, indexB)
+            this.linkCount += 1
+        }
+
+        getMap(): Map<number, InterUnitGraph.UnitPairEdges<Unit, VertexIndex, EdgeProps>[]> {
+            return this.map;
+        }
+    }
 }
 
 const emptyArray: any[] = [];

+ 11 - 12
src/mol-model-props/computed/interactions/common.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2019-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
@@ -12,15 +12,16 @@ import { Features } from './features'
 import { StructureElement } from '../../../mol-model/structure/structure'
 import { IntMap } from '../../../mol-data/int'
 
-type IntraProps = {
-    readonly type: ArrayLike<InteractionType>
-    readonly flag: AssignableArrayLike<InteractionFlag>
-}
 export { InteractionsIntraContacts }
-interface InteractionsIntraContacts extends IntAdjacencyGraph<Features.FeatureIndex, IntraProps> {
+interface InteractionsIntraContacts extends IntAdjacencyGraph<Features.FeatureIndex, InteractionsIntraContacts.Props> {
     readonly elementsIndex: InteractionsIntraContacts.ElementsIndex
 }
 namespace InteractionsIntraContacts {
+    export type Props = {
+        readonly type: ArrayLike<InteractionType>
+        readonly flag: AssignableArrayLike<InteractionFlag>
+    }
+
     /** maps unit elements to contacts, range for unit element i is offsets[i] to offsets[i + 1] */
     export type ElementsIndex = {
         /** intra contact indices */
@@ -32,7 +33,7 @@ namespace InteractionsIntraContacts {
     /**
      * Note: assumes that feature members of a contact are non-overlapping
      */
-    export function createElementsIndex(contacts: IntAdjacencyGraph<Features.FeatureIndex, IntraProps>, features: Features, elementsCount: number) {
+    export function createElementsIndex(contacts: IntAdjacencyGraph<Features.FeatureIndex, Props>, features: Features, elementsCount: number) {
         const offsets = new Int32Array(elementsCount + 1)
         const bucketFill = new Int32Array(elementsCount)
         const bucketSizes = new Int32Array(elementsCount)
@@ -83,8 +84,7 @@ namespace InteractionsIntraContacts {
 }
 
 export { InteractionsInterContacts }
-type InterProps = { type: InteractionType, flag: InteractionFlag }
-class InteractionsInterContacts extends InterUnitGraph<Unit, Features.FeatureIndex, InterProps> {
+class InteractionsInterContacts extends InterUnitGraph<Unit, Features.FeatureIndex, InteractionsInterContacts.Props> {
     private readonly elementKeyIndex: Map<string, number[]>
 
     getContactIndicesForElement(index: StructureElement.UnitIndex, unit: Unit): ReadonlyArray<number> {
@@ -95,7 +95,7 @@ class InteractionsInterContacts extends InterUnitGraph<Unit, Features.FeatureInd
         return `${index}|${unit.id}`
     }
 
-    constructor(map: Map<number, InterUnitGraph.UnitPairEdges<Unit, Features.FeatureIndex, InterProps>[]>, unitsFeatures: IntMap<Features>) {
+    constructor(map: Map<number, InterUnitGraph.UnitPairEdges<Unit, Features.FeatureIndex, InteractionsInterContacts.Props>[]>, unitsFeatures: IntMap<Features>) {
         super(map)
 
         let count = 0
@@ -135,8 +135,7 @@ class InteractionsInterContacts extends InterUnitGraph<Unit, Features.FeatureInd
     }
 }
 namespace InteractionsInterContacts {
-    export class Pair extends InterUnitGraph.UnitPairEdges<Unit, Features.FeatureIndex, InterProps> {}
-    export type Info = InterUnitGraph.EdgeInfo<Features.FeatureIndex, InterProps>
+    export type Props = { type: InteractionType, flag: InteractionFlag }
 }
 
 export const enum InteractionFlag {

+ 7 - 32
src/mol-model-props/computed/interactions/contacts-builder.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2019-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
@@ -8,9 +8,9 @@ import { Features } from './features';
 import { IntAdjacencyGraph } from '../../../mol-math/graph';
 import { InteractionType, InteractionsIntraContacts, InteractionsInterContacts, InteractionFlag } from './common';
 import { Unit } from '../../../mol-model/structure/structure';
-import { UniqueArray } from '../../../mol-data/generic';
 import { NumberArray } from '../../../mol-util/type-helpers';
 import { IntMap } from '../../../mol-data/int';
+import { InterUnitGraph } from '../../../mol-math/graph/inter-unit-graph';
 
 export { IntraContactsBuilder }
 
@@ -63,46 +63,21 @@ interface InterContactsBuilder {
 }
 
 namespace InterContactsBuilder {
-    function addMapEntry<A, B>(map: Map<A, B[]>, a: A, b: B) {
-        if (map.has(a)) map.get(a)!.push(b);
-        else map.set(a, [b]);
-    }
-
     export function create(): InterContactsBuilder {
-        let uA: Unit
-        let uB: Unit
-        let mapAB: Map<number, InteractionsInterContacts.Info[]>
-        let mapBA: Map<number, InteractionsInterContacts.Info[]>
-        let linkedA: UniqueArray<Features.FeatureIndex, Features.FeatureIndex>
-        let linkedB: UniqueArray<Features.FeatureIndex, Features.FeatureIndex>
-        let linkCount: number
-
-        const map = new Map<number, InteractionsInterContacts.Pair[]>();
+        const builder = new InterUnitGraph.Builder<Unit, Features.FeatureIndex, InteractionsInterContacts.Props>()
 
         return {
             startUnitPair(unitA: Unit, unitB: Unit) {
-                uA = unitA
-                uB = unitB
-                mapAB = new Map()
-                mapBA = new Map()
-                linkedA = UniqueArray.create()
-                linkedB = UniqueArray.create()
-                linkCount = 0
+                builder.startUnitPair(unitA, unitB)
             },
             finishUnitPair() {
-                if (linkCount === 0) return
-                addMapEntry(map, uA.id, new InteractionsInterContacts.Pair(uA, uB, linkCount, linkedA.array, mapAB))
-                addMapEntry(map, uB.id, new InteractionsInterContacts.Pair(uB, uA, linkCount, linkedB.array, mapBA))
+                builder.finishUnitPair()
             },
             add(indexA: Features.FeatureIndex, indexB: Features.FeatureIndex, type: InteractionType) {
-                addMapEntry(mapAB, indexA, { indexB, props: { type, flag: InteractionFlag.None } })
-                addMapEntry(mapBA, indexB, { indexB: indexA, props: { type, flag: InteractionFlag.None } })
-                UniqueArray.add(linkedA, indexA, indexA)
-                UniqueArray.add(linkedB, indexB, indexB)
-                linkCount += 1
+                builder.add(indexA, indexB, { type, flag: InteractionFlag.None })
             },
             getContacts(unitsFeatures: IntMap<Features>) {
-                return new InteractionsInterContacts(map, unitsFeatures);
+                return new InteractionsInterContacts(builder.getMap(), unitsFeatures);
             }
         }
     }

+ 1 - 1
src/mol-model/structure/structure/unit/bonds/data.ts

@@ -37,4 +37,4 @@ namespace InterUnitBonds {
     export type BondInfo = InterUnitGraph.EdgeInfo<StructureElement.UnitIndex, InterUnitEdgeProps>
 }
 
-export { IntraUnitBonds, InterUnitBonds }
+export { IntraUnitBonds, InterUnitBonds, InterUnitEdgeProps }

+ 25 - 57
src/mol-model/structure/structure/unit/bonds/inter-compute.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2017-2019 Mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2020 Mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
@@ -9,8 +9,7 @@ import { BondType } from '../../../model/types';
 import Structure from '../../structure';
 import Unit from '../../unit';
 import { getElementIdx, getElementPairThreshold, getElementThreshold, isHydrogen, BondComputationProps, MetalsSet, DefaultBondComputationProps } from './common';
-import { InterUnitBonds } from './data';
-import { UniqueArray } from '../../../../../mol-data/generic';
+import { InterUnitBonds, InterUnitEdgeProps } from './data';
 import { SortedArray } from '../../../../../mol-data/int';
 import { Vec3, Mat4 } from '../../../../../mol-math/linear-algebra';
 import StructureElement from '../../element';
@@ -18,29 +17,10 @@ import { StructConn } from '../../../../../mol-model-formats/structure/mmcif/bon
 import { ElementIndex } from '../../../model/indexing';
 import { getInterBondOrderFromTable } from '../../../model/properties/atomic/bonds';
 import { IndexPairBonds } from '../../../../../mol-model-formats/structure/mmcif/bonds/index-pair';
+import { InterUnitGraph } from '../../../../../mol-math/graph/inter-unit-graph';
 
 const MAX_RADIUS = 4;
 
-function addMapEntry<A, B>(map: Map<A, B[]>, a: A, b: B) {
-    if (map.has(a)) map.get(a)!.push(b);
-    else map.set(a, [b]);
-}
-
-interface PairState {
-    mapAB: Map<number, InterUnitBonds.BondInfo[]>,
-    mapBA: Map<number, InterUnitBonds.BondInfo[]>,
-    bondedA: UniqueArray<StructureElement.UnitIndex, StructureElement.UnitIndex>,
-    bondedB: UniqueArray<StructureElement.UnitIndex, StructureElement.UnitIndex>
-}
-
-function addBond(indexA: StructureElement.UnitIndex, indexB: StructureElement.UnitIndex, order: number, flag: BondType.Flag, state: PairState) {
-    addMapEntry(state.mapAB, indexA, { indexB, props: { order, flag } });
-    addMapEntry(state.mapBA, indexB, { indexB: indexA, props: { order, flag } });
-
-    UniqueArray.add(state.bondedA, indexA, indexA);
-    UniqueArray.add(state.bondedB, indexB, indexB);
-}
-
 const tmpDistVecA = Vec3()
 const tmpDistVecB = Vec3()
 function getDistance(unitA: Unit.Atomic, indexA: ElementIndex, unitB: Unit.Atomic, indexB: ElementIndex) {
@@ -51,17 +31,13 @@ function getDistance(unitA: Unit.Atomic, indexA: ElementIndex, unitB: Unit.Atomi
 
 const _imageTransform = Mat4.zero();
 
-function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComputationProps, map: Map<number, InterUnitBonds.UnitPairBonds[]>) {
-    const state: PairState = { mapAB: new Map(), mapBA: new Map(), bondedA: UniqueArray.create(), bondedB: UniqueArray.create() };
-    let bondCount = 0;
+function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComputationProps, builder: InterUnitGraph.Builder<Unit.Atomic, StructureElement.UnitIndex, InterUnitEdgeProps>) {
 
     const { elements: atomsA, residueIndex: residueIndexA } = unitA;
     const { x: xA, y: yA, z: zA } = unitA.model.atomicConformation;
     const { elements: atomsB, residueIndex: residueIndexB } = unitB;
     const atomCount = unitA.elements.length;
 
-    // const { type_symbol: type_symbolA, label_alt_id: label_alt_idA } = unitA.model.atomicHierarchy.atoms;
-    // const { type_symbol: type_symbolB, label_alt_id: label_alt_idB } = unitB.model.atomicHierarchy.atoms;
     const { type_symbol: type_symbolA, label_alt_id: label_alt_idA, label_atom_id: label_atom_idA } = unitA.model.atomicHierarchy.atoms;
     const { type_symbol: type_symbolB, label_alt_id: label_alt_idB, label_atom_id: label_atom_idB } = unitB.model.atomicHierarchy.atoms;
     const { label_comp_id: label_comp_idA, auth_seq_id: auth_seq_idA } = unitA.model.atomicHierarchy.residues;
@@ -82,6 +58,8 @@ function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComput
     const { center: bCenter, radius: bRadius } = lookup3d.boundary.sphere;
     const testDistanceSq = (bRadius + MAX_RADIUS) * (bRadius + MAX_RADIUS);
 
+    builder.startUnitPair(unitA, unitB)
+
     for (let _aI = 0 as StructureElement.UnitIndex; _aI < atomCount; _aI++) {
         const aI = atomsA[_aI];
         Vec3.set(imageA, xA[aI], yA[aI], zA[aI]);
@@ -92,8 +70,7 @@ function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComput
             for (let i = indexPairs.offset[aI], il = indexPairs.offset[aI + 1]; i < il; ++i) {
                 const _bI = SortedArray.indexOf(unitA.elements, indexPairs.b[i]) as StructureElement.UnitIndex;
                 if (_bI < 0) continue;
-                addBond(_aI, _bI, indexPairs.edgeProps.order[i], BondType.Flag.Covalent, state);
-                bondCount++;
+                builder.add(_aI, _bI, { order: indexPairs.edgeProps.order[i], flag: BondType.Flag.Covalent });
             }
             continue // assume `indexPairs` supplies all bonds
         }
@@ -107,8 +84,7 @@ function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComput
                     if (_bI < 0) continue;
                     // check if the bond is within MAX_RADIUS for this pair of units
                     if (getDistance(unitA, aI, unitB, p.atomIndex) > MAX_RADIUS) continue;
-                    addBond(_aI, _bI, se.order, se.flags, state);
-                    bondCount++;
+                    builder.add(_aI, _bI, { order: se.order, flag: se.flags });
                     added = true;
                 }
             }
@@ -136,8 +112,10 @@ function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComput
             const altB = label_alt_idB.value(bI);
             if (altA && altB && altA !== altB) continue;
 
-            // Do not include bonds between images of the same residue.
+            // Do not include bonds between images of the same residue with partial occupancy.
             // TODO: is this condition good enough?
+            // - It works for cases like 3WQJ (label_asym_id: I) which have partial occupancy.
+            // - Does NOT work for cases like 1RB8 (DC 7) with full occupancy.
             if (occupancyB.value(bI) < 1 && occA < 1)  {
                 if (auth_seq_idA.value(aI) === auth_seq_idB.value(bI)) {
                     continue;
@@ -155,13 +133,11 @@ function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComput
 
             if (isHa || isHb) {
                 if (dist < props.maxCovalentHydrogenBondingLength) {
-                    addBond(
-                        _aI, _bI,
-                        1, // covalent bonds involving a hydrogen are always of order 1
-                        BondType.Flag.Covalent | BondType.Flag.Computed,
-                        state
-                    );
-                    bondCount++;
+                    // covalent bonds involving a hydrogen are always of order 1
+                    builder.add(_aI, _bI, {
+                        order: 1,
+                        flag: BondType.Flag.Covalent | BondType.Flag.Computed
+                    });
                 }
                 continue;
             }
@@ -174,23 +150,15 @@ function findPairBonds(unitA: Unit.Atomic, unitB: Unit.Atomic, props: BondComput
             if (dist <= pairingThreshold) {
                 const atomIdB = label_atom_idB.value(bI);
                 const compIdB = label_comp_idB.value(residueIndexB[bI]);
-                addBond(
-                    _aI, _bI,
-                    getInterBondOrderFromTable(compIdA, compIdB, atomIdA, atomIdB),
-                    (isMetal ? BondType.Flag.MetallicCoordination : BondType.Flag.Covalent) | BondType.Flag.Computed,
-                    state
-                );
-                bondCount++;
+                builder.add(_aI, _bI, {
+                    order: getInterBondOrderFromTable(compIdA, compIdB, atomIdA, atomIdB),
+                    flag: (isMetal ? BondType.Flag.MetallicCoordination : BondType.Flag.Covalent) | BondType.Flag.Computed
+                });
             }
         }
     }
 
-    if (bondCount) {
-        addMapEntry(map, unitA.id, new InterUnitBonds.UnitPairBonds(unitA, unitB, bondCount, state.bondedA.array, state.mapAB));
-        addMapEntry(map, unitB.id, new InterUnitBonds.UnitPairBonds(unitB, unitA, bondCount, state.bondedB.array, state.mapBA));
-    }
-
-    return bondCount;
+    builder.finishUnitPair()
 }
 
 export interface InterBondComputationProps extends BondComputationProps {
@@ -198,24 +166,24 @@ export interface InterBondComputationProps extends BondComputationProps {
 }
 
 function findBonds(structure: Structure, props: InterBondComputationProps) {
-    const map = new Map<number, InterUnitBonds.UnitPairBonds[]>();
+    const builder = new InterUnitGraph.Builder<Unit.Atomic, StructureElement.UnitIndex, InterUnitEdgeProps>()
 
     if (props.noCompute) {
         // TODO add function that only adds bonds defined in structConn and avoids using
         //      structure.lookup and unit.lookup (expensive for large structure and not
         //      needed for archival files or files with an MD topology)
-        return new InterUnitBonds(map);
+        return new InterUnitBonds(builder.getMap());
     }
 
     Structure.eachUnitPair(structure, (unitA: Unit, unitB: Unit) => {
-        findPairBonds(unitA as Unit.Atomic, unitB as Unit.Atomic, props, map)
+        findPairBonds(unitA as Unit.Atomic, unitB as Unit.Atomic, props, builder)
     }, {
         maxRadius: MAX_RADIUS,
         validUnit: (unit: Unit) => Unit.isAtomic(unit),
         validUnitPair: (unitA: Unit, unitB: Unit) => props.validUnitPair(structure, unitA, unitB)
     })
 
-    return new InterUnitBonds(map);
+    return new InterUnitBonds(builder.getMap());
 }
 
 function computeInterUnitBonds(structure: Structure, props?: Partial<InterBondComputationProps>): InterUnitBonds {