added matrix object, principal axes calculation

src/mol-math/linear-algebra/matrix/matrix.ts

@@ -0,0 +1,75 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+interface Matrix { data: Helpers.NumberArray, size: number, cols: number, rows: number }
+
+namespace Matrix {
+    export function create(cols: number, rows: number, ctor: { new (size: number): Helpers.NumberArray } = Float32Array): Matrix {
+        const size = cols * rows
+        return { data: new ctor(size), size, cols, rows }
+    }
+
+    export function fromArray(data: Helpers.NumberArray, cols: number, rows: number): Matrix {
+        return { data, size: cols * rows, cols, rows }
+    }
+
+    export function transpose(out: Matrix, mat: Matrix): Matrix {
+        const nrows = mat.rows, ncols = mat.cols
+        const md = mat.data, mtd = out.data
+
+        for (let i = 0, mi = 0, mti = 0; i < nrows; mti += 1, mi += ncols, ++i) {
+            let ri = mti
+            for (let j = 0; j < ncols; ri += nrows, j++) mtd[ri] = md[mi + j]
+        }
+        return mat
+    }
+
+    /** out = matA * matB' */
+    export function multiplyABt (out: Matrix, matA: Matrix, matB: Matrix) {
+        const ncols = matA.cols, nrows = matA.rows, mrows = matB.rows
+        const ad = matA.data, bd = matB.data, cd = out.data
+
+        for (let i = 0, matAp = 0, outP = 0; i < nrows; matAp += ncols, i++) {
+            for (let pB = 0, j = 0; j < mrows; outP++, j++) {
+                let sum = 0.0
+                let pMatA = matAp
+                for (let k = 0; k < ncols; pMatA++, pB++, k++) {
+                    sum += ad[pMatA] * bd[pB]
+                }
+                cd[outP] = sum
+            }
+        }
+        return out
+    }
+
+    /** Get the mean of rows in `mat` */
+    export function meanRows (mat: Matrix) {
+        const nrows = mat.rows, ncols = mat.cols
+        const md = mat.data
+        const mean = new Array(ncols)
+
+        for (let j = 0; j < ncols; ++j) mean[ j ] = 0.0
+        for (let i = 0, p = 0; i < nrows; ++i) {
+            for (let j = 0; j < ncols; ++j, ++p) mean[ j ] += md[ p ]
+        }
+        for (let j = 0; j < ncols; ++j) mean[ j ] /= nrows
+
+        return mean
+    }
+
+    /** Subtract `row` from all rows in `mat` */
+    export function subRows (mat: Matrix, row: Helpers.NumberArray) {
+        const nrows = mat.rows, ncols = mat.cols
+        const md = mat.data
+
+        for (let i = 0, p = 0; i < nrows; ++i) {
+            for (let j = 0; j < ncols; ++j, ++p) md[ p ] -= row[ j ]
+        }
+        return mat
+    }
+}
+
+export default Matrix

src/mol-math/linear-algebra/matrix/principal-axes.ts

@@ -0,0 +1,191 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import Matrix from './matrix.js';
+import { Vec3 } from '../3d.js';
+// import { Vec3, Mat4 } from '../3d.js';
+import { svd } from './svd.js';
+
+// const negateVector = Vec3.create(-1, -1, -1)
+// const tmpMatrix = Mat4.identity()
+
+/**
+ * Principal axes
+ */
+class PrincipalAxes {
+    begA: Vec3
+    endA: Vec3
+    begB: Vec3
+    endB: Vec3
+    begC: Vec3
+    endC: Vec3
+
+    center: Vec3
+
+    vecA: Vec3
+    vecB: Vec3
+    vecC: Vec3
+
+    normVecA: Vec3
+    normVecB: Vec3
+    normVecC: Vec3
+
+    /**
+     * points is a 3xN matrix
+     */
+    constructor(points: Matrix) {
+        const n = points.rows
+        const n3 = n / 3
+        const pointsT = Matrix.create(n, 3)
+        const A = Matrix.create(3, 3)
+        const W = Matrix.create(1, 3)
+        const U = Matrix.create(3, 3)
+        const V = Matrix.create(3, 3)
+
+        // calculate
+        const mean = Matrix.meanRows(points)
+        Matrix.subRows(points, mean)
+        Matrix.transpose(pointsT, points)
+        Matrix.multiplyABt(A, pointsT, pointsT)
+        svd(A, W, U, V)
+
+        // center
+        const vm = Vec3.create(mean[0], mean[1], mean[2])
+
+        // normalized
+        const van = Vec3.create(U.data[0], U.data[3], U.data[6])
+        const vbn = Vec3.create(U.data[1], U.data[4], U.data[7])
+        const vcn = Vec3.create(U.data[2], U.data[5], U.data[8])
+
+        // scaled
+        const va = Vec3.scale(Vec3.zero(), van, Math.sqrt(W.data[0] / n3))
+        const vb = Vec3.scale(Vec3.zero(), vbn, Math.sqrt(W.data[1] / n3))
+        const vc = Vec3.scale(Vec3.zero(), vcn, Math.sqrt(W.data[2] / n3))
+        // const va = van.clone().multiplyScalar(Math.sqrt(W.data[0] / n3))
+        // const vb = vbn.clone().multiplyScalar(Math.sqrt(W.data[1] / n3))
+        // const vc = vcn.clone().multiplyScalar(Math.sqrt(W.data[2] / n3))
+
+        // points
+        this.begA = Vec3.sub(Vec3.clone(vm), vm, va)
+        this.endA = Vec3.add(Vec3.clone(vm), vm, va)
+        this.begB = Vec3.sub(Vec3.clone(vm), vm, vb)
+        this.endB = Vec3.add(Vec3.clone(vm), vm, vb)
+        this.begC = Vec3.sub(Vec3.clone(vm), vm, vc)
+        this.endC = Vec3.add(Vec3.clone(vm), vm, vc)
+        // this.begA = vm.clone().sub(va)
+        // this.endA = vm.clone().add(va)
+        // this.begB = vm.clone().sub(vb)
+        // this.endB = vm.clone().add(vb)
+        // this.begC = vm.clone().sub(vc)
+        // this.endC = vm.clone().add(vc)
+
+        //
+
+        this.center = vm
+
+        this.vecA = va
+        this.vecB = vb
+        this.vecC = vc
+
+        this.normVecA = van
+        this.normVecB = vbn
+        this.normVecC = vcn
+    }
+
+    // TODO
+    // /**
+    //  * Get the basis matrix descriping the axes
+    //  * @param  {Matrix4} [optionalTarget] - target object
+    //  * @return {Matrix4} the basis
+    //  */
+    // getBasisMatrix(optionalTarget = new Matrix4()) {
+    //     const basis = optionalTarget
+
+    //     basis.makeBasis(this.normVecB, this.normVecA, this.normVecC)
+    //     if (basis.determinant() < 0) {
+    //         basis.scale(negateVector)
+    //     }
+
+    //     return basis
+    // }
+
+    // TODO
+    // /**
+    //  * Get a quaternion descriping the axes rotation
+    //  * @param  {Quaternion} [optionalTarget] - target object
+    //  * @return {Quaternion} the rotation
+    //  */
+    // getRotationQuaternion(optionalTarget = new Quaternion()) {
+    //     const q = optionalTarget
+    //     q.setFromRotationMatrix(this.getBasisMatrix(tmpMatrix))
+
+    //     return q.inverse()
+    // }
+
+    // TODO
+    // /**
+    //  * Get the scale/length for each dimension for a box around the axes
+    //  * to enclose the atoms of a structure
+    //  * @param  {Structure|StructureView} structure - the structure
+    //  * @return {{d1a: Number, d2a: Number, d3a: Number, d1b: Number, d2b: Number, d3b: Number}} scale
+    //  */
+    // getProjectedScaleForAtoms(structure: Structure) {
+    //     let d1a = -Infinity
+    //     let d1b = -Infinity
+    //     let d2a = -Infinity
+    //     let d2b = -Infinity
+    //     let d3a = -Infinity
+    //     let d3b = -Infinity
+
+    //     const p = new Vector3()
+    //     const t = new Vector3()
+
+    //     const center = this.center
+    //     const ax1 = this.normVecA
+    //     const ax2 = this.normVecB
+    //     const ax3 = this.normVecC
+
+    //     structure.eachAtom(function (ap: AtomProxy) {
+    //         projectPointOnVector(p.copy(ap as any), ax1, center)  // TODO
+    //         const dp1 = t.subVectors(p, center).normalize().dot(ax1)
+    //         const dt1 = p.distanceTo(center)
+    //         if (dp1 > 0) {
+    //             if (dt1 > d1a) d1a = dt1
+    //         } else {
+    //             if (dt1 > d1b) d1b = dt1
+    //         }
+
+    //         projectPointOnVector(p.copy(ap as any), ax2, center)
+    //         const dp2 = t.subVectors(p, center).normalize().dot(ax2)
+    //         const dt2 = p.distanceTo(center)
+    //         if (dp2 > 0) {
+    //             if (dt2 > d2a) d2a = dt2
+    //         } else {
+    //             if (dt2 > d2b) d2b = dt2
+    //         }
+
+    //         projectPointOnVector(p.copy(ap as any), ax3, center)
+    //         const dp3 = t.subVectors(p, center).normalize().dot(ax3)
+    //         const dt3 = p.distanceTo(center)
+    //         if (dp3 > 0) {
+    //             if (dt3 > d3a) d3a = dt3
+    //         } else {
+    //             if (dt3 > d3b) d3b = dt3
+    //         }
+    //     })
+
+    //     return {
+    //         d1a: d1a,
+    //         d2a: d2a,
+    //         d3a: d3a,
+    //         d1b: -d1b,
+    //         d2b: -d2b,
+    //         d3b: -d3b
+    //     }
+    // }
+}
+
+export default PrincipalAxes

src/mol-math/linear-algebra/matrix/svd.ts

@@ -0,0 +1,292 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import Matrix from './matrix';
+
+// svd method adapted from http://inspirit.github.io/jsfeat/ MIT Eugene Zatepyakin
+
+export function swap(A: Helpers.NumberArray, i0: number, i1: number, t: number) {
+    t = A[i0]
+    A[i0] = A[i1]
+    A[i1] = t
+}
+
+export function hypot(a: number, b: number) {
+    a = Math.abs(a)
+    b = Math.abs(b)
+    if (a > b) {
+        b /= a
+        return a * Math.sqrt(1.0 + b * b)
+    }
+    if (b > 0) {
+        a /= b
+        return b * Math.sqrt(1.0 + a * a)
+    }
+    return 0.0
+}
+
+const EPSILON = 0.0000001192092896
+const FLT_MIN = 1E-37
+
+export function JacobiSVDImpl(At: Helpers.NumberArray, astep: number, _W: Helpers.NumberArray, Vt: Helpers.NumberArray, vstep: number, m: number, n: number, n1: number) {
+    const eps = EPSILON * 2.0
+    const minval = FLT_MIN
+    let i = 0
+    let j = 0
+    let k = 0
+    let iter = 0
+    const maxIter = Math.max(m, 30)
+    let Ai = 0
+    let Aj = 0
+    let Vi = 0
+    let Vj = 0
+    let changed = 0
+    let c = 0.0
+    let s = 0.0
+    let t = 0.0
+    let t0 = 0.0
+    let t1 = 0.0
+    let sd = 0.0
+    let beta = 0.0
+    let gamma = 0.0
+    let delta = 0.0
+    let a = 0.0
+    let p = 0.0
+    let b = 0.0
+    let seed = 0x1234
+    let val = 0.0
+    let val0 = 0.0
+    let asum = 0.0
+
+    const W = new Float64Array(n << 3)
+
+    for (; i < n; i++) {
+        for (k = 0, sd = 0; k < m; k++) {
+            t = At[i * astep + k]
+            sd += t * t
+        }
+        W[i] = sd
+
+        if (Vt) {
+            for (k = 0; k < n; k++) {
+                Vt[i * vstep + k] = 0
+            }
+            Vt[i * vstep + i] = 1
+        }
+    }
+
+    for (; iter < maxIter; iter++) {
+        changed = 0
+
+        for (i = 0; i < n - 1; i++) {
+            for (j = i + 1; j < n; j++) {
+                Ai = (i * astep) | 0
+                Aj = (j * astep) | 0
+                a = W[i]
+                p = 0
+                b = W[j]
+
+                k = 2
+                p += At[Ai] * At[Aj]
+                p += At[Ai + 1] * At[Aj + 1]
+
+                for (; k < m; k++) { p += At[Ai + k] * At[Aj + k] }
+
+                if (Math.abs(p) <= eps * Math.sqrt(a * b)) continue
+
+                p *= 2.0
+                beta = a - b
+                gamma = hypot(p, beta)
+                if (beta < 0) {
+                    delta = (gamma - beta) * 0.5
+                    s = Math.sqrt(delta / gamma)
+                    c = (p / (gamma * s * 2.0))
+                } else {
+                    c = Math.sqrt((gamma + beta) / (gamma * 2.0))
+                    s = (p / (gamma * c * 2.0))
+                }
+
+                a = 0.0
+                b = 0.0
+
+                k = 2 // unroll
+                t0 = c * At[Ai] + s * At[Aj]
+                t1 = -s * At[Ai] + c * At[Aj]
+                At[Ai] = t0; At[Aj] = t1
+                a += t0 * t0; b += t1 * t1
+
+                t0 = c * At[Ai + 1] + s * At[Aj + 1]
+                t1 = -s * At[Ai + 1] + c * At[Aj + 1]
+                At[Ai + 1] = t0; At[Aj + 1] = t1
+                a += t0 * t0; b += t1 * t1
+
+                for (; k < m; k++) {
+                    t0 = c * At[Ai + k] + s * At[Aj + k]
+                    t1 = -s * At[Ai + k] + c * At[Aj + k]
+                    At[Ai + k] = t0; At[Aj + k] = t1
+
+                    a += t0 * t0; b += t1 * t1
+                }
+
+                W[i] = a
+                W[j] = b
+
+                changed = 1
+
+                if (Vt) {
+                    Vi = (i * vstep) | 0
+                    Vj = (j * vstep) | 0
+
+                    k = 2
+                    t0 = c * Vt[Vi] + s * Vt[Vj]
+                    t1 = -s * Vt[Vi] + c * Vt[Vj]
+                    Vt[Vi] = t0; Vt[Vj] = t1
+
+                    t0 = c * Vt[Vi + 1] + s * Vt[Vj + 1]
+                    t1 = -s * Vt[Vi + 1] + c * Vt[Vj + 1]
+                    Vt[Vi + 1] = t0; Vt[Vj + 1] = t1
+
+                    for (; k < n; k++) {
+                        t0 = c * Vt[Vi + k] + s * Vt[Vj + k]
+                        t1 = -s * Vt[Vi + k] + c * Vt[Vj + k]
+                        Vt[Vi + k] = t0; Vt[Vj + k] = t1
+                    }
+                }
+            }
+        }
+        if (changed === 0) break
+    }
+
+    for (i = 0; i < n; i++) {
+        for (k = 0, sd = 0; k < m; k++) {
+            t = At[i * astep + k]
+            sd += t * t
+        }
+        W[i] = Math.sqrt(sd)
+    }
+
+    for (i = 0; i < n - 1; i++) {
+        j = i
+        for (k = i + 1; k < n; k++) {
+            if (W[j] < W[k]) { j = k }
+        }
+        if (i !== j) {
+            swap(W, i, j, sd)
+            if (Vt) {
+                for (k = 0; k < m; k++) {
+                    swap(At, i * astep + k, j * astep + k, t)
+                }
+
+                for (k = 0; k < n; k++) {
+                    swap(Vt, i * vstep + k, j * vstep + k, t)
+                }
+            }
+        }
+    }
+
+    for (i = 0; i < n; i++) {
+        _W[i] = W[i]
+    }
+
+    if (!Vt) {
+        return
+    }
+
+    for (i = 0; i < n1; i++) {
+        sd = i < n ? W[i] : 0
+
+        while (sd <= minval) {
+            // if we got a zero singular value, then in order to get the corresponding left singular vector
+            // we generate a random vector, project it to the previously computed left singular vectors,
+            // subtract the projection and normalize the difference.
+            val0 = (1.0 / m)
+            for (k = 0; k < m; k++) {
+                seed = (seed * 214013 + 2531011)
+                val = (((seed >> 16) & 0x7fff) & 256) !== 0 ? val0 : -val0
+                At[i * astep + k] = val
+            }
+            for (iter = 0; iter < 2; iter++) {
+                for (j = 0; j < i; j++) {
+                    sd = 0
+                    for (k = 0; k < m; k++) {
+                        sd += At[i * astep + k] * At[j * astep + k]
+                    }
+                    asum = 0.0
+                    for (k = 0; k < m; k++) {
+                        t = (At[i * astep + k] - sd * At[j * astep + k])
+                        At[i * astep + k] = t
+                        asum += Math.abs(t)
+                    }
+                    asum = asum ? 1.0 / asum : 0
+                    for (k = 0; k < m; k++) {
+                        At[i * astep + k] *= asum
+                    }
+                }
+            }
+            sd = 0
+            for (k = 0; k < m; k++) {
+                t = At[i * astep + k]
+                sd += t * t
+            }
+            sd = Math.sqrt(sd)
+        }
+
+        s = (1.0 / sd)
+        for (k = 0; k < m; k++) {
+            At[i * astep + k] *= s
+        }
+    }
+}
+
+export function svd(A: Matrix, W: Matrix, U: Matrix, V: Matrix) {
+    let at = 0
+    let i = 0
+    const _m = A.rows
+    const _n = A.cols
+    let m = _m
+    let n = _n
+
+    if (m < n) {
+        at = 1
+        i = m
+        m = n
+        n = i
+    }
+
+    const amt = Matrix.create(m, m)
+    const wmt = Matrix.create(1, n)
+    const vmt = Matrix.create(n, n)
+
+    if (at === 0) {
+        Matrix.transpose(amt, A)
+    } else {
+        for (i = 0; i < _n * _m; i++) {
+            amt.data[i] = A.data[i]
+        }
+        for (; i < n * m; i++) {
+            amt.data[i] = 0
+        }
+    }
+
+    JacobiSVDImpl(amt.data, m, wmt.data, vmt.data, n, m, n, m)
+
+    if (W) {
+        for (i = 0; i < n; i++) {
+            W.data[i] = wmt.data[i]
+        }
+        for (; i < _n; i++) {
+            W.data[i] = 0
+        }
+    }
+
+    if (at === 0) {
+        if (U) Matrix.transpose(U, amt)
+        if (V) Matrix.transpose(V, vmt)
+    } else {
+        if (U) Matrix.transpose(U, vmt)
+        if (V) Matrix.transpose(V, amt)
+    }
+}