SymmetryOperator in mol-math

src/mol-math/geometry/symmetry-operator.ts

+/**
+ * Copyright (c) 2017 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import { Vec3, Mat4 } from '../linear-algebra/3d'
+
+interface SymmetryOperator {
+    readonly name: string,
+    readonly hkl: Vec3,
+
+    readonly matrix: Mat4,
+    // cache the inverse of the transform
+    readonly inverse: Mat4,
+    // optimize the identity case
+    readonly isIdentity: boolean
+}
+
+namespace SymmetryOperator {
+    export const DefaultName = '1_555'
+    export const Default: SymmetryOperator = create(DefaultName, Mat4.identity());
+
+    export function create(name: string, matrix: Mat4, hkl?: number[]): SymmetryOperator {
+        const _hkl = hkl ? Vec3.create(hkl[0], hkl[1], hkl[2]) : Vec3.zero();
+        if (Mat4.isIdentity(matrix)) return { name, matrix, inverse: Mat4.identity(), isIdentity: true, hkl: _hkl };
+        if (!Mat4.isRotationAndTranslation(matrix)) throw new Error('Symmetry operator must be a composition of rotation and translation.');
+        return { name, matrix, inverse: Mat4.invert(Mat4.zero(), matrix), isIdentity: false, hkl: _hkl };
+    }
+
+    export interface CoordinateMapper { (index: number, slot: Vec3): Vec3 }
+    export interface ArrayMapping {
+        readonly invariantPosition: CoordinateMapper,
+        readonly position: CoordinateMapper,
+        x(index: number): number,
+        y(index: number): number,
+        z(index: number): number
+    }
+
+    export interface Coordinates { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number> }
+
+    export function createMapping(operator: SymmetryOperator, coords: Coordinates) {
+        const invariantPosition = SymmetryOperator.createCoordinateMapper(SymmetryOperator.Default, coords);
+        const position = operator.isIdentity ? invariantPosition : SymmetryOperator.createCoordinateMapper(operator, coords);
+        const { x, y, z } = createProjections(operator, coords);
+        return { invariantPosition, position, x, y, z };
+    }
+
+    export function createCoordinateMapper(t: SymmetryOperator, coords: Coordinates): CoordinateMapper {
+        if (t.isIdentity) return identityPosition(coords);
+        return generalPosition(t, coords);
+    }
+}
+
+export default SymmetryOperator
+
+interface Projections { x(index: number): number, y(index: number): number, z(index: number): number }
+
+function createProjections(t: SymmetryOperator, coords: SymmetryOperator.Coordinates): Projections {
+    if (t.isIdentity) return { x: projectCoord(coords.x), y: projectCoord(coords.y), z: projectCoord(coords.z) };
+    return { x: projectX(t, coords), y: projectY(t, coords), z: projectZ(t, coords) };
+}
+
+function projectCoord(xs: ArrayLike<number>) {
+    return (i: number) => xs[i];
+}
+
+function isW1(m: Mat4) {
+    return m[3] === 0 && m[7] === 0 && m[11] === 0 && m[15] === 1;
+}
+
+function projectX({ matrix: m }: SymmetryOperator, { x: xs, y: ys, z: zs}: SymmetryOperator.Coordinates) {
+    const xx = m[0], yy = m[4], zz = m[8], tx = m[12];
+
+    if (isW1(m)) {
+        // this should always be the case.
+        return (i: number) => xx * xs[i] + yy * ys[i] + zz * zs[i] + tx;
+    }
+
+    return (i: number) => {
+        const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
+        return (xx * x + yy * y + zz * z + tx) / w;
+    }
+}
+
+function projectY({ matrix: m }: SymmetryOperator, { x: xs, y: ys, z: zs}: SymmetryOperator.Coordinates) {
+    const xx = m[1], yy = m[5], zz = m[9], ty = m[13];
+
+    if (isW1(m)) {
+        // this should always be the case.
+        return (i: number) => xx * xs[i] + yy * ys[i] + zz * zs[i] + ty;
+    }
+
+    return (i: number) => {
+        const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
+        return (xx * x + yy * y + zz * z + ty) / w;
+    }
+}
+
+function projectZ({ matrix: m }: SymmetryOperator, { x: xs, y: ys, z: zs}: SymmetryOperator.Coordinates) {
+    const xx = m[2], yy = m[6], zz = m[10], tz = m[14];
+
+    if (isW1(m)) {
+        // this should always be the case.
+        return (i: number) => xx * xs[i] + yy * ys[i] + zz * zs[i] + tz;
+    }
+
+    return (i: number) => {
+        const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
+        return (xx * x + yy * y + zz * z + tz) / w;
+    }
+}
+
+function identityPosition({ x, y, z }: SymmetryOperator.Coordinates): SymmetryOperator.CoordinateMapper {
+    return (i, s) => {
+        s[0] = x[i];
+        s[1] = y[i];
+        s[2] = z[i];
+        return s;
+    }
+}
+
+function generalPosition({ matrix: m }: SymmetryOperator, { x: xs, y: ys, z: zs }: SymmetryOperator.Coordinates) {
+    if (isW1(m)) {
+        // this should always be the case.
+        return (i: number, r: Vec3): Vec3 => {
+            const x = xs[i], y = ys[i], z = zs[i];
+            r[0] = m[0] * x + m[4] * y + m[8] * z + m[12];
+            r[1] = m[1] * x + m[5] * y + m[9] * z + m[13];
+            r[2] = m[2] * x + m[6] * y + m[10] * z + m[14];
+            return r;
+        }
+    }
+    return (i: number, r: Vec3): Vec3 => {
+        r[0] = xs[i];
+        r[1] = ys[i];
+        r[2] = zs[i];
+        Vec3.transformMat4(r, r, m);
+        return r;
+    }
+}
\ No newline at end of file

src/mol-math/geometry/transform.ts

-/**
- * Copyright (c) 2017 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author David Sehnal <david.sehnal@gmail.com>
- */
-
-import { Mat4 } from '../linear-algebra/3d'
-
-interface Transform extends Readonly<{
-    transform: Mat4,
-    // cache the inverse of the transform
-    inverse: Mat4,
-    // optimize the identity case
-    isIdentity: boolean
-}> { }
-
-namespace Transform {
-    export function isIdentity(m: Mat4) {
-        throw 'nyi'
-    }
-
-    export function isRotationAndTranslation(m: Mat4) {
-        throw 'nyi'
-    }
-}
-
-export default Transform
\ No newline at end of file

src/mol-math/linear-algebra/3d.ts

@@ -152,7 +152,8 @@ export namespace Mat4 {
         let det = b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
 
         if (!det) {
-            return Number.NaN;
+            console.warn('non-invertible matrix.', a);
+            return out;
         }
         det = 1.0 / det;
 
@@ -430,6 +431,18 @@ export namespace Mat4 {
         // Calculate the determinant
         return b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
     }
+
+    export function isRotationAndTranslation(a: Mat4) {
+        const a00 = a[0], a01 = a[1], a02 = a[2], a03 = a[3],
+            a10 = a[4], a11 = a[5], a12 = a[6], a13 = a[7],
+            a20 = a[8], a21 = a[9], a22 = a[10], a23 = a[11],
+            /* a30 = a[12], a31 = a[13], a32 = a[14],*/ a33 = a[15];
+
+        if (a33 !== 1 || a03 !== 0 || a13 !== 0 || a23 !== 0) return false;
+        const det3x3 = a00 * (a11 * a22 - a21 * a21) - a01 * (a10 * a22 - a12 * a20) + a02 * (a10 * a21 - a11 * a20);
+        if (det3x3 < 1 - EPSILON.Value || det3x3 > 1 + EPSILON.Value) return false;
+        return true;
+    }
 }
 
 export namespace Vec3 {

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

@@ -70,9 +70,9 @@ function getConformation({ data }: mmCIF_Format, bounds: Interval): Conformation
         atomId: Column.window(atom_site.id, start, end),
         occupancy: Column.window(atom_site.occupancy, start, end),
         B_iso_or_equiv: Column.window(atom_site.B_iso_or_equiv, start, end),
-        __x: atom_site.Cartn_x.toArray({ array: Float32Array, start, end }),
-        __y: atom_site.Cartn_y.toArray({ array: Float32Array, start, end }),
-        __z: atom_site.Cartn_z.toArray({ array: Float32Array, start, end }),
+        x: atom_site.Cartn_x.toArray({ array: Float32Array, start, end }),
+        y: atom_site.Cartn_y.toArray({ array: Float32Array, start, end }),
+        z: atom_site.Cartn_z.toArray({ array: Float32Array, start, end }),
     }
 }
 

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

@@ -20,9 +20,9 @@ interface Conformation {
 
     // Coordinates. Generally, not to be accessed directly because the coordinate might be
     // transformed by an operator. Use Unit.getPosition instead.
-    __x: ArrayLike<number>,
-    __y: ArrayLike<number>,
-    __z: ArrayLike<number>
+    x: ArrayLike<number>,
+    y: ArrayLike<number>,
+    z: ArrayLike<number>
 }
 
 export default Conformation
\ No newline at end of file

src/mol-model/structure/structure.ts

@@ -7,7 +7,6 @@
 import Atom from './structure/atom'
 import AtomSet from './structure/atom/set'
 import Structure from './structure/structure'
-import Operator from './structure/operator'
 import Unit from './structure/unit'
 
-export { Atom, AtomSet, Structure, Operator, Unit }
+export { Atom, AtomSet, Structure, Unit }
\ No newline at end of file

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

-/**
- * Copyright (c) 2017 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author David Sehnal <david.sehnal@gmail.com>
- */
-
-import { Mat4 } from 'mol-math/linear-algebra'
-
-interface Operator extends Readonly<{
-    name: string,
-    hkl: number[], // defaults to [0, 0, 0] for non symmetry entries
-    matrix: Mat4,
-    // cache the inverse of the transform
-    inverse: Mat4,
-    // optimize the identity case
-    isIdentity: boolean
-}> { }
-
-namespace Operator {
-    export const Identity: Operator = { name: '1_555', hkl: [0, 0, 0], matrix: Mat4.identity(), inverse: Mat4.identity(), isIdentity: true };
-}
-
-export default Operator
\ No newline at end of file

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

@@ -6,9 +6,9 @@
 
 import { OrderedSet, Iterator } from 'mol-data/int'
 import UniqueArray from 'mol-data/util/unique-array'
+import SymmetryOperator from 'mol-math/geometry/symmetry-operator'
 import { Model, Format } from '../model'
 import Unit from './unit'
-import Operator from './operator'
 import AtomSet from './atom/set'
 import Atom from './atom'
 
@@ -35,7 +35,7 @@ namespace Structure {
         const builder = Builder();
 
         for (let c = 0; c < chains.count; c++) {
-            const unit = Unit.create(model, Operator.Identity);
+            const unit = Unit.create(model, SymmetryOperator.Default);
             builder.addUnit(unit);
             builder.addAtoms(unit.id, OrderedSet.ofBounds(chains.segments[c], chains.segments[c + 1]));
         }

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

@@ -4,45 +4,34 @@
  * @author David Sehnal <david.sehnal@gmail.com>
  */
 
+import SymmetryOperator from 'mol-math/geometry/symmetry-operator'
 import { Model } from '../model'
-import Operator from './operator'
-import { Vec3, Mat4 } from 'mol-math/linear-algebra'
 
-interface Unit extends Readonly<{
+interface Unit extends SymmetryOperator.ArrayMapping {
     // Structure-level unique identifier of the unit.
-    id: number,
+    readonly id: number,
 
     // Provides access to the underlying data.
-    model: Model,
+    readonly model: Model,
 
     // Determines the operation applied to this unit.
     // The transform and and inverse are baked into the "getPosition" function
-    operator: Operator,
+    readonly operator: SymmetryOperator,
 
     // Reference some commonly accessed things for faster access.
-    residueIndex: ArrayLike<number>,
-    chainIndex: ArrayLike<number>,
-    hierarchy: Model['hierarchy'],
-    conformation: Model['conformation']
-}> {
-    // returns the untransformed position. Used for spatial queries.
-    getInvariantPosition(atom: number, slot: Vec3): Vec3,
+    readonly residueIndex: ArrayLike<number>,
+    readonly chainIndex: ArrayLike<number>,
+    readonly hierarchy: Model['hierarchy'],
+    readonly conformation: Model['conformation']
 
-    // gets the transformed position of the specified atom.
-    getPosition(atom: number, slot: Vec3): Vec3,
-
-    // optimized x/y/z coordinate projections for your query needs.
-    x(atom: number): number,
-    y(atom: number): number,
-    z(atom: number): number
+    // TODO: add velocity?
 }
 
 namespace Unit {
-    export function create(model: Model, operator: Operator): Unit {
+    export function create(model: Model, operator: SymmetryOperator): Unit {
         const h = model.hierarchy;
-        const { __x: x, __y: y, __z: z } = model.conformation;
-        const getInvariantPosition = makeGetInvariantPosition(x, y, z);
-        const getPosition = operator.isIdentity ? getInvariantPosition : makeGetPosition(operator.matrix, x, y, z);
+        const { invariantPosition, position, x, y, z } = SymmetryOperator.createMapping(operator, model.conformation);
+
         return {
             id: nextUnitId(),
             model,
@@ -51,11 +40,9 @@ namespace Unit {
             chainIndex: h.chainSegments.segmentMap,
             hierarchy: model.hierarchy,
             conformation: model.conformation,
-            getInvariantPosition,
-            getPosition,
-            x: operator.isIdentity ? makeGetCoord(x) : makeX(operator.matrix, x, y, z),
-            y: operator.isIdentity ? makeGetCoord(y) : makeY(operator.matrix, x, y, z),
-            z: operator.isIdentity ? makeGetCoord(z) : makeZ(operator.matrix, x, y, z)
+            invariantPosition,
+            position,
+            x, y, z
         };
     }
 }
@@ -67,73 +54,4 @@ function nextUnitId() {
     const ret = _id;
     _id = (_id + 1) % 0x3fffffff;
     return ret;
-}
-
-function makeGetInvariantPosition(x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number>) {
-    return (a: number, r: Vec3): Vec3 => {
-        r[0] = x[a];
-        r[1] = y[a];
-        r[2] = z[a];
-        return r;
-    }
-}
-
-function makeGetCoord(xs: ArrayLike<number>) {
-    return (i: number) => xs[i];
-}
-
-function isWConst(m: Mat4) {
-    return m[3] === 0 && m[7] === 0 && m[11] === 0;
-}
-
-function makeX(m: Mat4, xs: ArrayLike<number>, ys: ArrayLike<number>, zs: ArrayLike<number>) {
-    const xx = m[0], yy = m[4], zz = m[8], ww = m[12];
-
-    if (isWConst(m)) {
-        const w = m[15] || 1.0;
-        return (i: number) => (xx * xs[i] + yy * ys[i] + zz * zs[i] + ww) / w;
-    }
-
-    return (i: number) => {
-        const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
-        return (xx * x + yy * y + zz * z + ww) / w;
-    }
-}
-
-function makeY(m: Mat4, xs: ArrayLike<number>, ys: ArrayLike<number>, zs: ArrayLike<number>) {
-    const xx = m[1], yy = m[5], zz = m[9], ww = m[13];
-
-    if (isWConst(m)) {
-        const w = m[15] || 1.0;
-        return (i: number) => (xx * xs[i] + yy * ys[i] + zz * zs[i] + ww) / w;
-    }
-
-    return (i: number) => {
-        const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
-        return (xx * x + yy * y + zz * z + ww) / w;
-    }
-}
-
-function makeZ(m: Mat4, xs: ArrayLike<number>, ys: ArrayLike<number>, zs: ArrayLike<number>) {
-    const xx = m[2], yy = m[6], zz = m[10], ww = m[14];
-
-    if (isWConst(m)) {
-        const w = m[15] || 1.0;
-        return (i: number) => (xx * xs[i] + yy * ys[i] + zz * zs[i] + ww) / w;
-    }
-
-    return (i: number) => {
-        const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
-        return (xx * x + yy * y + zz * z + ww) / w;
-    }
-}
-
-function makeGetPosition(m: Mat4, x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number>) {
-    return (a: number, r: Vec3): Vec3 => {
-        r[0] = x[a];
-        r[1] = y[a];
-        r[2] = z[a];
-        Vec3.transformMat4(r, r, m);
-        return r;
-    }
 }
\ No newline at end of file

src/perf-tests/structure.ts

@@ -242,8 +242,8 @@ export namespace PropertyAccess {
     }
 
     export async function run() {
-        //const { structures, models } = await readCIF('./examples/1cbs_full.bcif');
-        const { structures, models, mmcif } = await readCIF('e:/test/quick/3j3q_full.bcif');
+        const { structures, models, mmcif } = await readCIF('./examples/1cbs_full.bcif');
+        //const { structures, models, mmcif } = await readCIF('e:/test/quick/3j3q_full.bcif');
         //const { structures, models, mmcif } = await readCIF('e:/test/quick/1cbs_updated.cif');
 
         console.log(mmcif.pdbx_struct_oper_list.matrix.toArray());
@@ -272,7 +272,7 @@ export namespace PropertyAccess {
         console.log('atom.x', sumProperty(structures[0], Q.props.atom.x));
         console.timeEnd('atom.x');
         console.time('__x')
-        console.log('__x', sumProperty(structures[0], l => l.unit.conformation.__x[l.atom]));
+        console.log('__x', sumProperty(structures[0], l => l.unit.conformation.x[l.atom]));
         console.timeEnd('__x')
 
         //const authSeqId = Atom.property(l => l.unit.hierarchy.residues.auth_seq_id.value(l.unit.residueIndex[l.atom]));