wip, adding bond and sec struct related code from molql

src/apps/combine-mmcif/index.ts

@@ -208,8 +208,7 @@ export function getBirdBonds(mmcif: mmCIF_Database) {
 }
 
 export function getCcdBonds(mmcif: mmCIF_Database) {
-    const bonds: PartialStructConnRow[] = []
-
+    // const bonds: PartialStructConnRow[] = []
 }
 
 async function run(pdb: string, out?: string) {

src/apps/structure-info/index.ts

@@ -0,0 +1,63 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import * as argparse from 'argparse'
+import fetch from 'node-fetch'
+require('util.promisify').shim();
+
+import CIF from 'mol-io/reader/cif'
+import Computation from 'mol-util/computation'
+import { Model, Structure } from 'mol-model/structure'
+
+function showProgress(tag: string, p: Computation.Progress) {
+    console.log(`[${tag}] ${p.message} ${p.isIndeterminate ? '' : (p.current / p.max * 100).toFixed(2) + '% '}(${p.elapsedMs | 0}ms)`)
+}
+
+async function parseCif(data: string|Uint8Array) {
+    const comp = CIF.parse(data)
+    const ctx = Computation.observable({
+        updateRateMs: 250,
+        observer: p => showProgress(`cif parser ${typeof data === 'string' ? 'string' : 'binary'}`, p)
+    });
+    console.time('parse cif')
+    const parsed = await comp(ctx);
+    console.timeEnd('parse cif')
+    if (parsed.isError) throw parsed;
+    return parsed
+}
+
+export async function getPdb(pdb: string) {
+    console.log(`downloading ${pdb}...`)
+    const data = await fetch(`https://files.rcsb.org/download/${pdb}.cif`)
+    console.log(`done downloading ${pdb}`)
+
+    const parsed = await parseCif(await data.text())
+    return CIF.schema.mmCIF(parsed.result.blocks[0])
+}
+
+async function run(pdb: string, out?: string) {
+    const mmcif = await getPdb(pdb)
+
+    const models = Model.create({ kind: 'mmCIF', data: mmcif });
+    const structure = Structure.ofModel(models[0])
+
+    console.log(structure)
+    console.log(Model.bonds(models[0]))
+}
+
+const parser = new argparse.ArgumentParser({
+  addHelp: true,
+  description: 'Print info about a structure'
+});
+parser.addArgument([ '--pdb', '-p' ], {
+    help: 'Pdb entry id'
+});
+interface Args {
+    pdb: string
+}
+const args: Args = parser.parseArgs();
+
+run(args.pdb)

src/mol-math/geometry/grid-lookup.ts

@@ -1,2 +1,291 @@
+/*
+ * Copyright (c) 2017 MolQL contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
 // TODO: 3d grid lookup (use single cell for small sets), make bounding sphere part of the structure
-// TODO: assign radius to points?
+// TODO: assign radius to points?
+
+/**
+ * A "masked" 3D spatial lookup structure.
+ */
+
+import Mask from 'mol-util/mask'
+
+export type FindFunc = (x: number, y: number, z: number, radius: number) => Result
+export type CheckFunc = (x: number, y: number, z: number, radius: number) => boolean
+
+export interface Result {
+    readonly count: number,
+    readonly indices: number[],
+    readonly squaredDistances: number[]
+}
+
+export interface Positions { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number> }
+
+interface GridLookup { find: (mask: Mask) => FindFunc, check: (mask: Mask) => CheckFunc }
+
+function GridLookup(positions: Positions): GridLookup {
+    const tree = build(createInputData(positions));
+    return {
+        find(mask) {
+            const ctx = QueryContext.create(tree, mask, false);
+            return function (x: number, y: number, z: number, radius: number) {
+                QueryContext.update(ctx, x, y, z, radius);
+                nearest(ctx);
+                return ctx.buffer;
+            }
+        },
+        check(mask) {
+            const ctx = QueryContext.create(tree, mask, true);
+            return function (x: number, y: number, z: number, radius: number) {
+                QueryContext.update(ctx, x, y, z, radius);
+                return nearest(ctx);
+            }
+        }
+    }
+}
+
+interface InputData {
+    bounds: Box3D,
+    count: number,
+    positions: Positions
+}
+
+interface Box3D {
+    min: number[],
+    max: number[]
+}
+
+namespace Box3D {
+    export function createInfinite(): Box3D {
+        return {
+            min: [Number.MAX_VALUE, Number.MAX_VALUE, Number.MAX_VALUE],
+            max: [-Number.MAX_VALUE, -Number.MAX_VALUE, -Number.MAX_VALUE]
+        }
+    }
+}
+
+/**
+* Query context. Handles the actual querying.
+*/
+interface QueryContext<T> {
+    structure: T,
+    pivot: number[],
+    radius: number,
+    radiusSq: number,
+    buffer: QueryContext.Buffer,
+    mask: Mask,
+    isCheck: boolean
+}
+
+namespace QueryContext {
+    export interface Buffer extends Result {
+        count: number;
+        indices: any[];
+        squaredDistances: number[];
+    }
+
+    export function add<T>(ctx: QueryContext<T>, distSq: number, index: number) {
+        const buffer = ctx.buffer;
+        buffer.squaredDistances[buffer.count] = distSq;
+        buffer.indices[buffer.count++] = index;
+    }
+
+    function resetBuffer(buffer: Buffer) { buffer.count = 0; }
+
+    function createBuffer(): Buffer {
+        return {
+            indices: [],
+            count: 0,
+            squaredDistances: []
+        }
+    }
+
+    /**
+     * Query the tree and store the result to this.buffer. Overwrites the old result.
+     */
+    export function update<T>(ctx: QueryContext<T>, x: number, y: number, z: number, radius: number) {
+        ctx.pivot[0] = x;
+        ctx.pivot[1] = y;
+        ctx.pivot[2] = z;
+        ctx.radius = radius;
+        ctx.radiusSq = radius * radius;
+        resetBuffer(ctx.buffer);
+    }
+
+    export function create<T>(structure: T, mask: Mask, isCheck: boolean): QueryContext<T> {
+        return {
+            structure,
+            buffer: createBuffer(),
+            pivot: [0.1, 0.1, 0.1],
+            radius: 1.1,
+            radiusSq: 1.1 * 1.1,
+            mask,
+            isCheck
+        }
+    }
+}
+
+function createInputData(positions: { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number> }): InputData {
+    const { x, y, z } = positions;
+    const bounds = Box3D.createInfinite();
+    const count = x.length;
+    const { min, max } = bounds;
+
+    for (let i = 0; i < count; i++) {
+        min[0] = Math.min(x[i], min[0]);
+        min[1] = Math.min(y[i], min[1]);
+        min[2] = Math.min(z[i], min[2]);
+        max[0] = Math.max(x[i], max[0]);
+        max[1] = Math.max(y[i], max[1]);
+        max[2] = Math.max(z[i], max[2]);
+    }
+
+    return { positions, bounds, count };
+}
+
+/**
+ * Adapted from https://github.com/arose/ngl
+ * MIT License Copyright (C) 2014+ Alexander Rose
+ */
+
+interface SpatialHash {
+    size: number[],
+    min: number[],
+    grid: Uint32Array,
+    bucketOffset: Uint32Array,
+    bucketCounts: Int32Array,
+    bucketArray: Int32Array,
+    positions: Positions
+}
+
+interface State {
+    size: number[],
+    positions: Positions,
+    bounds: Box3D,
+    count: number
+}
+
+const enum Constants { Exp = 3 }
+
+function nearest(ctx: QueryContext<SpatialHash>): boolean {
+    const { min: [minX, minY, minZ], size: [sX, sY, sZ], bucketOffset, bucketCounts, bucketArray, grid, positions: { x: px, y: py, z: pz } } = ctx.structure;
+    const { radius: r, radiusSq: rSq, pivot: [x, y, z], isCheck, mask } = ctx;
+
+    const loX = Math.max(0, (x - r - minX) >> Constants.Exp);
+    const loY = Math.max(0, (y - r - minY) >> Constants.Exp);
+    const loZ = Math.max(0, (z - r - minZ) >> Constants.Exp);
+
+    const hiX = Math.min(sX, (x + r - minX) >> Constants.Exp);
+    const hiY = Math.min(sY, (y + r - minY) >> Constants.Exp);
+    const hiZ = Math.min(sZ, (z + r - minZ) >> Constants.Exp);
+
+    for (let ix = loX; ix <= hiX; ix++) {
+        for (let iy = loY; iy <= hiY; iy++) {
+            for (let iz = loZ; iz <= hiZ; iz++) {
+                const bucketIdx = grid[(((ix * sY) + iy) * sZ) + iz];
+
+                if (bucketIdx > 0) {
+                    const k = bucketIdx - 1;
+                    const offset = bucketOffset[k];
+                    const count = bucketCounts[k];
+                    const end = offset + count;
+
+                    for (let i = offset; i < end; i++) {
+                        const idx = bucketArray[i];
+                        if (!mask.has(idx)) continue;
+
+                        const dx = px[idx] - x;
+                        const dy = py[idx] - y;
+                        const dz = pz[idx] - z;
+                        const distSq = dx * dx + dy * dy + dz * dz;
+
+                        if (distSq <= rSq) {
+                            if (isCheck) return true;
+                            QueryContext.add(ctx, distSq, idx)
+                        }
+                    }
+                }
+            }
+        }
+    }
+    return ctx.buffer.count > 0;
+}
+
+function _build(state: State): SpatialHash {
+    const { bounds, size: [sX, sY, sZ], positions: { x: px, y: py, z: pz }, count } = state;
+    const n = sX * sY * sZ;
+    const { min: [minX, minY, minZ] } = bounds;
+
+    let bucketCount = 0;
+    const grid = new Uint32Array(n);
+    const bucketIndex = new Int32Array(count);
+    for (let i = 0; i < count; i++) {
+        const x = (px[i] - minX) >> Constants.Exp;
+        const y = (py[i] - minY) >> Constants.Exp;
+        const z = (pz[i] - minZ) >> Constants.Exp;
+        const idx = (((x * sY) + y) * sZ) + z;
+        if ((grid[idx] += 1) === 1) {
+            bucketCount += 1
+        }
+        bucketIndex[i] = idx;
+    }
+
+    const bucketCounts = new Int32Array(bucketCount);
+    for (let i = 0, j = 0; i < n; i++) {
+        const c = grid[i];
+        if (c > 0) {
+            grid[i] = j + 1;
+            bucketCounts[j] = c;
+            j += 1;
+        }
+    }
+
+    const bucketOffset = new Uint32Array(count);
+    for (let i = 1; i < count; ++i) {
+        bucketOffset[i] += bucketOffset[i - 1] + bucketCounts[i - 1];
+    }
+
+    const bucketFill = new Int32Array(bucketCount);
+    const bucketArray = new Int32Array(count);
+    for (let i = 0; i < count; i++) {
+        const bucketIdx = grid[bucketIndex[i]]
+        if (bucketIdx > 0) {
+            const k = bucketIdx - 1;
+            bucketArray[bucketOffset[k] + bucketFill[k]] = i;
+            bucketFill[k] += 1;
+        }
+    }
+
+    return {
+        size: state.size,
+        bucketArray,
+        bucketCounts,
+        bucketOffset,
+        grid,
+        min: state.bounds.min,
+        positions: state.positions
+    }
+}
+
+function build({ positions, bounds, count }: InputData): SpatialHash {
+    const size = [
+        ((bounds.max[0] - bounds.min[0]) >> Constants.Exp) + 1,
+        ((bounds.max[1] - bounds.min[1]) >> Constants.Exp) + 1,
+        ((bounds.max[2] - bounds.min[2]) >> Constants.Exp) + 1
+    ];
+
+    const state: State = {
+        size,
+        positions,
+        bounds,
+        count
+    }
+
+    return _build(state);
+}
+
+export default GridLookup;

src/mol-model/structure/model/formats/mmcif/bonds.ts

@@ -0,0 +1,258 @@
+/**
+ * Copyright (c) 2017-2018 MolQL contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import Model from '../../model'
+import Bonds from '../../properties/bonds'
+import { BondType } from '../../types'
+import { findEntityIdByAsymId, findAtomIndexByLabelName } from './util'
+import { Column } from 'mol-data/db'
+
+export class StructConn implements Bonds.StructConn {
+    private _residuePairIndex: Map<string, StructConn.Entry[]> | undefined = void 0;
+    private _atomIndex: Map<number, StructConn.Entry[]> | undefined = void 0;
+
+    private static _resKey(rA: number, rB: number) {
+        if (rA < rB) return `${rA}-${rB}`;
+        return `${rB}-${rA}`;
+    }
+
+    private getResiduePairIndex() {
+        if (this._residuePairIndex) return this._residuePairIndex;
+        this._residuePairIndex = new Map();
+        for (const e of this.entries) {
+            const ps = e.partners;
+            const l = ps.length;
+            for (let i = 0; i < l - 1; i++) {
+                for (let j = i + i; j < l; j++) {
+                    const key = StructConn._resKey(ps[i].residueIndex, ps[j].residueIndex);
+                    if (this._residuePairIndex.has(key)) {
+                        this._residuePairIndex.get(key)!.push(e);
+                    } else {
+                        this._residuePairIndex.set(key, [e]);
+                    }
+                }
+            }
+        }
+        return this._residuePairIndex;
+    }
+
+    private getAtomIndex() {
+        if (this._atomIndex) return this._atomIndex;
+        this._atomIndex = new Map();
+        for (const e of this.entries) {
+            for (const p of e.partners) {
+                const key = p.atomIndex;
+                if (this._atomIndex.has(key)) {
+                    this._atomIndex.get(key)!.push(e);
+                } else {
+                    this._atomIndex.set(key, [e]);
+                }
+            }
+        }
+        return this._atomIndex;
+    }
+
+    private static _emptyEntry = [];
+
+    getResidueEntries(residueAIndex: number, residueBIndex: number): ReadonlyArray<StructConn.Entry> {
+        return this.getResiduePairIndex().get(StructConn._resKey(residueAIndex, residueBIndex)) || StructConn._emptyEntry;
+    }
+
+    getAtomEntries(atomIndex: number): ReadonlyArray<StructConn.Entry> {
+        return this.getAtomIndex().get(atomIndex) || StructConn._emptyEntry;
+    }
+
+    constructor(public entries: StructConn.Entry[]) {
+    }
+}
+
+export namespace StructConn {
+    export interface Entry extends Bonds.StructConnEntry {
+        distance: number,
+        order: number,
+        flags: number,
+        partners: { residueIndex: number, atomIndex: number, symmetry: string }[]
+    }
+
+    type StructConnType =
+        | 'covale'
+        | 'covale_base'
+        | 'covale_phosphate'
+        | 'covale_sugar'
+        | 'disulf'
+        | 'hydrog'
+        | 'metalc'
+        | 'mismat'
+        | 'modres'
+        | 'saltbr'
+
+    export function create(model: Model): StructConn | undefined {
+        if (model.sourceData.kind !== 'mmCIF') return
+        const { struct_conn } = model.sourceData.data;
+        if (!struct_conn._rowCount) return void 0;
+
+        const { conn_type_id, pdbx_dist_value, pdbx_value_order } = struct_conn;
+        const p1 = {
+            label_asym_id: struct_conn.ptnr1_label_asym_id,
+            label_comp_id: struct_conn.ptnr1_label_comp_id,
+            label_seq_id: struct_conn.ptnr1_label_seq_id,
+            label_atom_id: struct_conn.ptnr1_label_atom_id,
+            label_alt_id: struct_conn.pdbx_ptnr1_label_alt_id,
+            ins_code: struct_conn.pdbx_ptnr1_PDB_ins_code,
+            symmetry: struct_conn.ptnr1_symmetry
+        };
+        const p2: typeof p1 = {
+            label_asym_id: struct_conn.ptnr2_label_asym_id,
+            label_comp_id: struct_conn.ptnr2_label_comp_id,
+            label_seq_id: struct_conn.ptnr2_label_seq_id,
+            label_atom_id: struct_conn.ptnr2_label_atom_id,
+            label_alt_id: struct_conn.pdbx_ptnr2_label_alt_id,
+            ins_code: struct_conn.pdbx_ptnr2_PDB_ins_code,
+            symmetry: struct_conn.ptnr2_symmetry
+        };
+
+        const _p = (row: number, ps: typeof p1) => {
+            if (ps.label_asym_id.valueKind(row) !== Column.ValueKind.Present) return void 0;
+            const asymId = ps.label_asym_id.value(row)
+            const residueIndex = model.hierarchy.findResidueKey(
+                findEntityIdByAsymId(model, asymId),
+                ps.label_comp_id.value(row),
+                asymId,
+                ps.label_seq_id.value(row),
+                ps.ins_code.value(row)
+            );
+            if (residueIndex < 0) return void 0;
+            const atomName = ps.label_atom_id.value(row);
+            // turns out "mismat" records might not have atom name value
+            if (!atomName) return void 0;
+            const atomIndex = findAtomIndexByLabelName(model, residueIndex, atomName, ps.label_alt_id.value(row));
+            if (atomIndex < 0) return void 0;
+            return { residueIndex, atomIndex, symmetry: ps.symmetry.value(row) || '1_555' };
+        }
+
+        const _ps = (row: number) => {
+            const ret = [];
+            let p = _p(row, p1);
+            if (p) ret.push(p);
+            p = _p(row, p2);
+            if (p) ret.push(p);
+            return ret;
+        }
+
+        const entries: StructConn.Entry[] = [];
+        for (let i = 0; i < struct_conn._rowCount; i++) {
+            const partners = _ps(i);
+            if (partners.length < 2) continue;
+
+            const type = conn_type_id.value(i)! as StructConnType;
+            const orderType = (pdbx_value_order.value(i) || '').toLowerCase();
+            let flags = BondType.Flag.None;
+            let order = 1;
+
+            switch (orderType) {
+                case 'sing': order = 1; break;
+                case 'doub': order = 2; break;
+                case 'trip': order = 3; break;
+                case 'quad': order = 4; break;
+            }
+
+            switch (type) {
+                case 'covale':
+                case 'covale_base':
+                case 'covale_phosphate':
+                case 'covale_sugar':
+                case 'modres':
+                    flags = BondType.Flag.Covalent;
+                    break;
+                case 'disulf': flags = BondType.Flag.Covalent | BondType.Flag.Sulfide; break;
+                case 'hydrog': flags = BondType.Flag.Hydrogen; break;
+                case 'metalc': flags = BondType.Flag.MetallicCoordination; break;
+                case 'saltbr': flags = BondType.Flag.Ion; break;
+            }
+
+            entries.push({ flags, order, distance: pdbx_dist_value.value(i), partners });
+        }
+
+        return new StructConn(entries);
+    }
+}
+
+export class ComponentBondInfo implements Bonds.ComponentBondInfo {
+    entries: Map<string, ComponentBondInfo.Entry> = new Map();
+
+    newEntry(id: string) {
+        let e = new ComponentBondInfo.Entry(id);
+        this.entries.set(id, e);
+        return e;
+    }
+}
+
+export namespace ComponentBondInfo {
+    export class Entry implements Bonds.ComponentBondInfoEntry {
+        map: Map<string, Map<string, { order: number, flags: number }>> = new Map();
+
+        add(a: string, b: string, order: number, flags: number, swap = true) {
+            let e = this.map.get(a);
+            if (e !== void 0) {
+                let f = e.get(b);
+                if (f === void 0) {
+                    e.set(b, { order, flags });
+                }
+            } else {
+                let map = new Map<string, { order: number, flags: number }>();
+                map.set(b, { order, flags });
+                this.map.set(a, map);
+            }
+
+            if (swap) this.add(b, a, order, flags, false);
+        }
+
+        constructor(public id: string) {
+        }
+    }
+
+    export function create(model: Model): ComponentBondInfo | undefined {
+        if (model.sourceData.kind !== 'mmCIF') return
+        const { chem_comp_bond } = model.sourceData.data;
+        if (!chem_comp_bond._rowCount) return void 0;
+
+        let info = new ComponentBondInfo();
+
+        const { comp_id, atom_id_1, atom_id_2, value_order, pdbx_aromatic_flag, _rowCount: rowCount } = chem_comp_bond;
+
+        let entry = info.newEntry(comp_id.value(0)!);
+
+        for (let i = 0; i < rowCount; i++) {
+
+            const id = comp_id.value(i)!;
+            const nameA = atom_id_1.value(i)!;
+            const nameB = atom_id_2.value(i)!;
+            const order = value_order.value(i)!;
+            const aromatic = pdbx_aromatic_flag.value(i) === 'Y';
+
+            if (entry.id !== id) {
+                entry = info.newEntry(id);
+            }
+
+            let flags: number = BondType.Flag.Covalent;
+            let ord = 1;
+            if (aromatic) flags |= BondType.Flag.Aromatic;
+            switch (order.toLowerCase()) {
+                case 'doub':
+                case 'delo':
+                    ord = 2;
+                    break;
+                case 'trip': ord = 3; break;
+                case 'quad': ord = 4; break;
+            }
+
+            entry.add(nameA, nameB, ord, flags);
+        }
+
+        return info;
+    }
+}

src/mol-model/structure/model/formats/mmcif/util.ts

@@ -0,0 +1,26 @@
+/**
+ * Copyright (c) 2017-2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import Model from '../../model'
+
+export function findEntityIdByAsymId(model: Model, asymId: string) {
+    if (model.sourceData.kind !== 'mmCIF') return ''
+    const { struct_asym } = model.sourceData.data
+    for (let i = 0, n = struct_asym._rowCount; i < n; ++i) {
+        if (struct_asym.id.value(i) === asymId) return struct_asym.entity_id.value(i)
+    }
+    return ''
+}
+
+export function findAtomIndexByLabelName(model: Model, residueIndex: number, atomName: string, altLoc: string | null) {
+    const { segmentMap, segments } = model.hierarchy.residueSegments
+    const idx = segmentMap[residueIndex]
+    const { label_atom_id, label_alt_id } = model.hierarchy.atoms;
+    for (let i = segments[idx], n = segments[idx + 1]; i <= n; ++i) {
+        if (label_atom_id.value(i) === atomName && (!altLoc || label_alt_id.value(i) === altLoc)) return i;
+    }
+    return -1;
+}

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

@@ -5,10 +5,14 @@
  */
 
 import UUID from 'mol-util/uuid'
+import GridLookup from 'mol-math/geometry/grid-lookup'
 import Format from './format'
 import Hierarchy from './properties/hierarchy'
 import Conformation from './properties/conformation'
 import Symmetry from './properties/symmetry'
+import Bonds from './properties/bonds'
+
+import computeBonds from './utils/compute-bonds'
 
 import from_gro from './formats/gro'
 import from_mmCIF from './formats/mmcif'
@@ -30,8 +34,11 @@ interface Model extends Readonly<{
     conformation: Conformation,
     symmetry: Symmetry,
 
-    atomCount: number
-}> { }
+    atomCount: number,
+}> {
+    '@spatialLookup'?: GridLookup,
+    '@bonds'?: Bonds
+} { }
 
 namespace Model {
     export function create(format: Format) {
@@ -40,6 +47,18 @@ namespace Model {
             case 'mmCIF': return from_mmCIF(format);
         }
     }
+    export function spatialLookup(model: Model): GridLookup {
+        if (model['@spatialLookup']) return model['@spatialLookup']!;
+        const lookup = GridLookup(model.conformation);
+        model['@spatialLookup'] = lookup;
+        return lookup;
+    }
+    export function bonds(model: Model): Bonds {
+        if (model['@bonds']) return model['@bonds']!;
+        const bonds = computeBonds(model);
+        model['@bonds'] = bonds;
+        return bonds;
+    }
 }
 
 export default Model

src/mol-model/structure/model/properties/bonds.ts

@@ -0,0 +1,51 @@
+/**
+ * Copyright (c) 2017-2018 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>
+ */
+
+import { BondType } from '../types'
+
+interface Bonds {
+    /**
+     * Where bonds for atom A start and end.
+     * Start offset at idx, end at idx + 1
+     */
+    offset: ArrayLike<number>,
+    neighbor: ArrayLike<number>,
+
+    order: ArrayLike<number>,
+    flags: ArrayLike<BondType.Flag>,
+
+    count: number
+}
+
+namespace Bonds {
+    export function createEmpty(): Bonds {
+        return { offset: [], neighbor: [], order: [], flags: [], count: 0 }
+    }
+    export function isCovalent(flags: number) {
+        return (flags & BondType.Flag.Covalent) !== 0;
+    }
+    export interface StructConnEntry {
+        flags: BondType.Flag,
+        order: number,
+        distance: number,
+        partners: { residueIndex: number, atomIndex: number, symmetry: string }[]
+    }
+    export interface StructConn {
+        getResidueEntries(residueAIndex: number, residueBIndex: number): ReadonlyArray<StructConnEntry>
+        getAtomEntries(atomIndex: number): ReadonlyArray<StructConnEntry>
+    }
+    export interface ComponentBondInfoEntry {
+        map: Map<string, Map<string, { order: number, flags: number }>>
+    }
+    export interface ComponentBondInfo {
+        entries: Map<string, ComponentBondInfoEntry>
+    }
+}
+
+
+
+export default Bonds

src/mol-model/structure/model/properties/seconday-structure.ts

@@ -0,0 +1,14 @@
+/**
+ * Copyright (c) 2017 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+/** Secondary structure "indexed" by residues. */
+export interface SecondaryStructure {
+    type: number[],
+    index: number[],
+    flags: number[],
+    /** unique value for each "element". This is because single sheet is speficied by multiple records. */
+    key: number[]
+}

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

@@ -337,4 +337,20 @@ export const VdwRadii = {
     'LV': 2.0,
     'UUH': 2.0
 }
-export const DefaultVdwRadius = 2.0
+export const DefaultVdwRadius = 2.0
+
+export interface BondType extends BitFlags<BondType.Flag> { }
+export namespace BondType {
+    export const is: (b: BondType, f: Flag) => boolean = BitFlags.has
+    export const enum Flag {
+        None                 = 0x0,
+        Covalent             = 0x1,
+        MetallicCoordination = 0x2,
+        Hydrogen             = 0x4,
+        Ion                  = 0x8,
+        Sulfide              = 0x10,
+        Aromatic             = 0x20,
+        Computed             = 0x40
+        // currently at most 16 flags are supported!!
+    }
+}

src/mol-util/mask.ts

@@ -0,0 +1,184 @@
+/*
+ * Copyright (c) 2017 MolQL contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+// TODO check if the removal of FastSet and the removal of the context object for forEach
+// have any performance implications
+
+function _ascSort(a: number, b: number) {
+    return a - b;
+}
+
+export function sortAsc<T extends ArrayLike<number>>(array: T): T {
+    Array.prototype.sort.call(array, _ascSort);
+    return array;
+}
+
+interface Mask {
+    '@type': 'mask'
+    size: number;
+    has(i: number): boolean;
+    /** in-order iteration of all "masked elements". */
+    forEach<Ctx>(f: (i: number, ctx?: Ctx) => void, ctx?: Ctx): Ctx | undefined;
+}
+
+namespace Mask {
+    class EmptyMask implements Mask {
+        '@type': 'mask'
+        size = 0;
+        has(i: number) { return false; }
+        forEach<Ctx>(f: (i: number, ctx?: Ctx) => void, ctx: Ctx) { return ctx; }
+        constructor() { }
+    }
+
+    class SingletonMask implements Mask {
+        '@type': 'mask'
+        size = 1;
+        has(i: number) { return i === this.idx; }
+        forEach<Ctx>(f: (i: number, ctx?: Ctx) => void, ctx?: Ctx) { f(this.idx, ctx); return ctx; }
+        constructor(private idx: number) { }
+    }
+
+    class BitMask implements Mask {
+        '@type': 'mask'
+        private length: number;
+        has(i: number) { return i < this.length && !!this.mask[i] as any; }
+
+        private _forEach<Ctx>(f: (i: number, ctx?: Ctx) => void, ctx: Ctx | undefined) {
+            for (let i = 0; i < this.length; i++) {
+                if (this.mask[i]) f(i, ctx);
+            }
+        }
+        forEach<Ctx>(f: (i: number, ctx?: Ctx) => void, ctx?: Ctx) {
+            this._forEach(f, ctx);
+            return ctx;
+        }
+        constructor(private mask: boolean[], public size: number) { this.length = mask.length;  }
+    }
+
+    class AllMask implements Mask {
+        '@type': 'mask'
+        has(i: number) { return true; }
+        private _forEach<Ctx>(f: (i: number, ctx?: Ctx) => void, ctx: Ctx | undefined) {
+            for (let i = 0; i < this.size; i++) {
+                f(i, ctx);
+            }
+        }
+        forEach<Ctx>(f: (i: number, ctx?: Ctx) => void, ctx?: Ctx) {
+            this._forEach(f, ctx);
+            return ctx;
+        }
+        constructor(public size: number) { }
+    }
+
+    class SetMask implements Mask {
+        '@type': 'mask'
+        private _flat: number[] | undefined = void 0;
+        size: number;
+        has(i: number) { return this.set.has(i); }
+
+        private _forEach<Ctx>(f: (i: number, ctx: Ctx) => void, ctx: Ctx) {
+            for (const idx of this.flatten()) {
+                f(idx, ctx);
+            }
+        }
+        private flatten() {
+            if (this._flat) return this._flat;
+            const indices = new Int32Array(this.size);
+            let offset = 0
+            this.set.forEach(i => indices[offset++] = i);
+            sortAsc(indices);
+            this._flat = indices as any as number[];
+            return this._flat;
+        }
+        forEach<Ctx>(f: (i: number, ctx: Ctx) => void, ctx: Ctx) {
+            this._forEach(f, ctx);
+            return ctx;
+        }
+        constructor(private set: Set<number>) {
+            this.size = set.size;
+        }
+    }
+
+    export function always(size: number) { return new AllMask(size); }
+    export const never = new EmptyMask();
+
+    export function ofSet(set: Set<number>): Mask {
+        return new SetMask(set);
+    }
+
+    export function singleton(i: number) {
+        return new SingletonMask(i);
+    }
+
+    export function ofUniqueIndices(indices: ArrayLike<number>): Mask {
+        const len = indices.length;
+        if (len === 0) return new EmptyMask();
+        if (len === 1) return new SingletonMask(indices[0]);
+
+        let max = 0;
+        for (const i of (indices as number[])) {
+            if (i > max) max = i;
+        }
+        if (len === max) return new AllMask(len);
+
+        const f = len / max;
+        if (f < 1 / 12) {
+            const set = new Set<number>();
+            for (const i of (indices as number[])) set.add(i);
+            return new SetMask(set);
+        }
+
+        const mask = new Int8Array(max + 1);
+        for (const i of (indices as number[])) {
+            mask[i] = 1;
+        }
+        return new BitMask(mask as any as boolean[], indices.length);
+    }
+
+    export function ofMask(mask: boolean[], size: number): Mask {
+        return new BitMask(mask, size);
+    }
+
+    export function hasAny(mask: Mask, xs: number[]) {
+        for (const x of xs) {
+            if (mask.has(x)) return true;
+        }
+        return false;
+    }
+
+    export function complement(mask: Mask, against: Mask) {
+        let count = 0
+        let max = 0
+        against.forEach(i => {
+            if (!mask.has(i)) {
+                count++;
+                if (i > max) max = i;
+            }
+        });
+
+        if (count / max < 1 / 12) {
+            // set based
+            const set = new Set<number>()
+            against.forEach(i => {
+                if (!mask.has(i)) {
+                    set.add(i);
+                }
+            });
+            return ofSet(set);
+        } else {
+            // mask based
+            const target = new Uint8Array(max + 1);
+            against.forEach(i => {
+                if (!mask.has(i)) {
+                    target[i] = 1;
+                }
+            });
+            return ofMask(target as any as boolean[], count);
+        }
+    }
+}
+
+export default Mask;

src/script.ts

@@ -1,207 +0,0 @@
-/**
- * Copyright (c) 2017 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Alexander Rose <alexander.rose@weirdbyte.de>
- * @author David Sehnal <david.sehnal@gmail.com>
- */
-
-import * as util from 'util'
-import * as fs from 'fs'
-
-require('util.promisify').shim();
-const readFileAsync = util.promisify(fs.readFile);
-
-import Gro from 'mol-io/reader/gro/parser'
-import Csv from 'mol-io/reader/csv/parser'
-import CIF from 'mol-io/reader/cif'
-
-import Computation from 'mol-util/computation'
-
-import { Model } from 'mol-model/structure'
-
-const file = '1crn.gro'
-// const file = 'water.gro'
-// const file = 'test.gro'
-// const file = 'md_1u19_trj.gro'
-
-function showProgress(tag: string, p: Computation.Progress) {
-    console.log(`[${tag}] ${p.message} ${p.isIndeterminate ? '' : (p.current / p.max * 100).toFixed(2) + '% '}(${p.elapsedMs | 0}ms)`)
-}
-
-async function runGro(input: string) {
-    console.time('parseGro');
-    const comp = Gro(input);
-
-    const ctx = Computation.observable({ updateRateMs: 150, observer: p => showProgress('GRO', p) });
-    const parsed = await comp(ctx);
-    console.timeEnd('parseGro');
-
-    if (parsed.isError) {
-        console.log(parsed);
-        return;
-    }
-
-    const groFile = parsed.result
-
-    console.log('structure count: ', groFile.structures.length);
-
-    const data = groFile.structures[0];
-
-    // const header = groFile.blocks[0].getCategory('header')
-    const { header, atoms } = data;
-    console.log(JSON.stringify(header, null, 2));
-    console.log('number of atoms:', atoms.count);
-
-    console.log(`'${atoms.residueNumber.value(1)}'`)
-    console.log(`'${atoms.residueName.value(1)}'`)
-    console.log(`'${atoms.atomName.value(1)}'`)
-    console.log(atoms.z.value(1))
-    console.log(`'${atoms.z.value(1)}'`)
-
-    const n = atoms.count;
-    console.log('rowCount', n)
-
-    console.time('getFloatArray x')
-    const x = atoms.x.toArray({ array: Float32Array })
-    console.timeEnd('getFloatArray x')
-    console.log(x.length, x[0], x[x.length - 1])
-
-    console.time('getFloatArray y')
-    const y = atoms.y.toArray({ array: Float32Array })
-    console.timeEnd('getFloatArray y')
-    console.log(y.length, y[0], y[y.length - 1])
-
-    console.time('getFloatArray z')
-    const z = atoms.z.toArray({ array: Float32Array })
-    console.timeEnd('getFloatArray z')
-    console.log(z.length, z[0], z[z.length - 1])
-
-    console.time('getIntArray residueNumber')
-    const residueNumber = atoms.residueNumber.toArray({ array: Int32Array })
-    console.timeEnd('getIntArray residueNumber')
-    console.log(residueNumber.length, residueNumber[0], residueNumber[residueNumber.length - 1])
-}
-
-export async function _gro() {
-    const input = await readFileAsync(`./examples/${file}`, 'utf8')
-    runGro(input)
-}
-
-// _gro()
-
-async function runCIF(input: string | Uint8Array) {
-    console.time('parseCIF');
-    const comp = typeof input === 'string' ? CIF.parseText(input) : CIF.parseBinary(input);
-
-    const ctx = Computation.observable({ updateRateMs: 250, observer: p => showProgress('CIF', p) });
-    const parsed = await comp(ctx);
-    console.timeEnd('parseCIF');
-    if (parsed.isError) {
-        console.log(parsed);
-        return;
-    }
-
-    const data = parsed.result.blocks[0];
-    const atom_site = data.categories._atom_site;
-    console.log(atom_site.getField('Cartn_x')!.float(0));
-    //console.log(atom_site.getField('label_atom_id')!.toStringArray());
-
-    const mmcif = CIF.schema.mmCIF(data);
-    console.log(mmcif.atom_site.Cartn_x.value(0));
-    console.log(mmcif.entity.type.toArray());
-    console.log(mmcif.pdbx_struct_oper_list.matrix.value(0));
-
-    console.time('createModels');
-    const models = Model.create({ kind: 'mmCIF', data: mmcif });
-    console.timeEnd('createModels');
-
-    for (let i = 0; i < models.length; i++) {
-        console.log(models[i].id, models[i].conformation.id);
-    }
-
-    // console.log(models[0].hierarchy.isMonotonous);
-    // console.log(models[0].hierarchy.atoms.type_symbol.value(0));
-    // console.log(models[0].hierarchy.residues.auth_comp_id.value(0));
-    // console.log(models[0].hierarchy.residues.auth_comp_id.value(1));
-    // console.log(models[0].hierarchy.chains.auth_asym_id.value(0));
-    // console.log(models[0].hierarchy.chains.auth_asym_id.value(1));
-    // console.log(models[0].hierarchy.chains.label_asym_id.value(1));
-    // console.log(models[0].conformation.x[0]);
-    // console.log(models[0].conformation.y[0]);
-    // console.log(models[0].conformation.z[0]);
-
-    // const schema = await _dic()
-    // if (schema) {
-    //     const mmcif2 = applySchema(schema, data)
-    //     // console.log(util.inspect(mmcif2.atom_site, {showHidden: false, depth: 3}))
-    //     console.log(mmcif2.atom_site.Cartn_x.value(0));
-    //     console.log(mmcif2.entity.type.toArray());
-    //     // console.log(mmcif2.pdbx_struct_oper_list.matrix.value(0)); // TODO
-    // } else {
-    //     console.log('error getting mmcif schema from dic')
-    // }
-}
-
-export async function _cif() {
-    let path = `./examples/1grm_updated.cif`;
-    // path = '../test/3j3q.cif'  // lets have a relative path for big test files
-    // path = 'e:/test/quick/3j3q_updated.cif';
-    const input = await readFileAsync(path, 'utf8')
-    console.log('------------------');
-    console.log('Text CIF:');
-    runCIF(input);
-
-    path = `./examples/1cbs_full.bcif`;
-
-    const input2 = await readFileAsync(path)
-    console.log('------------------');
-    console.log('BinaryCIF:');
-    const data = new Uint8Array(input2.byteLength);
-    for (let i = 0; i < input2.byteLength; i++) data[i] = input2[i];
-    runCIF(input2);
-}
-
-// _cif();
-
-const comp = Computation.create(async ctx => {
-    for (let i = 0; i < 0; i++) {
-        await new Promise(res => setTimeout(res, 500));
-        if (ctx.requiresUpdate) await ctx.update({ message: 'working', current: i, max: 2 });
-    }
-    return 42;
-});
-export async function testComp() {
-    const ctx = Computation.observable({ observer: p => showProgress('test', p) });
-    const ret = await comp(ctx);
-    console.log('computation returned', ret);
-}
-// testComp();
-
-
-const csvString = ` Year,Make,Model,Length
-1997,Ford,"E350
-
-MOIN",2.34
-2000,Mercury, Cougar,2.38`
-
-export async function testCsv () {
-    const parsed = await Csv(csvString)();
-
-    if (parsed.isError) {
-        console.log(parsed)
-        return;
-    }
-
-    const csvFile = parsed.result;
-    csvFile.table.columnNames.forEach(name => {
-        const col = csvFile.table.getColumn(name)
-        if (col) console.log(name, col.toStringArray())
-    })
-
-    const year = csvFile.table.getColumn('Year')
-    if (year) console.log('(int)Year', year.toIntArray())
-
-    const length = csvFile.table.getColumn('Length')
-    if (length) console.log('(float)Length', length.toFloatArray())
-}
-testCsv()