diff --git a/.github/workflows/typedoc.yml b/.github/workflows/typedoc.yml index ff91151..f7816a7 100644 --- a/.github/workflows/typedoc.yml +++ b/.github/workflows/typedoc.yml @@ -1,32 +1,15 @@ -name: Deploy TypeDoc on GitHub pages +name: TypeDoc on: workflow_dispatch: release: types: [published] -env: - NODE_VERSION: 22.x - ENTRY_FILE: 'src/index.js' - jobs: - deploy: - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v4 - - uses: actions/setup-node@v4 - with: - node-version: ${{ env.NODE_VERSION }} - - name: Install dependencies - run: npm install - - name: Build documentation - uses: zakodium/typedoc-action@v2 - with: - entry: ${{ env.ENTRY_FILE }} - - name: Deploy to GitHub pages - uses: JamesIves/github-pages-deploy-action@releases/v4 - with: - token: ${{ secrets.BOT_TOKEN }} - branch: gh-pages - folder: docs - clean: true + typedoc: + # Documentation: https://github.com/zakodium/workflows#typedoc + uses: zakodium/workflows/.github/workflows/typedoc.yml@typedoc-v1 + with: + entry: 'src/index.js' + secrets: + github-token: ${{ secrets.BOT_TOKEN }} diff --git a/benchmark/benchmarkLinePlot.js b/benchmark/benchmarkLinePlot.js new file mode 100644 index 0000000..175b042 --- /dev/null +++ b/benchmark/benchmarkLinePlot.js @@ -0,0 +1,55 @@ +import { bench, do_not_optimize, lineplot, run } from 'mitata'; +import { xSequentialFillFromStep } from 'ml-spectra-processing'; + +import { SparseMatrix } from '../src/index.js'; + +import { SparseMatrix as SparseMatrixOld } from './class/SparseMatrixOld.js'; +import { randomMatrix } from './utils/randomMatrix.js'; + +const density = 0.02; // Fixed density for this comparison; + +// Prepare matrices once +const sizes = Array.from( + xSequentialFillFromStep({ from: 4, step: 4, size: 13 }), +); +lineplot(() => { + bench('Sparse.mmul($size)', function* (ctx) { + const size = ctx.get('size'); + + // Prepare matrices once + const A = new SparseMatrix(randomMatrix(size, size, density)); + const B = new SparseMatrix(randomMatrix(size, size, density)); + // Benchmark the multiplication + yield () => do_not_optimize(A.mmul(B)); + }) + .gc('inner') + .args('size', sizes); // 16, 32, 64, 128, 256 + + bench('SparseOld.mmul($size)', function* (ctx) { + const size = ctx.get('size'); + const A = randomMatrix(size, size, density); + const B = randomMatrix(size, size, density); + const AOld = new SparseMatrixOld(A); + const BOld = new SparseMatrixOld(B); + + // Benchmark the multiplication + yield () => do_not_optimize(AOld.mmul(BOld)); + }) + .gc('inner') + .args('size', sizes); + + // bench('Dense.mmul($size)', function* (ctx) { + // const size = ctx.get('size'); + + // // Prepare matrices once + // const A = randomMatrix(size, size, density); + // const B = randomMatrix(size, size, density); + // const ADense = new Matrix(A); + // const BDense = new Matrix(B); + + // // Benchmark the multiplication + // yield () => do_not_optimize(ADense.mmul(BDense)); + // }).gc('inner').range('size', min, max, multiplier); +}); + +await run(); diff --git a/benchmark/class/SparseMatrixOld.js b/benchmark/class/SparseMatrixOld.js new file mode 100644 index 0000000..fca3b8e --- /dev/null +++ b/benchmark/class/SparseMatrixOld.js @@ -0,0 +1,1385 @@ +import HashTable from 'ml-hash-table'; + +/** @typedef {(row: number, column: number, value: number) => number | false} ForEachNonZeroCallback */ + +export class SparseMatrix { + constructor(rows, columns, options = {}) { + if (rows instanceof SparseMatrix) { + // clone + const other = rows; + this._init( + other.rows, + other.columns, + other.elements.clone(), + other.threshold, + ); + return; + } + + if (Array.isArray(rows)) { + const matrix = rows; + rows = matrix.length; + options = columns || {}; + columns = matrix[0].length; + this._init(rows, columns, new HashTable(options), options.threshold); + for (let i = 0; i < rows; i++) { + for (let j = 0; j < columns; j++) { + let value = matrix[i][j]; + if (this.threshold && Math.abs(value) < this.threshold) value = 0; + if (value !== 0) { + this.elements.set(i * columns + j, matrix[i][j]); + } + } + } + } else { + this._init(rows, columns, new HashTable(options), options.threshold); + } + } + + _init(rows, columns, elements, threshold) { + this.rows = rows; + this.columns = columns; + this.elements = elements; + this.threshold = threshold || 0; + } + + static eye(rows = 1, columns = rows) { + const min = Math.min(rows, columns); + const matrix = new SparseMatrix(rows, columns, { initialCapacity: min }); + for (let i = 0; i < min; i++) { + matrix.set(i, i, 1); + } + return matrix; + } + + clone() { + return new SparseMatrix(this); + } + + /** + * @returns {number[][]} + */ + to2DArray() { + const copy = new Array(this.rows); + for (let i = 0; i < this.rows; i++) { + copy[i] = new Array(this.columns); + for (let j = 0; j < this.columns; j++) { + copy[i][j] = this.get(i, j); + } + } + return copy; + } + + isSquare() { + return this.rows === this.columns; + } + + isSymmetric() { + if (!this.isSquare()) return false; + + let symmetric = true; + this.forEachNonZero((i, j, v) => { + if (this.get(j, i) !== v) { + symmetric = false; + return false; + } + return v; + }); + return symmetric; + } + + /** + * Search for the wither band in the main diagonals + * @returns {number} + */ + bandWidth() { + let min = this.columns; + let max = -1; + this.forEachNonZero((i, j, v) => { + let diff = i - j; + min = Math.min(min, diff); + max = Math.max(max, diff); + return v; + }); + return max - min; + } + + /** + * Test if a matrix is consider banded using a threshold + * @param {number} width + * @returns {boolean} + */ + isBanded(width) { + let bandWidth = this.bandWidth(); + return bandWidth <= width; + } + + /** + * @returns {number} + */ + get cardinality() { + return this.elements.size; + } + + /** + * @returns {number} + */ + get size() { + return this.rows * this.columns; + } + + /** + * @param {number} row + * @param {number} column + * @returns {number} + */ + get(row, column) { + return this.elements.get(row * this.columns + column); + } + + /** + * @param {number} row + * @param {number} column + * @param {number} value + * @returns {this} + */ + set(row, column, value) { + if (this.threshold && Math.abs(value) < this.threshold) value = 0; + if (value === 0) { + this.elements.remove(row * this.columns + column); + } else { + this.elements.set(row * this.columns + column, value); + } + return this; + } + + /** + * @param {SparseMatrix} other + * @returns {SparseMatrix} + */ + mmul(other) { + if (this.columns !== other.rows) { + console.warn( + 'Number of columns of left matrix are not equal to number of rows of right matrix.', + ); + } + + const m = this.rows; + const p = other.columns; + + const result = new SparseMatrix(m, p); + this.forEachNonZero((i, j, v1) => { + other.forEachNonZero((k, l, v2) => { + if (j === k) { + result.set(i, l, result.get(i, l) + v1 * v2); + } + return v2; + }); + return v1; + }); + return result; + } + + /** + * @param {SparseMatrix} other + * @returns {SparseMatrix} + */ + kroneckerProduct(other) { + const m = this.rows; + const n = this.columns; + const p = other.rows; + const q = other.columns; + + const result = new SparseMatrix(m * p, n * q, { + initialCapacity: this.cardinality * other.cardinality, + }); + this.forEachNonZero((i, j, v1) => { + other.forEachNonZero((k, l, v2) => { + result.set(p * i + k, q * j + l, v1 * v2); + return v2; + }); + return v1; + }); + return result; + } + + /** + * Calls `callback` for each value in the matrix that is not zero. + * The callback can return: + * - `false` to stop the iteration immediately + * - the same value to keep it unchanged + * - 0 to remove the value + * - something else to change the value + * @param {ForEachNonZeroCallback} callback + * @returns {this} + */ + forEachNonZero(callback) { + this.elements.forEachPair((key, value) => { + const i = (key / this.columns) | 0; + const j = key % this.columns; + let r = callback(i, j, value); + if (r === false) return false; // stop iteration + if (this.threshold && Math.abs(r) < this.threshold) r = 0; + if (r !== value) { + if (r === 0) { + this.elements.remove(key, true); + } else { + this.elements.set(key, r); + } + } + return true; + }); + this.elements.maybeShrinkCapacity(); + return this; + } + + getNonZeros() { + const cardinality = this.cardinality; + /** @type {number[]} */ + const rows = new Array(cardinality); + /** @type {number[]} */ + const columns = new Array(cardinality); + /** @type {number[]} */ + const values = new Array(cardinality); + let idx = 0; + this.forEachNonZero((i, j, value) => { + rows[idx] = i; + columns[idx] = j; + values[idx] = value; + idx++; + return value; + }); + return { rows, columns, values }; + } + + /** + * @param {number} newThreshold + * @returns {this} + */ + setThreshold(newThreshold) { + if (newThreshold !== 0 && newThreshold !== this.threshold) { + this.threshold = newThreshold; + this.forEachNonZero((i, j, v) => v); + } + return this; + } + + /** + * @returns {SparseMatrix} - New transposed sparse matrix + */ + transpose() { + let trans = new SparseMatrix(this.columns, this.rows, { + initialCapacity: this.cardinality, + }); + this.forEachNonZero((i, j, value) => { + trans.set(j, i, value); + return value; + }); + return trans; + } + + /** + * @returns {boolean} + */ + isEmpty() { + return this.rows === 0 || this.columns === 0; + } + + // The following was generated by `makeMathMethods.mjs` + //#region Math methods. + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + add(value) { + if (typeof value === 'number') return this.addS(value); + return this.addM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + addS(value) { + this.forEachNonZero((i, j, v) => v + value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + addM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) + v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static add(matrix, value) { + return new SparseMatrix(matrix).addS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + sub(value) { + if (typeof value === 'number') return this.subS(value); + return this.subM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + subS(value) { + this.forEachNonZero((i, j, v) => v - value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + subM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) - v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static sub(matrix, value) { + return new SparseMatrix(matrix).subS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + subtract(value) { + if (typeof value === 'number') return this.subtractS(value); + return this.subtractM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + subtractS(value) { + this.forEachNonZero((i, j, v) => v - value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + subtractM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) - v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static subtract(matrix, value) { + return new SparseMatrix(matrix).subtractS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + mul(value) { + if (typeof value === 'number') return this.mulS(value); + return this.mulM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + mulS(value) { + this.forEachNonZero((i, j, v) => v * value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + mulM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) * v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static mul(matrix, value) { + return new SparseMatrix(matrix).mulS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + multiply(value) { + if (typeof value === 'number') return this.multiplyS(value); + return this.multiplyM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + multiplyS(value) { + this.forEachNonZero((i, j, v) => v * value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + multiplyM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) * v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static multiply(matrix, value) { + return new SparseMatrix(matrix).multiplyS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + div(value) { + if (typeof value === 'number') return this.divS(value); + return this.divM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + divS(value) { + this.forEachNonZero((i, j, v) => v / value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + divM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) / v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static div(matrix, value) { + return new SparseMatrix(matrix).divS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + divide(value) { + if (typeof value === 'number') return this.divideS(value); + return this.divideM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + divideS(value) { + this.forEachNonZero((i, j, v) => v / value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + divideM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) / v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static divide(matrix, value) { + return new SparseMatrix(matrix).divideS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + mod(value) { + if (typeof value === 'number') return this.modS(value); + return this.modM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + modS(value) { + this.forEachNonZero((i, j, v) => v % value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + modM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) % v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static mod(matrix, value) { + return new SparseMatrix(matrix).modS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + modulus(value) { + if (typeof value === 'number') return this.modulusS(value); + return this.modulusM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + modulusS(value) { + this.forEachNonZero((i, j, v) => v % value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + modulusM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) % v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static modulus(matrix, value) { + return new SparseMatrix(matrix).modulusS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + and(value) { + if (typeof value === 'number') return this.andS(value); + return this.andM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + andS(value) { + this.forEachNonZero((i, j, v) => v & value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + andM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) & v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static and(matrix, value) { + return new SparseMatrix(matrix).andS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + or(value) { + if (typeof value === 'number') return this.orS(value); + return this.orM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + orS(value) { + this.forEachNonZero((i, j, v) => v | value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + orM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) | v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static or(matrix, value) { + return new SparseMatrix(matrix).orS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + xor(value) { + if (typeof value === 'number') return this.xorS(value); + return this.xorM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + xorS(value) { + this.forEachNonZero((i, j, v) => v ^ value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + xorM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) ^ v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static xor(matrix, value) { + return new SparseMatrix(matrix).xorS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + leftShift(value) { + if (typeof value === 'number') return this.leftShiftS(value); + return this.leftShiftM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + leftShiftS(value) { + this.forEachNonZero((i, j, v) => v << value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + leftShiftM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) << v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static leftShift(matrix, value) { + return new SparseMatrix(matrix).leftShiftS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + signPropagatingRightShift(value) { + if (typeof value === 'number') { + return this.signPropagatingRightShiftS(value); + } + return this.signPropagatingRightShiftM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + signPropagatingRightShiftS(value) { + this.forEachNonZero((i, j, v) => v >> value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + signPropagatingRightShiftM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) >> v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static signPropagatingRightShift(matrix, value) { + return new SparseMatrix(matrix).signPropagatingRightShiftS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + rightShift(value) { + if (typeof value === 'number') return this.rightShiftS(value); + return this.rightShiftM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + rightShiftS(value) { + this.forEachNonZero((i, j, v) => v >>> value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + rightShiftM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) >>> v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static rightShift(matrix, value) { + return new SparseMatrix(matrix).rightShiftS(value); + } + + /** + * @param {number|SparseMatrix} value + * @returns {this} + */ + zeroFillRightShift(value) { + if (typeof value === 'number') return this.zeroFillRightShiftS(value); + return this.zeroFillRightShiftM(value); + } + + /** + * @param {number} value + * @returns {this} + */ + zeroFillRightShiftS(value) { + this.forEachNonZero((i, j, v) => v >>> value); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {this} + */ + zeroFillRightShiftM(matrix) { + matrix.forEachNonZero((i, j, v) => { + this.set(i, j, this.get(i, j) >>> v); + return v; + }); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @param {number} value + * @returns {SparseMatrix} + */ + static zeroFillRightShift(matrix, value) { + return new SparseMatrix(matrix).zeroFillRightShiftS(value); + } + + /** + * @returns {this} + */ + not() { + this.forEachNonZero((i, j, v) => ~v); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static not(matrix) { + return new SparseMatrix(matrix).not(); + } + + /** + * @returns {this} + */ + abs() { + this.forEachNonZero((i, j, v) => Math.abs(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static abs(matrix) { + return new SparseMatrix(matrix).abs(); + } + + /** + * @returns {this} + */ + acos() { + this.forEachNonZero((i, j, v) => Math.acos(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static acos(matrix) { + return new SparseMatrix(matrix).acos(); + } + + /** + * @returns {this} + */ + acosh() { + this.forEachNonZero((i, j, v) => Math.acosh(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static acosh(matrix) { + return new SparseMatrix(matrix).acosh(); + } + + /** + * @returns {this} + */ + asin() { + this.forEachNonZero((i, j, v) => Math.asin(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static asin(matrix) { + return new SparseMatrix(matrix).asin(); + } + + /** + * @returns {this} + */ + asinh() { + this.forEachNonZero((i, j, v) => Math.asinh(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static asinh(matrix) { + return new SparseMatrix(matrix).asinh(); + } + + /** + * @returns {this} + */ + atan() { + this.forEachNonZero((i, j, v) => Math.atan(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static atan(matrix) { + return new SparseMatrix(matrix).atan(); + } + + /** + * @returns {this} + */ + atanh() { + this.forEachNonZero((i, j, v) => Math.atanh(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static atanh(matrix) { + return new SparseMatrix(matrix).atanh(); + } + + /** + * @returns {this} + */ + cbrt() { + this.forEachNonZero((i, j, v) => Math.cbrt(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static cbrt(matrix) { + return new SparseMatrix(matrix).cbrt(); + } + + /** + * @returns {this} + */ + ceil() { + this.forEachNonZero((i, j, v) => Math.ceil(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static ceil(matrix) { + return new SparseMatrix(matrix).ceil(); + } + + /** + * @returns {this} + */ + clz32() { + this.forEachNonZero((i, j, v) => Math.clz32(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static clz32(matrix) { + return new SparseMatrix(matrix).clz32(); + } + + /** + * @returns {this} + */ + cos() { + this.forEachNonZero((i, j, v) => Math.cos(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static cos(matrix) { + return new SparseMatrix(matrix).cos(); + } + + /** + * @returns {this} + */ + cosh() { + this.forEachNonZero((i, j, v) => Math.cosh(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static cosh(matrix) { + return new SparseMatrix(matrix).cosh(); + } + + /** + * @returns {this} + */ + exp() { + this.forEachNonZero((i, j, v) => Math.exp(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static exp(matrix) { + return new SparseMatrix(matrix).exp(); + } + + /** + * @returns {this} + */ + expm1() { + this.forEachNonZero((i, j, v) => Math.expm1(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static expm1(matrix) { + return new SparseMatrix(matrix).expm1(); + } + + /** + * @returns {this} + */ + floor() { + this.forEachNonZero((i, j, v) => Math.floor(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static floor(matrix) { + return new SparseMatrix(matrix).floor(); + } + + /** + * @returns {this} + */ + fround() { + this.forEachNonZero((i, j, v) => Math.fround(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static fround(matrix) { + return new SparseMatrix(matrix).fround(); + } + + /** + * @returns {this} + */ + log() { + this.forEachNonZero((i, j, v) => Math.log(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static log(matrix) { + return new SparseMatrix(matrix).log(); + } + + /** + * @returns {this} + */ + log1p() { + this.forEachNonZero((i, j, v) => Math.log1p(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static log1p(matrix) { + return new SparseMatrix(matrix).log1p(); + } + + /** + * @returns {this} + */ + log10() { + this.forEachNonZero((i, j, v) => Math.log10(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static log10(matrix) { + return new SparseMatrix(matrix).log10(); + } + + /** + * @returns {this} + */ + log2() { + this.forEachNonZero((i, j, v) => Math.log2(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static log2(matrix) { + return new SparseMatrix(matrix).log2(); + } + + /** + * @returns {this} + */ + round() { + this.forEachNonZero((i, j, v) => Math.round(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static round(matrix) { + return new SparseMatrix(matrix).round(); + } + + /** + * @returns {this} + */ + sign() { + this.forEachNonZero((i, j, v) => Math.sign(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static sign(matrix) { + return new SparseMatrix(matrix).sign(); + } + + /** + * @returns {this} + */ + sin() { + this.forEachNonZero((i, j, v) => Math.sin(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static sin(matrix) { + return new SparseMatrix(matrix).sin(); + } + + /** + * @returns {this} + */ + sinh() { + this.forEachNonZero((i, j, v) => Math.sinh(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static sinh(matrix) { + return new SparseMatrix(matrix).sinh(); + } + + /** + * @returns {this} + */ + sqrt() { + this.forEachNonZero((i, j, v) => Math.sqrt(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static sqrt(matrix) { + return new SparseMatrix(matrix).sqrt(); + } + + /** + * @returns {this} + */ + tan() { + this.forEachNonZero((i, j, v) => Math.tan(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static tan(matrix) { + return new SparseMatrix(matrix).tan(); + } + + /** + * @returns {this} + */ + tanh() { + this.forEachNonZero((i, j, v) => Math.tanh(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static tanh(matrix) { + return new SparseMatrix(matrix).tanh(); + } + + /** + * @returns {this} + */ + trunc() { + this.forEachNonZero((i, j, v) => Math.trunc(v)); + return this; + } + + /** + * @param {SparseMatrix} matrix + * @returns {SparseMatrix} + */ + static trunc(matrix) { + return new SparseMatrix(matrix).trunc(); + } + //#endregion +} + +SparseMatrix.prototype.klass = 'Matrix'; + +SparseMatrix.identity = SparseMatrix.eye; +SparseMatrix.prototype.tensorProduct = SparseMatrix.prototype.kroneckerProduct; diff --git a/benchmark/combinationOfCardinalityAndSize.js b/benchmark/combinationOfCardinalityAndSize.js new file mode 100644 index 0000000..e2b3907 --- /dev/null +++ b/benchmark/combinationOfCardinalityAndSize.js @@ -0,0 +1,164 @@ +import { writeFile } from 'node:fs/promises'; +import path from 'node:path'; + +import { run, bench, do_not_optimize, lineplot } from 'mitata'; +import { Matrix } from 'ml-matrix'; +import { xSequentialFillFromStep } from 'ml-spectra-processing'; + +import { SparseMatrix } from '../src/index.js'; +import { mmulLowDensity } from '../src/utils/mmulLowDensity.js'; +import { mmulMediumDensity } from '../src/utils/mmulMediumDensity.js'; +import { mmulSmall } from '../src/utils/mmulSmall.js'; + +import { randomMatrix } from './utils/randomMatrix.js'; + +/* eslint +func-names: 0 +camelcase: 0 +*/ +// Prepare matrices once +const cardinalities = Array.from( + xSequentialFillFromStep({ from: 10, step: 50, size: 9 }), +); + +// const dimensions = Array.from( +// xSequentialFillFromStep({ from: 30, step: 100, size: 20 }), +// ); + +const dimensions = [32]; +lineplot(() => { + bench('hibrid($cardinality,$dimension)', function* (ctx) { + const cardinality = ctx.get('cardinality'); + const size = ctx.get('dimension'); + // Prepare matrices once + let A = new SparseMatrix(randomMatrix(size, size, cardinality)); + let B = new SparseMatrix(randomMatrix(size, size, cardinality)); + A.mmul(B); + // Benchmark the multiplication + yield () => do_not_optimize(A.mmul(B)); + do_not_optimize(A); + do_not_optimize(B); + }) + .gc('inner') + .args('cardinality', cardinalities) //.range('size', 32, 1024, 2); //.args('size', sizes); + .args('dimension', dimensions); + + bench('small($cardinality,$dimension)', function* (ctx) { + const cardinality = ctx.get('cardinality'); + const size = ctx.get('dimension'); + // Prepare matrices once + let A = new SparseMatrix(randomMatrix(size, size, cardinality)); + let B = new SparseMatrix(randomMatrix(size, size, cardinality)); + mmulSmall(A, B); + // Benchmark the multiplication + yield () => do_not_optimize(mmulSmall(A, B)); + // Explicit cleanup + do_not_optimize(A); + do_not_optimize(B); + }) + .gc('inner') + .args('cardinality', cardinalities) //.range('size', 32, 1024, 2); //.args('size', sizes); + .args('dimension', dimensions); + + bench('low($cardinality,$dimension)', function* (ctx) { + const cardinality = ctx.get('cardinality'); + const size = ctx.get('dimension'); + // Prepare matrices once + let A = new SparseMatrix(randomMatrix(size, size, cardinality)); + let B = new SparseMatrix(randomMatrix(size, size, cardinality)); + mmulLowDensity(A, B); + // Benchmark the multiplication + yield () => do_not_optimize(mmulLowDensity(A, B)); + // Explicit cleanup + do_not_optimize(A); + do_not_optimize(B); + }) + .gc('inner') + .args('cardinality', cardinalities) //.range('size', 32, 1024, 2); //.args('size', sizes); + .args('dimension', dimensions); + + bench('dense($cardinality,$dimension)', function* (ctx) { + const cardinality = ctx.get('cardinality'); + const size = ctx.get('dimension'); + // Prepare matrices once + let A = new Matrix(randomMatrix(size, size, cardinality)); + let B = new Matrix(randomMatrix(size, size, cardinality)); + A.mmul(B); + // Benchmark the multiplication + yield () => do_not_optimize(A.mmul(B)); + // Explicit cleanup + do_not_optimize(A); + do_not_optimize(B); + }) + .gc('inner') + .args('cardinality', cardinalities) //.range('size', 32, 1024, 2); //.args('size', sizes); + .args('dimension', dimensions); + + bench('medium($cardinality,$dimension)', function* (ctx) { + const cardinality = ctx.get('cardinality'); + const size = ctx.get('dimension'); + // Prepare matrices once + let A = new SparseMatrix(randomMatrix(size, size, cardinality)); + let B = new SparseMatrix(randomMatrix(size, size, cardinality)); + mmulMediumDensity(A, B); + // Benchmark the multiplication + yield () => { + do_not_optimize(mmulMediumDensity(A, B)); + }; + + // Explicit cleanup + do_not_optimize(A); + do_not_optimize(B); + }) + .gc('inner') + .args('cardinality', cardinalities) //.range('size', 32, 1024, 2); //.args('size', sizes); + .args('dimension', dimensions); +}); + +// Run benchmarks and capture results +const results = await run({ + // Force GC between every benchmark + gc: true, + // More samples for statistical significance + min_samples: 20, + max_samples: 200, + // Longer warmup to stabilize CPU state + warmup_samples: 10, + warmup_threshold: 100, // ms + // Longer minimum time for stable measurements + min_cpu_time: 2000, // 2 seconds minimum + // Batch settings to reduce variance + batch_samples: 5, + batch_threshold: 10, // ms + // Enable colors + colors: true, +}); + +// Process and store results +const processedResults = []; + +for (const benchmark of results.benchmarks) { + for (const run of benchmark.runs) { + if (run.stats) { + processedResults.push({ + name: benchmark.alias, + cardinality: run.args.cardinality, + dimension: run.args.dimension, + avg: run.stats.avg, + min: run.stats.min, + max: run.stats.max, + p50: run.stats.p50, + p75: run.stats.p75, + p99: run.stats.p99, + samples: run.stats.samples.length, + ticks: run.stats.ticks, + }); + } + } +} + +// Save results to JSON file +await writeFile( + path.join(import.meta.dirname, `benchmark-results-${dimensions[0]}.json`), + JSON.stringify(processedResults, null, 2), +); diff --git a/benchmark/denseMatrix.js b/benchmark/denseMatrix.js deleted file mode 100644 index afab679..0000000 --- a/benchmark/denseMatrix.js +++ /dev/null @@ -1,18 +0,0 @@ -import { Matrix } from 'ml-matrix'; - -let size = 200; - -console.time('dense'); -const matrix = new Matrix(size, size).fill(1); -const matrix2 = new Matrix(size, size).fill(1); - -const product = matrix.mmul(matrix2); -// sum of all the elements -let sum = 0; -for (let i = 0; i < size; i++) { - for (let j = 0; j < size; j++) { - sum += product.get(i, j); - } -} -console.log(sum); -console.timeEnd('dense'); diff --git a/benchmark/densityDecreaseWhenSizeIncrease.js b/benchmark/densityDecreaseWhenSizeIncrease.js new file mode 100644 index 0000000..8ed6101 --- /dev/null +++ b/benchmark/densityDecreaseWhenSizeIncrease.js @@ -0,0 +1,33 @@ +import { bench, do_not_optimize, group, run } from 'mitata'; + +// import { Matrix } from 'ml-matrix'; +import { SparseMatrix } from '../src/index.js'; + +import { SparseMatrix as SparseMatrixOld } from './class/SparseMatrixOld.js'; +import { randomMatrix } from './utils/randomMatrix.js'; + +const sizes = [8, 16, 32, 256, 512, 1024]; +const densities = [0.125, 0.0625, 0.03125, 0.0039, 0.00197, 0.001]; + +for (let i = 0; i < sizes.length; i++) { + const size = sizes[i]; + const density = densities[i]; + const denseA = randomMatrix(size, size, density); + const denseB = randomMatrix(size, size, density); + const AOld = new SparseMatrixOld(denseA); + const BOld = new SparseMatrixOld(denseB); + + const A = new SparseMatrix(denseA); + const B = new SparseMatrix(denseB); + + group(`size:${size}-density:${density}`, () => { + bench('mmulNew', () => { + do_not_optimize(A.mmul(B)); + }).gc('inner'); + bench('mmul', () => { + do_not_optimize(AOld.mmul(BOld)); + }).gc('inner'); + }); +} + +await run({ silent: false }); diff --git a/benchmark/plot.html b/benchmark/plot.html new file mode 100644 index 0000000..4710a4d --- /dev/null +++ b/benchmark/plot.html @@ -0,0 +1,206 @@ + + + + + Benchmark Results Plot + + + + +

Benchmark Results: avg vs cardinality

+ +
+
+ + + diff --git a/benchmark/sparseMatrix.byhand.js b/benchmark/sparseMatrix.byhand.js deleted file mode 100644 index 8261900..0000000 --- a/benchmark/sparseMatrix.byhand.js +++ /dev/null @@ -1,41 +0,0 @@ -import { SparseMatrix } from '../lib/index.js'; - -let size = 200; - -console.time('dense'); - -const matrix1 = new SparseMatrix(size, size); -fillSparseMatrix(matrix1); - -const matrix2 = new SparseMatrix(size, size); -fillSparseMatrix(matrix2); - -const product = new SparseMatrix(size, size); -for (let i = 0; i < size; i++) { - for (let j = 0; j < size; j++) { - for (let k = 0; k < size; k++) { - product.set( - i, - j, - product.get(i, j) + matrix1.get(i, k) * matrix2.get(k, j), - ); - } - } -} - -let sum = 0; -for (let i = 0; i < size; i++) { - for (let j = 0; j < size; j++) { - sum += product.get(i, j); - } -} -console.timeEnd('dense'); -console.log(sum); - -function fillSparseMatrix(matrix) { - for (let row = 0; row < matrix.rows; row++) { - for (let column = 0; column < matrix.columns; column++) { - matrix.set(row, column, 1); - } - } -} diff --git a/benchmark/sparseMatrix.js b/benchmark/sparseMatrix.js deleted file mode 100644 index 53b117e..0000000 --- a/benchmark/sparseMatrix.js +++ /dev/null @@ -1,30 +0,0 @@ -import { SparseMatrix } from '../lib/index.js'; - -let size = 200; - -console.time('dense'); - -const matrix = new SparseMatrix(size, size); -fillSparseMatrix(matrix); - -const matrix2 = new SparseMatrix(size, size); -fillSparseMatrix(matrix2); - -const product = matrix.mmul(matrix2); -// sum of all the elements -let sum = 0; -for (let i = 0; i < size; i++) { - for (let j = 0; j < size; j++) { - sum += product.get(i, j); - } -} -console.timeEnd('dense'); -console.log(sum); - -function fillSparseMatrix(matrix) { - for (let row = 0; row < matrix.rows; row++) { - for (let column = 0; column < matrix.columns; column++) { - matrix.set(row, column, 1); - } - } -} diff --git a/benchmark/utils/randomMatrix.js b/benchmark/utils/randomMatrix.js new file mode 100644 index 0000000..41f4923 --- /dev/null +++ b/benchmark/utils/randomMatrix.js @@ -0,0 +1,21 @@ +export function randomMatrix(rows, cols, cardinality) { + const total = rows * cols; + const positions = new Set(); + + // Generate unique random positions + while (positions.size < cardinality) { + positions.add(Math.floor(Math.random() * total)); + } + + // Build the matrix with zeros + const matrix = Array.from({ length: rows }, () => new Float64Array(cols)); + + // Assign random values to the selected positions + for (const pos of positions) { + const i = Math.floor(pos / cols); + const j = pos % cols; + matrix[i][j] = Math.random() * 10; + } + + return matrix; +} diff --git a/benchmark/utils/randomSparseMatrix.js b/benchmark/utils/randomSparseMatrix.js new file mode 100644 index 0000000..61fe737 --- /dev/null +++ b/benchmark/utils/randomSparseMatrix.js @@ -0,0 +1,24 @@ +import { SparseMatrix } from '../../src/index.js'; + +export function randomSparseMatrix(rows, cols, density = 0.01) { + const total = rows * cols; + const cardinality = Math.round(total * density); + const positions = new Set(); + + // Generate unique random positions + while (positions.size < cardinality) { + positions.add(Math.floor(Math.random() * total)); + } + + // Build the matrix with zeros + const matrix = Array.from({ length: rows }, () => new Float64Array(cols)); + + // Assign random values to the selected positions + for (const pos of positions) { + const i = Math.floor(pos / cols); + const j = pos % cols; + matrix[i][j] = Math.random() * 10; + } + + return new SparseMatrix(matrix); +} diff --git a/eslint.config.mjs b/eslint.config.js similarity index 55% rename from eslint.config.mjs rename to eslint.config.js index 428560c..3ca4395 100644 --- a/eslint.config.mjs +++ b/eslint.config.js @@ -1,9 +1,12 @@ import { defineConfig, globalIgnores } from 'eslint/config'; import ts from 'eslint-config-cheminfo-typescript/base'; -export default defineConfig(globalIgnores(['lib']), ts, { +export default defineConfig(globalIgnores(['benchmark/class/*', 'lib']), ts, { files: ['benchmark/**'], rules: { + camelcase: 'off', + 'func-names': 'off', 'no-console': 'off', + 'no-unused-vars': 'off', }, }); diff --git a/package.json b/package.json index df9f8e2..0801258 100644 --- a/package.json +++ b/package.json @@ -40,7 +40,9 @@ "@vitest/coverage-v8": "^3.1.4", "eslint": "^9.27.0", "eslint-config-cheminfo-typescript": "^18.0.1", + "mitata": "^1.0.34", "ml-matrix": "^6.12.1", + "ml-spectra-processing": "^14.12.0", "prettier": "^3.5.3", "rimraf": "^6.0.1", "typescript": "^5.8.3", diff --git a/src/__tests__/getNonZeros.test.js b/src/__tests__/getNonZeros.test.js new file mode 100644 index 0000000..0520c30 --- /dev/null +++ b/src/__tests__/getNonZeros.test.js @@ -0,0 +1,32 @@ +import { describe, expect, it } from 'vitest'; + +import { SparseMatrix } from '../index.js'; + +describe('Sparse Matrix', () => { + it('getNonZeros', () => { + let m2 = new SparseMatrix( + [ + [0, 0, 0, 0, 0, 0], + [1, 0, 0, 2, 1, 1], + [0, 3, 0, 0, 5, 5], + [0, 0, 0, 0, 0, 0], + [1, 0, 0, 1, 9, 9], + ], + { initialCapacity: 12 }, + ); + + // Default (coordinate list) + expect(m2.getNonZeros()).toEqual({ + rows: Float64Array.from([4, 1, 4, 1, 1, 1, 4, 2, 4, 2, 2]), + columns: Float64Array.from([0, 0, 3, 3, 4, 5, 5, 1, 4, 4, 5]), + values: Float64Array.from([1, 1, 1, 2, 1, 1, 9, 3, 9, 5, 5]), + }); + + // CSR format + expect(m2.getNonZeros({ csr: true })).toEqual({ + rows: Float64Array.from([0, 0, 4, 7, 7, 11]), + columns: Float64Array.from([0, 3, 4, 5, 1, 4, 5, 0, 3, 4, 5]), + values: Float64Array.from([1, 2, 1, 1, 3, 5, 5, 1, 1, 9, 9]), + }); + }); +}); diff --git a/src/__tests__/index.test.js b/src/__tests__/index.test.js index 0843eba..81ad9b3 100644 --- a/src/__tests__/index.test.js +++ b/src/__tests__/index.test.js @@ -1,8 +1,27 @@ +import { Matrix } from 'ml-matrix'; import { describe, expect, it } from 'vitest'; import { SparseMatrix } from '../index.js'; describe('Sparse Matrix', () => { + it('add', () => { + let m1 = new SparseMatrix([ + [2, 0, 1], + [0, 0, 3], + [2, 0, 1], + ]); + let m2 = new SparseMatrix([ + [0, 1, 5], + [2, 0, 0], + [-2, 0, -1], + ]); + let m3 = m1.add(m2).to2DArray(); + expect(m3).toStrictEqual([ + [2, 1, 6], + [2, 0, 3], + [0, 0, 0], + ]); + }); it('mmul', () => { let m1 = new SparseMatrix([ [2, 0, 1], @@ -20,6 +39,46 @@ describe('Sparse Matrix', () => { expect(m3.cardinality).toBe(1); expect(m3.get(0, 1)).toBe(2); + expectMatrixClose(m3.to2DArray(), [ + [0, 2], + [0, 0], + ]); + + // Compare with dense multiplication + const denseM1 = new Matrix(m1.to2DArray()); + const denseM2 = new Matrix(m2.to2DArray()); + const expectedDense = denseM1.mmul(denseM2); + expectMatrixClose(m3.to2DArray(), expectedDense.to2DArray()); + }); + + it('mmul', () => { + const size = 32; + const density = 0.1; + const A = randomMatrix(size, size, density * size ** 2); + const B = randomMatrix(size, size, density * size ** 2); + const m1 = new SparseMatrix(A); + const m2 = new SparseMatrix(B); + const m3 = m1.mmul(m2); + + const denseM1 = new Matrix(A); + const denseM2 = new Matrix(B); + const expectedDense = denseM1.mmul(denseM2); + expectMatrixClose(m3.to2DArray(), expectedDense.to2DArray()); + }); + + it('mmul with low density', () => { + const size = 128; + const cardinality = 64; + const A = randomMatrix(size, size, cardinality); + const B = randomMatrix(size, size, cardinality); + const m1 = new SparseMatrix(A); + const m2 = new SparseMatrix(B); + const m3 = m1.mmul(m2); + + const denseM1 = new Matrix(A); + const denseM2 = new Matrix(B); + const expectedDense = denseM1.mmul(denseM2); + expectMatrixClose(m3.to2DArray(), expectedDense.to2DArray()); }); it('kronecker', () => { @@ -122,3 +181,38 @@ describe('Banded matrices', () => { expect(matrix4.isBanded(1)).toBe(true); }); }); + +/** + * Helper to compare two 2D arrays element-wise using toBeCloseTo. + */ +function expectMatrixClose(received, expected, precision = 6) { + expect(received.length).toBe(expected.length); + for (let i = 0; i < received.length; i++) { + expect(received[i].length).toBe(expected[i].length); + for (let j = 0; j < received[i].length; j++) { + expect(received[i][j]).toBeCloseTo(expected[i][j], precision); + } + } +} + +function randomMatrix(rows, cols, cardinality) { + const total = rows * cols; + const positions = new Set(); + + // Generate unique random positions + while (positions.size < cardinality) { + positions.add(Math.floor(Math.random() * total)); + } + + // Build the matrix with zeros + const matrix = Array.from({ length: rows }, () => new Float64Array(cols)); + + // Assign random values to the selected positions + for (const pos of positions) { + const i = Math.floor(pos / cols); + const j = pos % cols; + matrix[i][j] = Math.random() * 10; + } + + return matrix; +} diff --git a/src/index.js b/src/index.js index f770b91..6c23c82 100644 --- a/src/index.js +++ b/src/index.js @@ -1,5 +1,11 @@ import HashTable from 'ml-hash-table'; +import { cooToCsr } from './utils/cooToCsr.js'; +import { mmulLowDensity } from './utils/mmulLowDensity.js'; +import { mmulMediumDensity } from './utils/mmulMediumDensity.js'; +import { mmulSmall } from './utils/mmulSmall.js'; + +/** @typedef {(row: number, column: number, value: number) => void} WithEachNonZeroCallback */ /** @typedef {(row: number, column: number, value: number) => number | false} ForEachNonZeroCallback */ export class SparseMatrix { @@ -18,11 +24,11 @@ export class SparseMatrix { if (Array.isArray(rows)) { const matrix = rows; - rows = matrix.length; + const nbRows = matrix.length; options = columns || {}; columns = matrix[0].length; - this._init(rows, columns, new HashTable(options), options.threshold); - for (let i = 0; i < rows; i++) { + this._init(nbRows, columns, new HashTable(options), options.threshold); + for (let i = 0; i < nbRows; i++) { for (let j = 0; j < columns; j++) { let value = matrix[i][j]; if (this.threshold && Math.abs(value) < this.threshold) value = 0; @@ -31,6 +37,7 @@ export class SparseMatrix { } } } + this.elements.maybeShrinkCapacity(); } else { this._init(rows, columns, new HashTable(options), options.threshold); } @@ -150,35 +157,30 @@ export class SparseMatrix { } else { this.elements.set(row * this.columns + column, value); } + return this; } /** + * Matrix multiplication, does not modify the current instance. * @param {SparseMatrix} other - * @returns {SparseMatrix} + * @returns {SparseMatrix} returns a new matrix instance. + * @throws {Error} If the number of columns of this matrix does not match the number of rows of the other matrix. */ mmul(other) { if (this.columns !== other.rows) { - // eslint-disable-next-line no-console - console.warn( + throw new RangeError( 'Number of columns of left matrix are not equal to number of rows of right matrix.', ); } - const m = this.rows; - const p = other.columns; + if (this.cardinality < 42 && other.cardinality < 42) { + return mmulSmall(this, other); + } else if (other.rows > 100 && other.cardinality < 100) { + return mmulLowDensity(this, other); + } - const result = new SparseMatrix(m, p); - this.forEachNonZero((i, j, v1) => { - other.forEachNonZero((k, l, v2) => { - if (j === k) { - result.set(i, l, result.get(i, l) + v1 * v2); - } - return v2; - }); - return v1; - }); - return result; + return mmulMediumDensity(this, other); } /** @@ -192,18 +194,63 @@ export class SparseMatrix { const q = other.columns; const result = new SparseMatrix(m * p, n * q, { - initialCapacity: this.cardinality * other.cardinality, - }); - this.forEachNonZero((i, j, v1) => { - other.forEachNonZero((k, l, v2) => { - result.set(p * i + k, q * j + l, v1 * v2); - return v2; - }); - return v1; + initialCapacity: this.cardinality * other.cardinality + 10, }); + + const { + columns: otherCols, + rows: otherRows, + values: otherValues, + } = other.getNonZeros(); + const { + columns: thisCols, + rows: thisRows, + values: thisValues, + } = this.getNonZeros(); + + const nbThisActive = thisCols.length; + const nbOtherActive = otherCols.length; + for (let t = 0; t < nbThisActive; t++) { + const pi = p * thisRows[t]; + const qj = q * thisCols[t]; + for (let o = 0; o < nbOtherActive; o++) { + result.set( + pi + otherRows[o], + qj + otherCols[o], + otherValues[o] * thisValues[t], + ); + } + } + return result; } + /** + * Calls `callback` for each value in the matrix that is not zero. + * Unlike `forEachNonZero`, the callback's return value has no effect. + * @param {WithEachNonZeroCallback} callback + * @param {boolean} [sort] + * @returns {void} + */ + withEachNonZero(callback, sort = false) { + const { state, table, values } = this.elements; + const nbStates = state.length; + const activeIndex = new Float64Array(this.cardinality); + for (let i = 0, j = 0; i < nbStates; i++) { + if (state[i] === 1) activeIndex[j++] = i; + } + + if (sort) { + activeIndex.sort((a, b) => table[a] - table[b]); + } + + const columns = this.columns; + for (const i of activeIndex) { + const key = table[i]; + callback((key / columns) | 0, key % columns, values[i]); + } + } + /** * Calls `callback` for each value in the matrix that is not zero. * The callback can return: @@ -234,23 +281,44 @@ export class SparseMatrix { return this; } - getNonZeros() { + /** + * Returns the non-zero elements of the matrix in coordinates (COO) or CSR format. + * + * **COO (Coordinate) format:** + * Stores the non-zero elements as three arrays: `rows`, `columns`, and `values`, where each index corresponds to a non-zero entry at (row, column) with the given value. + * + * **CSR (Compressed Sparse Row) format:** + * Stores the matrix using three arrays: + * - `rows`: Row pointer array of length `numRows + 1`, where each entry indicates the start of a row in the `columns` and `values` arrays. + * - `columns`: Column indices of non-zero elements. + * - `values`: Non-zero values. + * This format is efficient for row slicing and matrix-vector multiplication. + * + * @param {Object} [options={}] - Options for output format and sorting. + * @param {boolean} [options.csr] - If true, returns the result in CSR format. Otherwise, returns COO format. + * @param {boolean} [options.sort] - If true, sorts the non-zero elements by their indices. + * @returns {Object} If `csr` is not specified, returns an object with Float64Array `rows`, `columns`, and `values` (COO format). + * If `csr` is true, returns `{ rows, columns, values }` in CSR format. + */ + getNonZeros(options = {}) { const cardinality = this.cardinality; - /** @type {number[]} */ - const rows = new Array(cardinality); - /** @type {number[]} */ - const columns = new Array(cardinality); - /** @type {number[]} */ - const values = new Array(cardinality); + const rows = new Float64Array(cardinality); + const columns = new Float64Array(cardinality); + const values = new Float64Array(cardinality); + + const { csr, sort } = options; + let idx = 0; - this.forEachNonZero((i, j, value) => { + this.withEachNonZero((i, j, value) => { rows[idx] = i; columns[idx] = j; values[idx] = value; idx++; - return value; - }); - return { rows, columns, values }; + }, sort || csr); + + return csr + ? cooToCsr({ rows, columns, values }, this.rows) + : { rows, columns, values }; } /** @@ -272,9 +340,8 @@ export class SparseMatrix { let trans = new SparseMatrix(this.columns, this.rows, { initialCapacity: this.cardinality, }); - this.forEachNonZero((i, j, value) => { + this.withEachNonZero((i, j, value) => { trans.set(j, i, value); - return value; }); return trans; } diff --git a/src/utils/cooToCsr.js b/src/utils/cooToCsr.js new file mode 100644 index 0000000..0f4a2dd --- /dev/null +++ b/src/utils/cooToCsr.js @@ -0,0 +1,23 @@ +/** + * Converts a matrix from Coordinate (COO) format to Compressed Sparse Row (CSR) format. + * @param {Object} cooMatrix - The matrix in COO format with properties: values, columns, rows. + * @param {number} nbRows - The number of rows in the matrix. + * @returns {{rows: Float64Array, columns: Float64Array, values: Float64Array}} The matrix in CSR format. + */ +export function cooToCsr(cooMatrix, nbRows) { + const { values, columns, rows } = cooMatrix; + const csrRowPtr = new Float64Array(nbRows + 1); + + // Count non-zeros per row + const numberOfNonZeros = rows.length; + for (let i = 0; i < numberOfNonZeros; i++) { + csrRowPtr[rows[i] + 1]++; + } + + // Compute cumulative sum + for (let i = 1; i <= nbRows; i++) { + csrRowPtr[i] += csrRowPtr[i - 1]; + } + + return { rows: csrRowPtr, columns, values }; +} diff --git a/src/utils/mmulLowDensity.js b/src/utils/mmulLowDensity.js new file mode 100644 index 0000000..0310357 --- /dev/null +++ b/src/utils/mmulLowDensity.js @@ -0,0 +1,41 @@ +import { SparseMatrix } from '../index.js'; + +/** + * Multiplies two sparse matrices, optimized for cases where the right-hand side matrix + * has low density (number of rows > 100 and cardinality < 100). + * + * @private + * @param {SparseMatrix} left - The left-hand side matrix. + * @param {SparseMatrix} right - The right-hand side matrix to multiply with. + * @returns {SparseMatrix} The resulting matrix after multiplication. + */ +export function mmulLowDensity(left, right) { + const { + columns: otherCols, + rows: otherRows, + values: otherValues, + } = right.getNonZeros(); + const { + columns: thisCols, + rows: thisRows, + values: thisValues, + } = left.getNonZeros(); + + const m = left.rows; + const p = right.columns; + const output = new SparseMatrix(m, p); + + const nbOtherActive = otherCols.length; + const nbThisActive = thisCols.length; + for (let t = 0; t < nbThisActive; t++) { + const i = thisRows[t]; + const j = thisCols[t]; + for (let o = 0; o < nbOtherActive; o++) { + if (j === otherRows[o]) { + const l = otherCols[o]; + output.set(i, l, output.get(i, l) + otherValues[o] * thisValues[t]); + } + } + } + return output; +} diff --git a/src/utils/mmulMediumDensity.js b/src/utils/mmulMediumDensity.js new file mode 100644 index 0000000..2bee6f8 --- /dev/null +++ b/src/utils/mmulMediumDensity.js @@ -0,0 +1,41 @@ +import { SparseMatrix } from '../index.js'; + +/** + * Multiplies two sparse matrices where the right-hand side matrix is of medium density. + * Uses CSR format for the right-hand side for efficient multiplication. + * + * @private + * @param {SparseMatrix} left - The left-hand side matrix. + * @param {SparseMatrix} right - The right-hand side matrix to multiply with (in CSR format). + * @returns {SparseMatrix} The resulting matrix after multiplication. + */ + +export function mmulMediumDensity(left, right) { + const m = left.rows; + const p = right.columns; + const result = new SparseMatrix(m, p); + const { + columns: thisCols, + rows: thisRows, + values: thisValues, + } = left.getNonZeros(); + const { + columns: otherCols, + rows: otherRows, + values: otherValues, + } = right.getNonZeros({ csr: true }); + + const nbThisActive = thisCols.length; + for (let t = 0; t < nbThisActive; t++) { + const i = thisRows[t]; + const j = thisCols[t]; + const oStart = otherRows[j]; + const oEnd = otherRows[j + 1]; + for (let k = oStart; k < oEnd; k++) { + const l = otherCols[k]; + result.set(i, l, result.get(i, l) + otherValues[k] * thisValues[t]); + } + } + + return result; +} diff --git a/src/utils/mmulSmall.js b/src/utils/mmulSmall.js new file mode 100644 index 0000000..ba23481 --- /dev/null +++ b/src/utils/mmulSmall.js @@ -0,0 +1,33 @@ +import { SparseMatrix } from '../index.js'; + +/** + * Multiplies two very small sparse matrices (both with cardinalities < 42). + * + * @private + * @param {SparseMatrix} left - The left-hand side matrix. + * @param {SparseMatrix} right - The right-hand side matrix to multiply with. + * @returns {SparseMatrix} The resulting matrix after multiplication. + */ +export function mmulSmall(left, right) { + const { + columns: otherCols, + rows: otherRows, + values: otherValues, + } = right.getNonZeros(); + + const nbOtherActive = otherCols.length; + + const m = left.rows; + const p = right.columns; + const output = new SparseMatrix(m, p); + + left.withEachNonZero((i, j, v1) => { + for (let o = 0; o < nbOtherActive; o++) { + if (j === otherRows[o]) { + const l = otherCols[o]; + output.set(i, l, output.get(i, l) + otherValues[o] * v1); + } + } + }); + return output; +} diff --git a/vitest.config.js b/vitest.config.js new file mode 100644 index 0000000..27eb510 --- /dev/null +++ b/vitest.config.js @@ -0,0 +1,7 @@ +import { defineConfig } from 'vitest/config'; + +export default defineConfig({ + test: { + coverage: { include: ['src/**'] }, + }, +});