From 5aff896e00a78c155fabf6670f42af976200e4e9 Mon Sep 17 00:00:00 2001 From: "Rogelio J. Baucells" Date: Wed, 29 Apr 2015 17:54:17 -0400 Subject: [PATCH] element wise operations --- lib/function/arithmetic/add.js | 200 +----- lib/function/arithmetic/subtract.js | 196 +----- lib/type/matrix/util/elementWiseOperations.js | 590 ++++++++++++++++++ test/function/arithmetic/add.test.js | 36 +- test/function/arithmetic/subtract.test.js | 36 +- 5 files changed, 669 insertions(+), 389 deletions(-) create mode 100644 lib/type/matrix/util/elementWiseOperations.js diff --git a/lib/function/arithmetic/add.js b/lib/function/arithmetic/add.js index ec4192b56..0748da2ed 100644 --- a/lib/function/arithmetic/add.js +++ b/lib/function/arithmetic/add.js @@ -1,6 +1,5 @@ 'use strict'; -var clone = require('../../util/object').clone; var DimensionError = require('../../error/DimensionError'); function factory (type, config, load, typed) { @@ -8,12 +7,10 @@ function factory (type, config, load, typed) { var collection = load(require('../../type/collection')); var matrix = load(require('../construction/matrix')); var equal = load(require('../relational/equal')); - var sparseScatter = load(require('./sparseScatter')); var addScalar = load(require('./addScalar')); - var multiplyScalar = load(require('./multiplyScalar')); + var elementWiseOperations = load(require('../../type/matrix/util/elementWiseOperations')); - var DenseMatrix = type.DenseMatrix, - SparseMatrix = type.SparseMatrix; + var SparseMatrix = type.SparseMatrix; /** * Add two values, `x + y`. @@ -67,21 +64,22 @@ function factory (type, config, load, typed) { switch (y.storage()) { case 'sparse': // sparse + sparse - c = _addSparseMatrixSparseMatrix(x, y, xsize, ysize); + c = elementWiseOperations.algorithm4(x, y, addScalar, false); break; default: - c = _addSparseMatrixMatrix(x, y.valueOf(), xsize, ysize); + // sparse + dense + c = elementWiseOperations.algorithm1(y, x, addScalar, true); break; } break; default: switch (y.storage()) { case 'sparse': - // sparse + sparse - c = _addMatrixSparseMatrix(x.valueOf(), y, xsize, ysize); + // dense + sparse + c = elementWiseOperations.algorithm1(x, y, addScalar, false); break; default: - c = _addMatrixMatrix(x.valueOf(), y.valueOf(), x.storage()); + c = elementWiseOperations.algorithm0(x, y, addScalar, false); break; } break; @@ -197,188 +195,6 @@ function factory (type, config, load, typed) { return a.clone(); }; - /** - * C = A + B - * - * @param {Matrix} a SparseMatrix (MxN) - * @param {Matrix} b SparseMatrix (MxN) - * - * @return {Matrix} SparseMatrix (MxN) - */ - var _addSparseMatrixSparseMatrix = function (a, b, asize, bsize) { - // check dimensions - if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) - throw new RangeError('Dimension mismatch in add. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); - // rows and columns - var m = asize[0]; - var n = asize[1]; - // a arrays - var avalues = a._values; - var adt = a._datatype; - // b arrays - var bvalues = b._values; - var bdt = b._datatype; - // process data types - var dt = adt && bdt && adt === bdt ? adt : undefined; - // multiply scalar implementation - var mf = dt ? multiplyScalar.signatures[dt + ',' + dt] : multiplyScalar; - var af = dt ? addScalar.signatures[dt + ',' + dt] : addScalar; - // flag indicating both matrices (a & b) contain data - var values = avalues && bvalues; - // c arrays - var cvalues = values ? [] : undefined; - var cindex = []; - var cptr = new Array(n); - // c matrix - var c = new SparseMatrix({ - values: cvalues, - index: cindex, - ptr: cptr, - size: [m, n] - }); - // column vector (store matrix values) - var x = values ? new Array(m) : undefined; - // column vector to signal row values in column j - var w = new Array(m); - // loop columns - for (var j = 0; j < n; j++) { - // init ptr for j - cptr[j] = cindex.length; - // process column j of a and write it to x - sparseScatter(a, j, 1, w, x, j + 1, c, mf, af); - // process column j of b and write it to x - sparseScatter(b, j, 1, w, x, j + 1, c, mf, af); - // check matrix contains values (pattern matrix) - if (values) { - // loop column values in C - for (var p0 = cptr[j], p1 = cindex.length, p = p0; p < p1; p++) { - // copy x[i] to c[i, j] - cvalues.push(x[cindex[p]]); - } - } - } - // finish cptr - cptr[n] = cindex.length; - // return matrix - return c; - }; - - /** - * C = A + B - * - * @param {Matrix} a SparseMatrix (MxN) - * @param {Matrix} b DenseMatrix (MxN) - * - * @return {Matrix} SparseMatrix (MxN) - */ - var _addSparseMatrixMatrix = function (a, b, asize, bsize) { - // check dimensions - if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) - throw new RangeError('Dimension mismatch in add. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); - // rows and columns - var m = asize[0]; - var n = asize[1]; - // b array - var data = b; - // c arrays - var cvalues = []; - var cindex = []; - var cptr = new Array(n); - // c matrix - var c = new SparseMatrix({ - values: cvalues, - index: cindex, - ptr: cptr, - size: [m, n] - }); - // column vector (store matrix values) - var x = new Array(m); - // column vector to signal row values in column j - var w = new Array(m); - // loop columns - for (var j = 0; j < n; j++) { - // init ptr for j - cptr[j] = cindex.length; - // copy matrix b column to x - for (var i = 0; i < m; i++) { - // value - var v = data[i][j]; - // check for zero - if (!equal(v, 0)) { - x[i] = v; - w[i] = j + 1; - cindex.push(i); - } - } - // process column j of a and write it to x - sparseScatter(a, j, 1, w, x, j + 1, c, multiplyScalar, addScalar); - // loop column values in C - for (var p0 = cptr[j], p1 = cindex.length, p = p0; p < p1; p++) { - // copy x[i] to c[i, j] - cvalues.push(x[cindex[p]]); - } - } - // finish cptr - cptr[n] = cindex.length; - // return matrix - return c; - }; - - /** - * C = A + B - * - * @param {Matrix} a DenseMatrix (MxN) - * @param {Matrix} b SparseMatrix (MxN) - * - * @return {Matrix} DenseMatrix (MxN) - */ - var _addMatrixSparseMatrix = function (a, b, asize, bsize) { - // check dimensions - if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) - throw new RangeError('Dimension mismatch in add. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); - // rows and columns - var m = asize[0]; - var n = asize[1]; - // a array - var data = a; - // b arrays - var bvalues = b._values; - var bindex = b._index; - var bptr = b._ptr; - // c arrays - var cdata = clone(data); - // c matrix - var c = new DenseMatrix({ - data: cdata, - size: [m, n] - }); - // loop columns - for (var j = 0; j < n; j++) { - // loop values for column j - for (var k0 = bptr[j], k1 = bptr[j + 1], k = k0; k < k1; k++) { - // row - var i = bindex[k]; - // aggregate value - cdata[i][j] = addScalar(cdata[i][j], bvalues[k]); - } - } - // return matrix - return c; - }; - - /** - * C = A + B - * - * @param {Matrix} a DenseMatrix (MxN) - * @param {Matrix} b DenseMatrix (MxN) - * - * @return {Matrix} DenseMatrix (MxN) - */ - var _addMatrixMatrix = function (a, b, format) { - // TODO: find a better implementation - return matrix(collection.deepMap2(a, b, add), format); - }; - return add; } diff --git a/lib/function/arithmetic/subtract.js b/lib/function/arithmetic/subtract.js index a8df94982..ec9ce8b81 100644 --- a/lib/function/arithmetic/subtract.js +++ b/lib/function/arithmetic/subtract.js @@ -1,21 +1,16 @@ 'use strict'; -var clone = require('../../util/object').clone; var DimensionError = require('../../error/DimensionError'); function factory (type, config, load, typed) { - var unaryMinus = load(require('./unaryMinus')); var matrix = load(require('../construction/matrix')); var equal = load(require('../relational/equal')); - var sparseScatter = load(require('./sparseScatter')); - var addScalar = load(require('./addScalar')); - var multiplyScalar = load(require('./multiplyScalar')); + var elementWiseOperations = load(require('../../type/matrix/util/elementWiseOperations')); var collection = load(require('../../type/collection')); - var DenseMatrix = type.DenseMatrix, - SparseMatrix = type.SparseMatrix; + var SparseMatrix = type.SparseMatrix; /** * Subtract two values, `x - y`. @@ -105,21 +100,23 @@ function factory (type, config, load, typed) { switch (y.storage()) { case 'sparse': // sparse - sparse - c = _subtractSparseMatrixSparseMatrix(x, y, xsize, ysize); + c = elementWiseOperations.algorithm5(x, y, subtract); break; default: - c = _subtractSparseMatrixMatrix(x, y.valueOf(), xsize, ysize); + // sparse - dense + c = elementWiseOperations.algorithm3(y, x, subtract, true); break; } break; default: switch (y.storage()) { case 'sparse': - // sparse - sparse - c = _subtractMatrixSparseMatrix(x.valueOf(), y, xsize, ysize); + // dense - sparse + c = elementWiseOperations.algorithm1(x, y, subtract, false); break; default: - c = _subtractMatrixMatrix(x.valueOf(), y.valueOf(), x.storage()); + // dense - dense + c = elementWiseOperations.algorithm0(x, y, subtract, false); break; } break; @@ -289,181 +286,6 @@ function factory (type, config, load, typed) { return b.clone(); }; - /** - * C = A - B - * - * @param {Matrix} a SparseMatrix (MxN) - * @param {Matrix} b SparseMatrix (MxN) - * - * @return {Matrix} SparseMatrix (MxN) - */ - var _subtractSparseMatrixSparseMatrix = function (a, b, asize, bsize) { - // check dimensions - if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) - throw new RangeError('Dimension mismatch in add. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); - // rows and columns - var m = asize[0]; - var n = asize[1]; - // a arrays - var avalues = a._values; - // b arrays - var bvalues = b._values; - // flag indicating both matrices (a & b) contain data - var values = avalues && bvalues; - // c arrays - var cvalues = values ? [] : undefined; - var cindex = []; - var cptr = new Array(n); - // c matrix - var c = new SparseMatrix({ - values: cvalues, - index: cindex, - ptr: cptr, - size: [m, n] - }); - // column vector (store matrix values) - var x = values ? new Array(m) : undefined; - // column vector to signal row values in column j - var w = new Array(m); - // loop columns - for (var j = 0; j < n; j++) { - // init ptr for j - cptr[j] = cindex.length; - // process column j of a and write it to x - sparseScatter(a, j, 1, w, x, j + 1, c, multiplyScalar, addScalar); - // process column j of b and write it to x (multiply value by negative one) - sparseScatter(b, j, -1, w, x, j + 1, c, multiplyScalar, addScalar); - // check matrix contains values (pattern matrix) - if (values) { - // loop column values in C - for (var p0 = cptr[j], p1 = cindex.length, p = p0; p < p1; p++) { - // copy x[i] to c[i, j] - cvalues.push(x[cindex[p]]); - } - } - } - // finish cptr - cptr[n] = cindex.length; - // return matrix - return c; - }; - - /** - * C = A - B - * - * @param {Matrix} a SparseMatrix (MxN) - * @param {Matrix} b DenseMatrix (MxN) - * - * @return {Matrix} SparseMatrix (MxN) - */ - var _subtractSparseMatrixMatrix = function (a, b, asize, bsize) { - // check dimensions - if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) - throw new RangeError('Dimension mismatch in add. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); - // rows and columns - var m = asize[0]; - var n = asize[1]; - // b array - var data = b; - // c arrays - var cvalues = []; - var cindex = []; - var cptr = new Array(n); - // c matrix - var c = new SparseMatrix({ - values: cvalues, - index: cindex, - ptr: cptr, - size: [m, n] - }); - // column vector (store matrix values) - var x = new Array(m); - // column vector to signal row values in column j - var w = new Array(m); - // loop columns - for (var j = 0; j < n; j++) { - // init ptr for j - cptr[j] = cindex.length; - // copy matrix b column to x - for (var i = 0; i < m; i++) { - // -value - var v = unaryMinus(data[i][j]); - // check for zero - if (!equal(v, 0)) { - x[i] = v; - w[i] = j + 1; - cindex.push(i); - } - } - // process column j of a and write it to x - sparseScatter(a, j, 1, w, x, j + 1, c, multiplyScalar, addScalar); - // loop column values in C - for (var p0 = cptr[j], p1 = cindex.length, p = p0; p < p1; p++) { - // copy x[i] to c[i, j] - cvalues.push(x[cindex[p]]); - } - } - // finish cptr - cptr[n] = cindex.length; - // return matrix - return c; - }; - - /** - * C = A - B - * - * @param {Matrix} a DenseMatrix (MxN) - * @param {Matrix} b SparseMatrix (MxN) - * - * @return {Matrix} DenseMatrix (MxN) - */ - var _subtractMatrixSparseMatrix = function (a, b, asize, bsize) { - // check dimensions - if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) - throw new RangeError('Dimension mismatch in add. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); - // rows and columns - var m = asize[0]; - var n = asize[1]; - // a array - var data = a; - // b arrays - var bvalues = b._values; - var bindex = b._index; - var bptr = b._ptr; - // c arrays - var cdata = clone(data); - // c matrix - var c = new DenseMatrix({ - data: cdata, - size: [m, n] - }); - // loop columns - for (var j = 0; j < n; j++) { - // loop values for column j - for (var k0 = bptr[j], k1 = bptr[j + 1], k = k0; k < k1; k++) { - // row - var i = bindex[k]; - // subtract value - cdata[i][j] = subtract(cdata[i][j], bvalues[k]); - } - } - // return matrix - return c; - }; - - /** - * C = A - B - * - * @param {Matrix} a DenseMatrix (MxN) - * @param {Matrix} b DenseMatrix (MxN) - * - * @return {Matrix} DenseMatrix (MxN) - */ - var _subtractMatrixMatrix = function (a, b, format) { - // TODO: find a better implementation - return matrix(collection.deepMap2(a, b, subtract), format); - }; - return subtract; } diff --git a/lib/type/matrix/util/elementWiseOperations.js b/lib/type/matrix/util/elementWiseOperations.js new file mode 100644 index 000000000..73fa7327d --- /dev/null +++ b/lib/type/matrix/util/elementWiseOperations.js @@ -0,0 +1,590 @@ +'use strict'; + +var util = require('../../../util/index'); +var DimensionError = require('../../../error/DimensionError'); + +var object = util.object; + +function factory (type, config, load) { + + var equal = load(require('../../../function/relational/equal')); + var collection = load(require('../../collection')); + + var DenseMatrix = type.DenseMatrix, + SparseMatrix = type.SparseMatrix; + + /** + * Iterates over DenseMatrix items and invokes the callback function f(Aij, Bij). + * Callback function invoked MxN times. + * + * + * C(i,j) = f(Aij, Bij) + * + * + * @param {Matrix} a The DenseMatrix instance (A) + * @param {Matrix} b The DenseMatrix instance (B) + * @param {function} callback The f(Aij,Bij) operation to invoke + * @param {boolean} inverse A true value indicates callback should be invoked f(Bij,Aij) + * + * @return {Matrix} DenseMatrix (C) + * + * see https://github.com/josdejong/mathjs/pull/346#issuecomment-97477571 + */ + var algorithm0 = function (a, b, callback) { + // TODO: look for a better algorithm here! implement "inverse" option + return collection.deepMap2(a, b, callback); + }; + + /** + * Iterates over SparseMatrix nonzero items and invokes the callback function f(Dij, Sij). + * Callback function invoked NNZ times (number of nonzero items in SparseMatrix). + * + * + * ┌ f(Dij, Sij) ; S(i,j) !== 0 + * C(i,j) = ┤ + * └ Dij ; otherwise + * + * + * @param {Matrix} denseMatrix The DenseMatrix instance (D) + * @param {Matrix} sparseMatrix The SparseMatrix instance (S) + * @param {function} callback The f(Dij,Sij) operation to invoke, where Dij = DenseMatrix(i,j) and Sij = SparseMatrix(i,j) + * @param {boolean} inverse A true value indicates callback should be invoked f(Sij,Dij) + * + * @return {Matrix} DenseMatrix (C) + * + * see https://github.com/josdejong/mathjs/pull/346#issuecomment-97477571 + */ + var algorithm1 = function (denseMatrix, sparseMatrix, callback, inverse) { + // dense matrix arrays + var adata = denseMatrix._data; + var asize = denseMatrix._size; + var adt = denseMatrix._datatype; + // sparse matrix arrays + var bvalues = sparseMatrix._values; + var bindex = sparseMatrix._index; + var bptr = sparseMatrix._ptr; + var bsize = sparseMatrix._size; + var bdt = sparseMatrix._datatype; + + // validate dimensions + if (asize.length !== bsize.length) + throw new DimensionError(asize.length, bsize.length); + + // check rows & columns + if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) + throw new RangeError('Dimension mismatch. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); + + // sparse matrix cannot be a Pattern matrix + if (!bvalues) + throw new Error('Cannot perform operation on Dense Matrix and Pattern Sparse Matrix'); + + // rows & columns + var rows = asize[0]; + var columns = asize[1]; + + // process data types + var dt = adt && bdt && adt === bdt ? adt : undefined; + // callback implementation + var cf = dt && callback.signatures ? callback.signatures[dt + ',' + dt] || callback : callback; + + // result (DenseMatrix) + var cdata = object.clone(adata); + + // loop columns in b + for (var j = 0; j < columns; j++) { + // values in column j + for (var k0 = bptr[j], k1 = bptr[j + 1], k = k0; k < k1; k++) { + // row + var i = bindex[k]; + // update C(i,j) + cdata[i][j] = inverse ? cf(bvalues[k], cdata[i][j]) : cf(cdata[i][j], bvalues[k]); + } + } + + // return dense matrix + return new DenseMatrix({ + data: cdata, + size: [rows, columns], + datatype: dt + }); + }; + + /** + * Iterates over SparseMatrix nonzero items and invokes the callback function f(Dij, Sij). + * Callback function invoked NNZ times (number of nonzero items in SparseMatrix). + * + * + * ┌ f(Dij, Sij) ; S(i,j) !== 0 + * C(i,j) = ┤ + * └ 0 ; otherwise + * + * + * @param {Matrix} denseMatrix The DenseMatrix instance (D) + * @param {Matrix} sparseMatrix The SparseMatrix instance (S) + * @param {function} callback The f(Dij,Sij) operation to invoke, where Dij = DenseMatrix(i,j) and Sij = SparseMatrix(i,j) + * @param {boolean} inverse A true value indicates callback should be invoked f(Sij,Dij) + * + * @return {Matrix} SparseMatrix (C) + * + * see https://github.com/josdejong/mathjs/pull/346#issuecomment-97477571 + */ + var algorithm2 = function (denseMatrix, sparseMatrix, callback, inverse) { + // dense matrix arrays + var adata = denseMatrix._data; + var asize = denseMatrix._size; + var adt = denseMatrix._datatype; + // sparse matrix arrays + var bvalues = sparseMatrix._values; + var bindex = sparseMatrix._index; + var bptr = sparseMatrix._ptr; + var bsize = sparseMatrix._size; + var bdt = sparseMatrix._datatype; + + // validate dimensions + if (asize.length !== bsize.length) + throw new DimensionError(asize.length, bsize.length); + + // check rows & columns + if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) + throw new RangeError('Dimension mismatch. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); + + // sparse matrix cannot be a Pattern matrix + if (!bvalues) + throw new Error('Cannot perform operation on Dense Matrix and Pattern Sparse Matrix'); + + // rows & columns + var rows = asize[0]; + var columns = asize[1]; + + // process data types + var dt = adt && bdt && adt === bdt ? adt : undefined; + // callback implementation + var cf = dt && callback.signatures ? callback.signatures[dt + ',' + dt] || callback : callback; + + // result (SparseMatrix) + var cvalues = []; + var cindex = []; + var cptr = new Array(columns + 1); + + // loop columns in b + for (var j = 0; j < columns; j++) { + // update cptr + cptr[j] = cindex.length; + // values in column j + for (var k0 = bptr[j], k1 = bptr[j + 1], k = k0; k < k1; k++) { + // row + var i = bindex[k]; + // update C(i,j) + var cij = inverse ? cf(bvalues[k], adata[i][j]) : cf(adata[i][j], bvalues[k]); + // check for nonzero + if (!equal(cij, 0)) { + // push i & v + cindex.push(i); + cvalues.push(cij); + } + } + } + // update cptr + cptr[columns] = cindex.length; + + // return sparse matrix + return new SparseMatrix({ + values: cvalues, + index: cindex, + ptr: cptr, + size: [rows, columns], + datatype: dt + }); + }; + + /** + * Iterates over SparseMatrix items and invokes the callback function f(Dij, Sij). + * Callback function invoked M*N times. + * + * + * ┌ f(Dij, Sij) ; S(i,j) !== 0 + * C(i,j) = ┤ + * └ f(Dij, 0) ; otherwise + * + * + * @param {Matrix} denseMatrix The DenseMatrix instance (D) + * @param {Matrix} sparseMatrix The SparseMatrix instance (C) + * @param {function} callback The f(Dij,Sij) operation to invoke, where Dij = DenseMatrix(i,j) and Sij = SparseMatrix(i,j) + * @param {boolean} inverse A true value indicates callback should be invoked f(Sij,Dij) + * + * @return {Matrix} DenseMatrix (C) + * + * see https://github.com/josdejong/mathjs/pull/346#issuecomment-97477571 + */ + var algorithm3 = function (denseMatrix, sparseMatrix, callback, inverse) { + // dense matrix arrays + var adata = denseMatrix._data; + var asize = denseMatrix._size; + var adt = denseMatrix._datatype; + // sparse matrix arrays + var bvalues = sparseMatrix._values; + var bindex = sparseMatrix._index; + var bptr = sparseMatrix._ptr; + var bsize = sparseMatrix._size; + var bdt = sparseMatrix._datatype; + var bzero = sparseMatrix._zero; + + // validate dimensions + if (asize.length !== bsize.length) + throw new DimensionError(asize.length, bsize.length); + + // check rows & columns + if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) + throw new RangeError('Dimension mismatch. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); + + // sparse matrix cannot be a Pattern matrix + if (!bvalues) + throw new Error('Cannot perform operation on Dense Matrix and Pattern Sparse Matrix'); + + // rows & columns + var rows = asize[0]; + var columns = asize[1]; + + // process data types (zero value is required!) + var dt = adt && bdt && typeof(bzero) !== 'undefined' && bzero !== null && adt === bdt ? adt : undefined; + // callback implementation + var cf = dt && callback.signatures ? callback.signatures[dt + ',' + dt] || callback : callback; + + // sparse matrix zero value + var zero = dt ? bzero : 0; + + // result (DenseMatrix) + var cdata = new Array(rows); + + // initialize dense matrix + for (var z = 0; z < rows; z++) { + // initialize row + cdata[z] = new Array(columns); + } + + // workspace + var x = new Array(rows); + // marks indicating we have a value in x for a given column + var w = new Array(rows); + + // loop columns in b + for (var j = 0; j < columns; j++) { + // column mark + var mark = j + 1; + // values in column j + for (var k0 = bptr[j], k1 = bptr[j + 1], k = k0; k < k1; k++) { + // row + var i = bindex[k]; + // update workspace + x[i] = inverse ? cf(bvalues[k], adata[i][j]) : cf(adata[i][j], bvalues[k]); + w[i] = mark; + } + // process workspace + for (var y = 0; y < rows; y++) { + // check we have a calculated value for current row + if (w[y] === mark) { + // use calculated value + cdata[y][j] = x[y]; + } + else { + // calculate value + cdata[y][j] = inverse ? cf(zero, adata[y][j]) : cf(adata[y][j], zero); + } + } + } + + // return dense matrix + return new DenseMatrix({ + data: cdata, + size: [rows, columns], + datatype: dt + }); + }; + + /** + * Iterates over SparseMatrix A and SparseMatrix B nonzero items and invokes the callback function f(Aij, Bij). + * Callback function invoked MAX(NNZA, NNZB) times + * + * + * ┌ f(Aij, Bij) ; A(i,j) !== 0 && B(i,j) !== 0 + * C(i,j) = ┤ A(i,j) ; A(i,j) !== 0 + * └ B(i,j) ; B(i,j) !== 0 + * + * + * @param {Matrix} a The SparseMatrix instance (A) + * @param {Matrix} b The SparseMatrix instance (B) + * @param {function} callback The f(Aij,Bij) operation to invoke + * @param {boolean} inverse A true value indicates callback should be invoked f(Bij,Aij) + * + * @return {Matrix} SparseMatrix (C) + * + * see https://github.com/josdejong/mathjs/pull/346#issuecomment-97477571 + */ + var algorithm4 = function (a, b, callback, inverse) { + // sparse matrix arrays + var avalues = a._values; + var asize = a._size; + var adt = a._datatype; + var azero = a._zero; + // sparse matrix arrays + var bvalues = b._values; + var bsize = b._size; + var bdt = b._datatype; + var bzero = b._zero; + + // validate dimensions + if (asize.length !== bsize.length) + throw new DimensionError(asize.length, bsize.length); + + // check rows & columns + if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) + throw new RangeError('Dimension mismatch. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); + + // rows & columns + var rows = asize[0]; + var columns = asize[1]; + + // process data types + var dt = adt && bdt && typeof(azero) !== 'undefined' && azero !== null && typeof(bzero) !== 'undefined' && bzero !== null && adt === bdt ? adt : undefined; + // callback implementation + var cf = dt && callback.signatures ? callback.signatures[dt + ',' + dt] || callback : callback; + // equal + var ef = dt ? equal.signatures[dt + ',' + dt] || equal : equal; + + // sparse matrix zero value + var zero = dt ? azero : 0; + + // result arrays + var cvalues = avalues && bvalues ? [] : undefined; + var cindex = []; + var cptr = new Array(columns + 1); + // matrix + var c = new SparseMatrix({ + values: cvalues, + index: cindex, + ptr: cptr, + size: [rows, columns], + datatype: dt + }); + + // workspace + var x = avalues && bvalues ? new Array(rows) : undefined; + // marks indicating we have a value in x for a given column + var w = new Array(rows); + + // loop columns + for (var j = 0; j < columns; j++) { + // update cptr + cptr[j] = cindex.length; + // columns mark + var mark = j + 1; + // scatter the values of A(:,j) into workspace + _scatter(a, j, w, x, mark, c, cf, inverse); + // scatter the values of B(:,j) into workspace + _scatter(b, j, w, x, mark, c, cf, inverse); + // check we need to process values (non pattern matrix) + if (x) { + // initialize first index in j + var k = cptr[j]; + // loop index in j + while (k < cindex.length) { + // row + var i = cindex[k]; + // value @ i + var v = x[i]; + // check for zero value + if (!ef(v, zero)) { + // push value + cvalues.push(v); + // increment pointer + k++; + } + else { + // remove value @ i, do not increment pointer + cindex.splice(k, 1); + } + } + } + } + // update cptr + cptr[columns] = cindex.length; + + // return sparse matrix + return c; + }; + + /** + * Iterates over SparseMatrix A and SparseMatrix B nonzero items and invokes the callback function f(Aij, Bij). + * Callback function invoked MAX(NNZA, NNZB) times + * + * + * ┌ f(Aij, Bij) ; A(i,j) !== 0 || B(i,j) !== 0 + * C(i,j) = ┤ + * └ 0 ; otherwise + * + * + * @param {Matrix} a The SparseMatrix instance (A) + * @param {Matrix} b The SparseMatrix instance (B) + * @param {function} callback The f(Aij,Bij) operation to invoke + * + * @return {Matrix} SparseMatrix (C) + * + * see https://github.com/josdejong/mathjs/pull/346#issuecomment-97477571 + */ + var algorithm5 = function (a, b, callback) { + // sparse matrix arrays + var avalues = a._values; + var asize = a._size; + var adt = a._datatype; + var azero = a._zero; + // sparse matrix arrays + var bvalues = b._values; + var bsize = b._size; + var bdt = b._datatype; + var bzero = b._zero; + + // validate dimensions + if (asize.length !== bsize.length) + throw new DimensionError(asize.length, bsize.length); + + // check rows & columns + if (asize[0] !== bsize[0] || asize[1] !== bsize[1]) + throw new RangeError('Dimension mismatch. Matrix A (' + asize + ') must match Matrix B (' + bsize + ')'); + + // rows & columns + var rows = asize[0]; + var columns = asize[1]; + + // process data types (zero value is required!) + var dt = adt && bdt && typeof(azero) !== 'undefined' && azero !== null && typeof(bzero) !== 'undefined' && bzero !== null && adt === bdt ? adt : undefined; + // callback implementation + var cf = dt && callback.signatures ? callback.signatures[dt + ',' + dt] || callback : callback; + // equal + var ef = dt ? equal.signatures[dt + ',' + dt] || equal : equal; + + // sparse matrix zero value + var zero = dt ? azero : 0; + + // result arrays + var cvalues = avalues && bvalues ? [] : undefined; + var cindex = []; + var cptr = new Array(columns + 1); + // matrix + var c = new SparseMatrix({ + values: cvalues, + index: cindex, + ptr: cptr, + size: [rows, columns], + datatype: dt, + zero: zero + }); + + // workspaces + var x = cvalues ? new Array(rows) : undefined; + // marks indicating we have a value in x for a given column + var w = new Array(rows); + + // loop columns + for (var j = 0; j < columns; j++) { + // update cptr + cptr[j] = cindex.length; + // columns mark + var mark = j + 1; + // scatter the values of A(:,j) into workspace + _scatter(a, j, w, x, mark, c); + // scatter the values of B(:,j) into workspace + _scatter(b, j, w, x, mark, c, cf, false, true, zero); + // check we need to process values (non pattern matrix) + if (x) { + // initialize first index in j + var k = cptr[j]; + // loop index in j + while (k < cindex.length) { + // row + var i = cindex[k]; + // value @ i + var v = x[i]; + // check for zero value + if (!ef(v, zero)) { + // push value + cvalues.push(v); + // increment pointer + k++; + } + else { + // remove value @ i, do not increment pointer + cindex.splice(k, 1); + } + } + } + } + // update cptr + cptr[columns] = cindex.length; + + // return sparse matrix + return c; + }; + + var _scatter = function (a, j, w, x, mark, c, f, inverse, update, value) { + // a arrays + var avalues = a._values; + var aindex = a._index; + var aptr = a._ptr; + // c arrays + var cindex = c._index; + + // vars + var k, k0, k1, i; + + // check we need to process values (pattern matrix) + if (x) { + // values in j + for (k0 = aptr[j], k1 = aptr[j + 1], k = k0; k < k1; k++) { + // row + i = aindex[k]; + // check value exists in current j + if (w[i] !== mark) { + // i is new entry in j + w[i] = mark; + // add i to pattern of C + cindex.push(i); + // x(i) = A + x[i] = update ? (inverse ? f(avalues[k], value) : f(value, avalues[k])) : avalues[k]; + } + else { + // i exists in C already + x[i] = inverse ? f(avalues[k], x[i]) : f(x[i], avalues[k]); + } + } + } + else { + // values in j + for (k0 = aptr[j], k1 = aptr[j + 1], k = k0; k < k1; k++) { + // row + i = aindex[k]; + // check value exists in current j + if (w[i] !== mark) { + // i is new entry in j + w[i] = mark; + // add i to pattern of C + cindex.push(i); + } + } + } + // number of nonzero elements in C + return cindex.length; + }; + + // return algorithms + return { + algorithm0: algorithm0, + algorithm1: algorithm1, + algorithm2: algorithm2, + algorithm3: algorithm3, + algorithm4: algorithm4, + algorithm5: algorithm5 + }; +} + +exports.name = 'elementWiseOperations'; +exports.factory = factory; diff --git a/test/function/arithmetic/add.test.js b/test/function/arithmetic/add.test.js index 30404da11..e0e27da4a 100644 --- a/test/function/arithmetic/add.test.js +++ b/test/function/arithmetic/add.test.js @@ -117,7 +117,7 @@ describe('add', function() { assert.deepEqual(add([3,4], 2), [5,6]); }); - it('should add array and matrix correctly', function() { + it('should add array and dense matrix correctly', function() { var a = [1,2,3]; var b = math.matrix([3,2,1]); var c = add(a, b); @@ -125,6 +125,15 @@ describe('add', function() { assert.ok(c instanceof math.type.Matrix); assert.deepEqual(c, math.matrix([4,4,4])); }); + + it('should add array and sparse matrix correctly', function() { + var a = [[1,2,3],[4,5,6]]; + var b = math.sparse([[6, 5, 4],[ 3, 2, 1]]); + var c = add(a, b); + + assert.ok(c instanceof math.type.Matrix); + assert.deepEqual(c, math.matrix([[7,7,7],[7,7,7]])); + }); }); describe('DenseMatrix', function () { @@ -156,17 +165,25 @@ describe('add', function() { assert.ok(c instanceof math.type.Matrix); assert.deepEqual(c, math.matrix([4,4,4])); }); + + it('should add dense and sparse matrices correctly', function() { + var a = math.matrix([[1,2,3],[1,0,0]]); + var b = math.sparse([[3,2,1],[0,0,1]]); + var c = add(a, b); + + assert.ok(c instanceof math.type.Matrix); + assert.deepEqual(c, math.matrix([[4,4,4],[1,0,1]])); + }); }); describe('SparseMatrix', function () { it('should add matrices correctly', function() { var a2 = math.matrix([[1,2],[3,4]], 'sparse'); - var a3 = math.matrix([[5,6],[7,8]], 'sparse'); + var a3 = math.matrix([[5,-2],[7,-4]], 'sparse'); var a4 = add(a2, a3); assert.ok(a4 instanceof math.type.Matrix); - assert.deepEqual(a4.size(), [2,2]); - assert.deepEqual(a4.valueOf(), [[6,8],[10,12]]); + assert.deepEqual(a4, math.sparse([[6,0],[10,0]])); }); it('should add a scalar and a matrix correctly', function() { @@ -180,7 +197,16 @@ describe('add', function() { var c = add(a, b); assert.ok(c instanceof math.type.Matrix); - assert.deepEqual(c, math.matrix([[4,4,4],[1,0,1]], 'sparse')); + assert.deepEqual(c, math.matrix([[4,4,4],[1,0,1]])); + }); + + it('should add sparse and dense matrices correctly', function() { + var a = math.sparse([[1,2,3],[1,0,0]]); + var b = math.matrix([[3,2,1],[0,0,1]]); + var c = add(a, b); + + assert.ok(c instanceof math.type.Matrix); + assert.deepEqual(c, math.matrix([[4,4,4],[1,0,1]])); }); it('should add two pattern matrices correctly', function() { diff --git a/test/function/arithmetic/subtract.test.js b/test/function/arithmetic/subtract.test.js index 2db289d06..95f856add 100644 --- a/test/function/arithmetic/subtract.test.js +++ b/test/function/arithmetic/subtract.test.js @@ -125,7 +125,7 @@ describe('subtract', function() { assert.deepEqual(subtract([3,4], 2), [1,2]); }); - it('should subtract array and matrix correctly', function() { + it('should subtract array and dense matrix correctly', function() { var a = [1,2,3]; var b = math.matrix([3,2,1]); var c = subtract(a, b); @@ -133,6 +133,15 @@ describe('subtract', function() { assert.ok(c instanceof math.type.Matrix); assert.deepEqual(c, math.matrix([-2,0,2])); }); + + it('should subtract array and dense matrix correctly', function() { + var a = [[1,2,3],[4,5,6]]; + var b = math.sparse([[6,5,4],[ 3, 2, 1]]); + var c = subtract(a, b); + + assert.ok(c instanceof math.type.Matrix); + assert.deepEqual(c, math.matrix([[-5,-3,-1],[1,3,5]])); + }); }); describe('DenseMatrix', function () { @@ -159,17 +168,25 @@ describe('subtract', function() { assert.ok(c instanceof math.type.Matrix); assert.deepEqual(c, math.matrix([-2,0,2])); }); + + it('should subtract dense and sparse matrices correctly', function() { + var a = math.matrix([[1,2,3],[1,0,0]]); + var b = math.sparse([[3,2,1],[0,0,1]]); + var c = subtract(a, b); + + assert.ok(c instanceof math.type.Matrix); + assert.deepEqual(c, math.matrix([[-2,0,2],[1,0,-1]])); + }); }); describe('SparseMatrix', function () { it('should subtract matrices correctly', function() { - var a2 = math.matrix([[10,20],[30,40]], 'sparse'); - var a3 = math.matrix([[5,6],[7,8]], 'sparse'); + var a2 = math.matrix([[10,20],[30,0]], 'sparse'); + var a3 = math.matrix([[5,6],[30,8]], 'sparse'); var a4 = subtract(a2, a3); assert.ok(a4 instanceof math.type.Matrix); - assert.deepEqual(a4.size(), [2,2]); - assert.deepEqual(a4.valueOf(), [[5,14],[23,32]]); + assert.deepEqual(a4, math.sparse([[5,14],[0,-8]])); }); it('should subtract a scalar and a matrix correctly', function() { @@ -185,6 +202,15 @@ describe('subtract', function() { assert.ok(c instanceof math.type.Matrix); assert.deepEqual(c.valueOf(), [[-2,0,2],[1,0,-1]]); }); + + it('should subtract sparse and dense matrices correctly', function() { + var a = math.sparse([[1,2,3],[1,0,0]]); + var b = math.matrix([[3,2,1],[0,0,1]]); + var c = subtract(a, b); + + assert.ok(c instanceof math.type.Matrix); + assert.deepEqual(c, math.matrix([[-2,0,2],[1,0,-1]])); + }); }); it('should throw an error in case of invalid number of arguments', function() {