diff --git a/src/number/type.js b/src/number/type.js index d8e6806..0d48348 100644 --- a/src/number/type.js +++ b/src/number/type.js @@ -1,8 +1,10 @@ import {plain} from './helpers.js' -import {BooleanT} from '#boolean/BooleanT.js' -import {Returns} from '#core/Type.js' +import {NumberT} from './NumberT.js' + +import {Returns, ReturnType} from '#core/Type.js' import {match} from '#core/TypePatterns.js' -import {NumberT} from '#number/NumberT.js' +import {BooleanT} from '#boolean/BooleanT.js' +import {Complex} from '#complex/Complex.js' const num = f => Returns(NumberT, f) @@ -10,5 +12,17 @@ export const number = [ plain(a => a), // conversions from Boolean should be consistent with one and zero: match(BooleanT, num(p => p ? NumberT.one : NumberT.zero)), - match([], num(() => 0)) + match([], num(() => 0)), + match(Complex, (math, C) => { + const im = math.im.resolve(C) + const re = math.re.resolve(C) + const compNum = math.number.resolve(ReturnType(re)) + const isZ = math.isZero.resolve(ReturnType(im)) + return num(z => { + if (!isZ(im(z))) { + throw new RangeError(`can't convert Complex ${z} to number`) + } + return compNum(re(z)) + }) + }) ] diff --git a/src/vector/__test__/arithmetic.spec.js b/src/vector/__test__/arithmetic.spec.js index 5dcc3ca..ce888d0 100644 --- a/src/vector/__test__/arithmetic.spec.js +++ b/src/vector/__test__/arithmetic.spec.js @@ -24,6 +24,30 @@ describe('Vector arithmetic functions', () => { assert.deepStrictEqual( add([[1, 2], [4, 2]], [0, -1]), [[1, 1], [4, 1]]) }) + it('(pseudo)inverts matrices', () => { + const inv = math.invert + // inverses + assert.deepStrictEqual(inv([3, 4, 5]), [3/50, 2/25, 1/10]) + assert.deepStrictEqual(inv([[4]]), [[0.25]]) + assert.deepStrictEqual(inv([[5, 2], [-7, -3]]), [[3, 2], [-7, -5]]) + assert(math.equal( + inv([[3, 0, 2], [2, 1, 0], [1, 4, 2]]), + [[1/10, 2/5, -1/10], [-1/5, 1/5, 1/5], [7/20, -3/5, 3/20]])) + // pseudoinverses + assert.deepStrictEqual(inv([[1, 0], [1, 0]]), [[1/2, 1/2], [0, 0]]) + assert.deepStrictEqual( + inv([[1, 0], [0, 1], [0, 1]]), [[1, 0, 0], [0, 1/2, 1/2]]) + assert.deepStrictEqual( + inv([[1, 0, 0, 0, 2], + [0, 0, 3, 0, 0], + [0, 0, 0, 0, 0], + [0, 4, 0, 0, 0]]), + [[1/5, 0, 0, 0], + [ 0, 0, 0, 1/4], + [ 0, 1/3, 0, 0], + [ 0, 0, 0, 0], + [2/5, 0, 0, 0]]) + }) it('multiplies vectors and matrices', () => { const mult = math.multiply const pyth = [3, 4, 5] @@ -38,6 +62,10 @@ describe('Vector arithmetic functions', () => { assert.deepStrictEqual( mult(mat32, [[1, 2], [3, 4]]), [[-8, -10], [-4, -4], [0, 2]]) + assert(math.equal( + mult([[3, 0, 2], [2, 1, 0], [1, 4, 2]], + [[1/10, 2/5, -1/10], [-1/5, 1/5, 1/5], [7/20, -3/5, 3/20]]), + [[1, 0, 0], [0, 1, 0], [0, 0, 1]])) }) it('negates a vector', () => { assert.deepStrictEqual(math.negate([-3, 4, -5]), [3, -4, 5]) diff --git a/src/vector/__test__/type.spec.js b/src/vector/__test__/type.spec.js index 369033a..1208ffc 100644 --- a/src/vector/__test__/type.spec.js +++ b/src/vector/__test__/type.spec.js @@ -22,4 +22,24 @@ describe('Vector type functions', () => { assert.deepStrictEqual( tsp([[1, 2, 3], [4, 5, 6]]), [[1, 4], [2, 5], [3, 6]]) }) + it('can take adjoint (conjugate transpose) of a matrix', () => { + const cx = math.complex + assert.deepStrictEqual( + math.adjoint([[cx(1, 1), cx(2, 2)], [cx(3, 3), cx(4, 4)]]), + [[cx(1, -1), cx(3, -3)], [cx(2, -2), cx(4, -4)]]) + }) + it('generates identity from an example matrix or a number of rows', () => { + const id = math.identity + const cx = math.complex + assert.deepStrictEqual(id(2), [[1, 0], [0, 1]]) + assert.deepStrictEqual(id(cx(2)), [[cx(1), cx(0)], [cx(0), cx(1)]]) + assert.deepStrictEqual( + id([[1, 2, 3]]), + [[1, 0 , 0], [0, 1, 0], [0, 0, 1]]) + }) + it('takes the determinant of a matrix', () => { + assert.strictEqual( + math.determinant([[6, 1, 1], [4, -2, 5], [2, 8, 7]]), + -306) + }) }) diff --git a/src/vector/arithmetic.js b/src/vector/arithmetic.js index 1687138..71332c2 100644 --- a/src/vector/arithmetic.js +++ b/src/vector/arithmetic.js @@ -3,7 +3,7 @@ import { } from './helpers.js' import {Vector} from './Vector.js' -import {ReturnType} from '#core/Type.js' +import {Returns, ReturnType} from '#core/Type.js' import {match} from '#core/TypePatterns.js' import {ReturnsAs} from '#generic/helpers.js' @@ -19,6 +19,201 @@ export const abs = promoteUnary('abs') export const add = promoteBinary('add') export const dotMultiply = promoteBinary('multiply') +export const invert = match(Vector, (math, V, strategy) => { + if (V.vectorDepth > 2) { + throw new TypeError( + 'invert not implemented for arrays of dimension > 2') + } + const normsq = math.normsq.resolve(V, strategy) + const NormT = ReturnType(normsq) + const zNorm = math.zero(NormT) + const isRealZ = math.isZero.resolve(NormT, strategy) + if (V.vectorDepth === 1) { + const invNorm = math.invert.resolve(NormT, strategy) + const scalarMult = math.multiply.resolve( + [V, ReturnType(invNorm)], strategy) + return ReturnsAs(scalarMult, v => { + const nsq = normsq(v) + if (isRealZ(nsq)) return Array(v.length).fill(zNorm) + return scalarMult(v, invNorm(nsq)) + }) + } + // usual matrix situation, want to find a matrix whose product with v + // is the identity, or as close as we can get to that if the rank is + // deficient. We use the Moore-Penrose pseudoinverse. + const clone = math.clone.resolve(V, strategy) + const VComp = V.Component + const Elt = VComp.Component + const invElt = math.invert.resolve(Elt, strategy) + const det = math.determinant.resolve(V, strategy) + const neg = math.negate.resolve(Elt, strategy) + const multMM = math.multiply.resolve([V, V], strategy) + const multMS = math.multiply.resolve([V, ReturnType(invElt)], strategy) + const multVS = math.multiply.resolve([VComp, ReturnType(invElt)], strategy) + const multSS = math.multiply.resolve([Elt, Elt], strategy) + const sub = math.subtract.resolve([Elt, Elt], strategy) + const id = math.identity.resolve(V, strategy) + const abs = math.abs.resolve(Elt, strategy) + const gt = math.larger.resolve([NormT, NormT], strategy) + const isEltZ = math.isZero.resolve(Elt, strategy) + const zElt = math.zero(Elt) + const adj = math.adjoint.resolve(V, strategy) + const methods = { + invElt, det, isRealZ, neg, multMM, multMS, multVS, multSS, + id, abs, isEltZ, zElt, gt, sub, adj, clone + } + if (ReturnType(abs) !== NormT) { + throw new TypeError('type inconsistency in matrix invert') + } + return Returns(V, m => { + const rows = m.length + const cols = m[0].length + const nsq = normsq(m) + if (isRealZ(nsq)) { // all-zero matrix + const retval = [] + for (let ix = 0; ix < cols; ++ix) { + retval.push(Array(rows).fill(zNorm)) + } + return retval + } + if (rows == cols) { + // the inv helper will return falsy if not invertible + const retval = inv(m, rows, methods) + if (retval) return retval + } + + return pinv(m, rows, cols, methods) + }) +}) + +// Returns the inverse of a rows×rows matrix or false if not invertible +// Note: destroys m in the inversion process. +function inv( + origm, + rows, + { + invElt, det, isRealZ, neg, multMS, multVS, multSS, + id, abs, isEltZ, zElt, gt, sub, clone + } +) { + switch (rows) { + case 1: return [[invElt(origm[0][0])]] + case 2: { + const dt = det(origm) + if (isRealZ(dt)) return false + const divisor = invElt(dt) + const [[a, b], [c, d]] = origm + return multMS([[d, neg(b)], [neg(c), a]], divisor) + } + default: { // Gauss-Jordan elimination + const m = clone(origm) + const B = id(m) + const cols = rows + // Loop over columns, performing row reductions: + for (let c = 0; c < cols; ++c) { + // Pivot: Find row r that has the largest entry in column c, and + // swap row c and r: + let colMax = abs(m[c][c]) + let rMax = c + for (let r = c + 1; r < rows; ++r) { + const mag = abs(m[r][c]) + if (gt(mag, colMax)) { + colMax = mag + rMax = r + } + } + if (isRealZ(colMax)) return false + if (rMax !== c) { + [m[c], m[rMax], B[c], B[rMax]] + = [m[rMax], m[c], B[rMax], B[c]] + } + // Normalize the cth row: + const normalizer = invElt(m[c][c]) + const mc = multVS(m[c], normalizer) + m[c] = mc + const Bc = multVS(B[c], normalizer) + B[c] = Bc + + // Eliminate nonzero values on other rows at column c + for (let r = 0; r < rows; ++r) { + if (r === c) continue + const mr = m[r] + const Br = B[r] + const mrc = mr[c] + if (!isEltZ(mr[c])) { + // Subtract Arc times row c from row r to eliminate A[r][c] + mr[c] = zElt + for (let s = c + 1; s < cols; ++s) { + mr[s] = sub(mr[s], multSS(mrc, mc[s])) + } + for (let s = 0; s < cols; ++s) { + Br[s] = sub(Br[s], multSS(mrc, Bc[s])) + } + } + } + } + + return B + }} +} + +// Calculates Moore-Penrose pseudoinverse +// uses rank factorization per mathjs; SVD appears to be considered better +// but not worth the effort to implement for this prototype +function pinv(m, rows, cols, methods) { + const {C, F} = rankFactorization(m, rows, cols, methods) + const {multMM, adj} = methods + const Cstar = adj(C) + const Cpinv = multMM(inv(multMM(Cstar, C), Cstar.length, methods), Cstar) + const Fstar = adj(F) + const Fpinv = multMM(Fstar, inv(multMM(F, Fstar), F.length, methods)) + return multMM(Fpinv, Cpinv) +} + +// warning: destroys m in computing the row-reduced echelon form +// TODO: this code should be merged with inv to the extent possible. It's +// a very similar process. +function rref(origm, rows, cols, methods) { + const {isEltZ, invElt, multVS, zElt, sub, multSS, clone} = methods + const m = clone(origm) + let lead = -1 + for (let r = 0; r < rows && ++lead < cols; ++r) { + if (cols <= lead) return m + let i = r + while (isEltZ(m[i][lead])) { + if (++i === rows) { + i = r + if (++lead === cols) return m + } + } + + if (i !== r) [m[i], m[r]] = [m[r], m[i]] + + let normalizer = invElt(m[r][lead]) + const mr = multVS(m[r], normalizer) + m[r] = mr + for (let i = 0; i < rows; ++i) { + if (i === r) continue + const mi = m[i] + const toRemove = mi[lead] + mi[lead] = zElt + for (let j = lead + 1; j < cols; ++j) { + mi[j] = sub(mi[j], multSS(toRemove, mr[j])) + } + } + } + return m +} + +function rankFactorization(m, rows, cols, methods) { + const RREF = rref(m, rows, cols, methods) + const {isEltZ} = methods + const rankRows = RREF.map(row => row.some(elt => !isEltZ(elt))) + const C = m.map(row => row.filter((_, j) => j < rows && rankRows[j])) + const F = RREF.filter((_, i) => rankRows[i]) + return {C, F} +} + export const multiply = [ distributeFirst('multiply'), distributeSecond('multiply'), diff --git a/src/vector/type.js b/src/vector/type.js index 88cbe3f..ad549d5 100644 --- a/src/vector/type.js +++ b/src/vector/type.js @@ -1,5 +1,5 @@ import {Vector} from './Vector.js' -import {Returns, Unknown} from '#core/Type.js' +import {OneOf, Returns, ReturnType, Unknown} from '#core/Type.js' import {Any, Multiple, match} from '#core/TypePatterns.js' export const vector = match(Multiple(Any), (math, [TV]) => { @@ -9,16 +9,125 @@ export const vector = match(Multiple(Any), (math, [TV]) => { return Returns(Vector(CompType), v => v) }) -export const transpose = match(Vector, (_math, V) => { - const wrapV = V.vectorDepth === 1 - const Mat = wrapV ? Vector(V) : V - return Returns(Mat, v => { - if (wrapV) v = [v] +export const determinant = match(Vector(Vector), (math, M, strategy) => { + const Elt = M.Component.Component + const cloneElt = math.clone.resolve(Elt, strategy) + const mult = math.multiply.resolve([Elt, Elt], strategy) + const sub = math.subtract.resolve( + [ReturnType(mult), ReturnType(mult)], strategy) + const isZ = math.isZero.resolve(Elt, strategy) + const zElt = math.zero(Elt) + const clone = math.clone.resolve(M, strategy) + const div = math.divide.resolve([ReturnType(sub), Elt], strategy) + const neg = math.negate.resolve(Elt, strategy) + return Returns(OneOf(ReturnType(clone), ReturnType(sub), Elt), origm => { + const rows = origm.length + switch (rows) { + case 1: return cloneElt(origm[0][0]) + case 2: { + const [[a, b], [c, d]] = origm + return sub(mult(a, d), mult(b, c)) + } + default: { // Bareiss algorithm + const m = clone(origm) + let negated = false + const rowIndices = [...Array(rows).keys()] // track row indices + // because the algorithm may swap rows + for (let k = 0; k < rows; ++k) { + let k_ = rowIndices[k] + if (isZ(m[k_][k])) { + let _k = k + while (++_k < rows) { + if (!isZ(m[rowIndices[_k]][k])) { + k_ = rowIndices[_k] + rowIndices[_k] = rowIndices[k] + rowIndices[k] = k_ + negated = !negated + break + } + } + if (_k === rows) return zElt + } + const piv = m[k_][k] // we now know nonzero + const piv_ = k === 0 ? 1 : m[rowIndices[k-1]][k-1] + for (let i = k + 1; i < rows; ++i) { + const i_ = rowIndices[i] + for (let j = k + 1; j < rows; ++j) { + m[i_][j] = div( + sub(mult(m[i_][j], piv), mult(m[i_][k], m[k_][j])), + piv_) + } + } + } + const det = m[rowIndices[rows - 1]][rows - 1] + return negated ? neg(det) : det + }} + }) +}) + +function identitizer(cols, zero, one) { + const retval = [] + for (let ix = 0; ix < cols; ++ix) { + const row = Array(cols).fill(zero) + row[ix] = one + retval.push(row) + } + return retval +} + +export const identity = [ + match(Any, (math, V) => { + const toNum = math.number.resolve(V) + const zero = math.zero(V) + const one = math.one(V) + return Returns(Vector(Vector(V)), n => identitizer(toNum(n), zero, one)) + }), + match(Vector, (math, V) => { + switch (V.vectorDepth) { + case 1: { + const Elt = V.Component + const one = math.one(Elt) + return Returns(math.typeOf(one), () => one) + } + case 2: { + const Elt = V.Component.Component + const zero = math.zero(Elt) + const one = math.one(Elt) + return Returns(V, m => identitizer( + m.length ? m[0].length : 0, zero, one)) + } + default: + throw new RangeError( + `'identity' not implemented on ${V.vectorDepth} dimensional arrays`) + } + }) +] + +// transposes a 2D matrix +function transposer(wrap, eltFun) { + return v => { + if (wrap) v = [v] const cols = v.length ? v[0].length : 0 const retval = [] for (let ix = 0; ix < cols; ++ix) { - retval.push(v.map(row => row[ix])) + retval.push(v.map(row => eltFun(row[ix]))) } return retval - }) + } +} + +export const transpose = match(Vector, (_math, V) => { + const wrapV = V.vectorDepth === 1 + const Mat = wrapV ? Vector(V) : V + return Returns(Mat, transposer(wrapV, elt => elt)) +}) + +// or with conjugation: +export const adjoint = match(Vector, (math, V, strategy) => { + const wrapV = V.vectorDepth === 1 + const VComp = V.Component + const Elt = wrapV ? VComp : VComp.Component + const conj = math.conj.resolve(Elt, strategy) + const Mat = Vector(Vector(ReturnType(conj))) + return Returns(Mat, transposer(wrapV, conj)) })