[Pkg-javascript-commits] [science.js] 28/87: Add new vector module and test for lin.tridag.
bhuvan krishna
bhuvan-guest at moszumanska.debian.org
Thu Dec 8 06:11:54 UTC 2016
This is an automated email from the git hooks/post-receive script.
bhuvan-guest pushed a commit to branch master
in repository science.js.
commit bb57b8e623d79dcf48679a678111dfec48744954
Author: Jason Davies <jason at jasondavies.com>
Date: Fri Aug 26 09:14:27 2011 +0100
Add new vector module and test for lin.tridag.
---
Makefile | 13 ++++
science.js | 153 ++++++++++++++++++++++++++++++++++++++++++++++
science.lin.js | 14 ++---
science.lin.min.js | 2 +-
science.min.js | 2 +-
src/lin/tridag.js | 14 ++---
src/vector/cross.js | 9 +++
src/vector/determinant.js | 17 ++++++
src/vector/dot.js | 7 +++
src/vector/gaussjordan.js | 55 +++++++++++++++++
src/vector/inverse.js | 26 ++++++++
src/vector/length.js | 3 +
src/vector/multiply.js | 19 ++++++
src/vector/normalize.js | 4 ++
src/vector/transpose.js | 12 ++++
src/vector/vector.js | 1 +
test/lin/tridag-test.js | 40 ++++++++++++
17 files changed, 373 insertions(+), 18 deletions(-)
diff --git a/Makefile b/Makefile
index 182c9b2..a65ed9d 100644
--- a/Makefile
+++ b/Makefile
@@ -13,6 +13,7 @@ all: \
.INTERMEDIATE science.js: \
src/start.js \
science.core.js \
+ science.vector.js \
src/end.js
science.core.js: \
@@ -21,6 +22,18 @@ science.core.js: \
src/core/hypot.js \
src/core/zeroes.js
+science.vector.js: \
+ src/vector/vector.js \
+ src/vector/cross.js \
+ src/vector/dot.js \
+ src/vector/length.js \
+ src/vector/normalize.js \
+ src/vector/determinant.js \
+ src/vector/gaussjordan.js \
+ src/vector/inverse.js \
+ src/vector/multiply.js \
+ src/vector/transpose.js
+
science.lin.js: \
src/start.js \
src/lin/lin.js \
diff --git a/science.js b/science.js
index 84589d3..058aa85 100644
--- a/science.js
+++ b/science.js
@@ -27,4 +27,157 @@ science.zeroes = function(n) {
this, Array.prototype.slice.call(arguments, 1));
return a;
};
+science.vector = {};
+science.vector.cross = function(a, b) {
+ // TODO how to handle non-3D vectors?
+ // TODO handle 7D vectors?
+ return [
+ a[1] * b[2] - a[2] * b[1],
+ a[2] * b[0] - a[0] * b[2],
+ a[0] * b[1] - a[1] * b[0]
+ ];
+};
+science.vector.dot = function(a, b) {
+ var s = 0,
+ i = -1,
+ n = Math.min(a.length, b.length);
+ while (++i < n) s += a[i] * b[i];
+ return s;
+};
+science.vector.length = function(p) {
+ return Math.sqrt(science.vector.dot(p, p));
+};
+science.vector.normalize = function(p) {
+ var length = science.vector.length(p);
+ return p.map(function(d) { return d / length; });
+};
+// 4x4 matrix determinant.
+science.vector.determinant = function(matrix) {
+ var m = matrix[0].concat(matrix[1]).concat(matrix[2]).concat(matrix[3]);
+ return (
+ m[12] * m[9] * m[6] * m[3] - m[8] * m[13] * m[6] * m[3] -
+ m[12] * m[5] * m[10] * m[3] + m[4] * m[13] * m[10] * m[3] +
+ m[8] * m[5] * m[14] * m[3] - m[4] * m[9] * m[14] * m[3] -
+ m[12] * m[9] * m[2] * m[7] + m[8] * m[13] * m[2] * m[7] +
+ m[12] * m[1] * m[10] * m[7] - m[0] * m[13] * m[10] * m[7] -
+ m[8] * m[1] * m[14] * m[7] + m[0] * m[9] * m[14] * m[7] +
+ m[12] * m[5] * m[2] * m[11] - m[4] * m[13] * m[2] * m[11] -
+ m[12] * m[1] * m[6] * m[11] + m[0] * m[13] * m[6] * m[11] +
+ m[4] * m[1] * m[14] * m[11] - m[0] * m[5] * m[14] * m[11] -
+ m[8] * m[5] * m[2] * m[15] + m[4] * m[9] * m[2] * m[15] +
+ m[8] * m[1] * m[6] * m[15] - m[0] * m[9] * m[6] * m[15] -
+ m[4] * m[1] * m[10] * m[15] + m[0] * m[5] * m[10] * m[15]);
+};
+// Performs in-place Gauss-Jordan elimination.
+//
+// Based on Jarno Elonen's Python version (public domain):
+// http://elonen.iki.fi/code/misc-notes/python-gaussj/index.html
+science.vector.gaussjordan = function(m, eps) {
+ if (!eps) eps = 1e-10;
+
+ var h = m.length,
+ w = m[0].length,
+ y = -1,
+ y2,
+ x;
+
+ while (++y < h) {
+ var maxrow = y;
+
+ // Find max pivot.
+ y2 = y; while (++y2 < h) {
+ if (Math.abs(m[y2][y]) > Math.abs(m[maxrow][y]))
+ maxrow = y2;
+ }
+
+ // Swap.
+ var tmp = m[y];
+ m[y] = m[maxrow];
+ m[maxrow] = tmp;
+
+ // Singular?
+ if (Math.abs(m[y][y]) <= eps) return false;
+
+ // Eliminate column y.
+ y2 = y; while (++y2 < h) {
+ var c = m[y2][y] / m[y][y];
+ x = y - 1; while (++x < w) {
+ m[y2][x] -= m[y][x] * c;
+ }
+ }
+ }
+
+ // Backsubstitute.
+ y = h; while (--y >= 0) {
+ var c = m[y][y];
+ y2 = -1; while (++y2 < y) {
+ x = w; while (--x >= y) {
+ m[y2][x] -= m[y][x] * m[y2][y] / c;
+ }
+ }
+ m[y][y] /= c;
+ // Normalize row y.
+ x = h - 1; while (++x < w) {
+ m[y][x] /= c;
+ }
+ }
+ return true;
+};
+// Find matrix inverse using Gauss-Jordan.
+science.vector.inverse = function(m) {
+ var n = m.length
+ i = -1;
+
+ // Check if the matrix is square.
+ if (n !== m[0].length) return;
+
+ // Augment with identity matrix I to get AI.
+ m = m.map(function(row, i) {
+ var identity = new Array(n),
+ j = -1;
+ while (++j < n) identity[j] = i === j ? 1 : 0;
+ return row.concat(identity);
+ });
+
+ // Compute IA^-1.
+ science.vector.gaussjordan(m);
+
+ // Remove identity matrix I to get A^-1.
+ while (++i < n) {
+ m[i] = m[i].slice(n);
+ }
+
+ return m;
+};
+science.vector.multiply = function(a, b) {
+ var m = a.length,
+ n = b[0].length,
+ p = b.length,
+ i = -1,
+ j,
+ k;
+ if (p !== a[0].length) throw {"error": "columns(a) != rows(b); " + a[0].length + " != " + p};
+ var ab = new Array(m);
+ while (++i < m) {
+ ab[i] = new Array(n);
+ j = -1; while(++j < n) {
+ var s = 0;
+ k = -1; while (++k < p) s += a[i][k] * b[k][j];
+ ab[i][j] = s;
+ }
+ }
+ return ab;
+};
+science.vector.transpose = function(a) {
+ var m = a.length,
+ n = a[0].length,
+ i = -1,
+ j,
+ b = new Array(n);
+ while (++i < n) {
+ b[i] = new Array(m);
+ j = -1; while (++j < m) b[i][j] = a[j][i];
+ }
+ return b;
+};
})()
\ No newline at end of file
diff --git a/science.lin.js b/science.lin.js
index e29a079..cbdc2ac 100644
--- a/science.lin.js
+++ b/science.lin.js
@@ -12,18 +12,16 @@
* @param {number} n
*/
science.lin.tridag = function(a, b, c, d, x, n) {
- c[0] /= b[0];
- d[0] /= b[0];
var i,
- id;
+ m;
for (i = 1; i < n; i++) {
- id = 1 / (b[i] - c[i-1] * a[i]);
- c[i] *= id;
- d[i] = (d[i] - d[i - 1] * a[i]) * id;
+ m = a[i] / b[i - 1];
+ b[i] -= m * c[i - 1];
+ d[i] -= m * d[i - 1];
}
- x[n] = d[n];
+ x[n - 1] = d[n - 1] / b[n - 1];
for (i = n - 2; i >= 0; i--) {
- x[i] = d[i] - c[i] * x[i + 1];
+ x[i] = (d[i] - c[i] * x[i + 1]) / b[i];
}
};
})()
\ No newline at end of file
diff --git a/science.lin.min.js b/science.lin.min.js
index 53458fd..a470da4 100644
--- a/science.lin.min.js
+++ b/science.lin.min.js
@@ -1 +1 @@
-(function(){science.lin={},science.lin.tridag=function(a,b,c,d,e,f){c[0]/=b[0],d[0]/=b[0];var g,h;for(g=1;g<f;g++)h=1/(b[g]-c[g-1]*a[g]),c[g]*=h,d[g]=(d[g]-d[g-1]*a[g])*h;e[f]=d[f];for(g=f-2;g>=0;g--)e[g]=d[g]-c[g]*e[g+1]}})()
\ No newline at end of file
+(function(){science.lin={},science.lin.tridag=function(a,b,c,d,e,f){var g,h;for(g=1;g<f;g++)h=a[g]/b[g-1],b[g]-=h*c[g-1],d[g]-=h*d[g-1];e[f-1]=d[f-1]/b[f-1];for(g=f-2;g>=0;g--)e[g]=(d[g]-c[g]*e[g+1])/b[g]}})()
\ No newline at end of file
diff --git a/science.min.js b/science.min.js
index e8dd89d..a68af5b 100644
--- a/science.min.js
+++ b/science.min.js
@@ -1 +1 @@
-(function(){science={version:"1.5.0"},science.functor=function(a){return typeof a=="function"?a:function(){return a}},science.hypot=function(a,b){a=Math.abs(a),b=Math.abs(b);var c,d;a>b?(c=a,d=b):(c=b,d=a);var e=d/c;return c*Math.sqrt(1+e*e)},science.zeroes=function(a){var b=-1,c=[];if(arguments.length===1)while(++b<a)c[b]=0;else while(++b<a)c[b]=science.zeroes.apply(this,Array.prototype.slice.call(arguments,1));return c}})()
\ No newline at end of file
+(function(){science={version:"1.5.0"},science.functor=function(a){return typeof a=="function"?a:function(){return a}},science.hypot=function(a,b){a=Math.abs(a),b=Math.abs(b);var c,d;a>b?(c=a,d=b):(c=b,d=a);var e=d/c;return c*Math.sqrt(1+e*e)},science.zeroes=function(a){var b=-1,c=[];if(arguments.length===1)while(++b<a)c[b]=0;else while(++b<a)c[b]=science.zeroes.apply(this,Array.prototype.slice.call(arguments,1));return c},science.vector={},science.vector.cross=function(a,b){return[a[1]*b [...]
\ No newline at end of file
diff --git a/src/lin/tridag.js b/src/lin/tridag.js
index 7b58072..a24a95f 100644
--- a/src/lin/tridag.js
+++ b/src/lin/tridag.js
@@ -11,17 +11,15 @@
* @param {number} n
*/
science.lin.tridag = function(a, b, c, d, x, n) {
- c[0] /= b[0];
- d[0] /= b[0];
var i,
- id;
+ m;
for (i = 1; i < n; i++) {
- id = 1 / (b[i] - c[i-1] * a[i]);
- c[i] *= id;
- d[i] = (d[i] - d[i - 1] * a[i]) * id;
+ m = a[i] / b[i - 1];
+ b[i] -= m * c[i - 1];
+ d[i] -= m * d[i - 1];
}
- x[n] = d[n];
+ x[n - 1] = d[n - 1] / b[n - 1];
for (i = n - 2; i >= 0; i--) {
- x[i] = d[i] - c[i] * x[i + 1];
+ x[i] = (d[i] - c[i] * x[i + 1]) / b[i];
}
};
diff --git a/src/vector/cross.js b/src/vector/cross.js
new file mode 100644
index 0000000..97cd71f
--- /dev/null
+++ b/src/vector/cross.js
@@ -0,0 +1,9 @@
+science.vector.cross = function(a, b) {
+ // TODO how to handle non-3D vectors?
+ // TODO handle 7D vectors?
+ return [
+ a[1] * b[2] - a[2] * b[1],
+ a[2] * b[0] - a[0] * b[2],
+ a[0] * b[1] - a[1] * b[0]
+ ];
+};
diff --git a/src/vector/determinant.js b/src/vector/determinant.js
new file mode 100644
index 0000000..b79e76a
--- /dev/null
+++ b/src/vector/determinant.js
@@ -0,0 +1,17 @@
+// 4x4 matrix determinant.
+science.vector.determinant = function(matrix) {
+ var m = matrix[0].concat(matrix[1]).concat(matrix[2]).concat(matrix[3]);
+ return (
+ m[12] * m[9] * m[6] * m[3] - m[8] * m[13] * m[6] * m[3] -
+ m[12] * m[5] * m[10] * m[3] + m[4] * m[13] * m[10] * m[3] +
+ m[8] * m[5] * m[14] * m[3] - m[4] * m[9] * m[14] * m[3] -
+ m[12] * m[9] * m[2] * m[7] + m[8] * m[13] * m[2] * m[7] +
+ m[12] * m[1] * m[10] * m[7] - m[0] * m[13] * m[10] * m[7] -
+ m[8] * m[1] * m[14] * m[7] + m[0] * m[9] * m[14] * m[7] +
+ m[12] * m[5] * m[2] * m[11] - m[4] * m[13] * m[2] * m[11] -
+ m[12] * m[1] * m[6] * m[11] + m[0] * m[13] * m[6] * m[11] +
+ m[4] * m[1] * m[14] * m[11] - m[0] * m[5] * m[14] * m[11] -
+ m[8] * m[5] * m[2] * m[15] + m[4] * m[9] * m[2] * m[15] +
+ m[8] * m[1] * m[6] * m[15] - m[0] * m[9] * m[6] * m[15] -
+ m[4] * m[1] * m[10] * m[15] + m[0] * m[5] * m[10] * m[15]);
+};
diff --git a/src/vector/dot.js b/src/vector/dot.js
new file mode 100644
index 0000000..74dfd00
--- /dev/null
+++ b/src/vector/dot.js
@@ -0,0 +1,7 @@
+science.vector.dot = function(a, b) {
+ var s = 0,
+ i = -1,
+ n = Math.min(a.length, b.length);
+ while (++i < n) s += a[i] * b[i];
+ return s;
+};
diff --git a/src/vector/gaussjordan.js b/src/vector/gaussjordan.js
new file mode 100644
index 0000000..08a8ca1
--- /dev/null
+++ b/src/vector/gaussjordan.js
@@ -0,0 +1,55 @@
+// Performs in-place Gauss-Jordan elimination.
+//
+// Based on Jarno Elonen's Python version (public domain):
+// http://elonen.iki.fi/code/misc-notes/python-gaussj/index.html
+science.vector.gaussjordan = function(m, eps) {
+ if (!eps) eps = 1e-10;
+
+ var h = m.length,
+ w = m[0].length,
+ y = -1,
+ y2,
+ x;
+
+ while (++y < h) {
+ var maxrow = y;
+
+ // Find max pivot.
+ y2 = y; while (++y2 < h) {
+ if (Math.abs(m[y2][y]) > Math.abs(m[maxrow][y]))
+ maxrow = y2;
+ }
+
+ // Swap.
+ var tmp = m[y];
+ m[y] = m[maxrow];
+ m[maxrow] = tmp;
+
+ // Singular?
+ if (Math.abs(m[y][y]) <= eps) return false;
+
+ // Eliminate column y.
+ y2 = y; while (++y2 < h) {
+ var c = m[y2][y] / m[y][y];
+ x = y - 1; while (++x < w) {
+ m[y2][x] -= m[y][x] * c;
+ }
+ }
+ }
+
+ // Backsubstitute.
+ y = h; while (--y >= 0) {
+ var c = m[y][y];
+ y2 = -1; while (++y2 < y) {
+ x = w; while (--x >= y) {
+ m[y2][x] -= m[y][x] * m[y2][y] / c;
+ }
+ }
+ m[y][y] /= c;
+ // Normalize row y.
+ x = h - 1; while (++x < w) {
+ m[y][x] /= c;
+ }
+ }
+ return true;
+};
diff --git a/src/vector/inverse.js b/src/vector/inverse.js
new file mode 100644
index 0000000..0a8bb5e
--- /dev/null
+++ b/src/vector/inverse.js
@@ -0,0 +1,26 @@
+// Find matrix inverse using Gauss-Jordan.
+science.vector.inverse = function(m) {
+ var n = m.length
+ i = -1;
+
+ // Check if the matrix is square.
+ if (n !== m[0].length) return;
+
+ // Augment with identity matrix I to get AI.
+ m = m.map(function(row, i) {
+ var identity = new Array(n),
+ j = -1;
+ while (++j < n) identity[j] = i === j ? 1 : 0;
+ return row.concat(identity);
+ });
+
+ // Compute IA^-1.
+ science.vector.gaussjordan(m);
+
+ // Remove identity matrix I to get A^-1.
+ while (++i < n) {
+ m[i] = m[i].slice(n);
+ }
+
+ return m;
+};
diff --git a/src/vector/length.js b/src/vector/length.js
new file mode 100644
index 0000000..49618f5
--- /dev/null
+++ b/src/vector/length.js
@@ -0,0 +1,3 @@
+science.vector.length = function(p) {
+ return Math.sqrt(science.vector.dot(p, p));
+};
diff --git a/src/vector/multiply.js b/src/vector/multiply.js
new file mode 100644
index 0000000..dc29173
--- /dev/null
+++ b/src/vector/multiply.js
@@ -0,0 +1,19 @@
+science.vector.multiply = function(a, b) {
+ var m = a.length,
+ n = b[0].length,
+ p = b.length,
+ i = -1,
+ j,
+ k;
+ if (p !== a[0].length) throw {"error": "columns(a) != rows(b); " + a[0].length + " != " + p};
+ var ab = new Array(m);
+ while (++i < m) {
+ ab[i] = new Array(n);
+ j = -1; while(++j < n) {
+ var s = 0;
+ k = -1; while (++k < p) s += a[i][k] * b[k][j];
+ ab[i][j] = s;
+ }
+ }
+ return ab;
+};
diff --git a/src/vector/normalize.js b/src/vector/normalize.js
new file mode 100644
index 0000000..3f6bc36
--- /dev/null
+++ b/src/vector/normalize.js
@@ -0,0 +1,4 @@
+science.vector.normalize = function(p) {
+ var length = science.vector.length(p);
+ return p.map(function(d) { return d / length; });
+};
diff --git a/src/vector/transpose.js b/src/vector/transpose.js
new file mode 100644
index 0000000..48a61db
--- /dev/null
+++ b/src/vector/transpose.js
@@ -0,0 +1,12 @@
+science.vector.transpose = function(a) {
+ var m = a.length,
+ n = a[0].length,
+ i = -1,
+ j,
+ b = new Array(n);
+ while (++i < n) {
+ b[i] = new Array(m);
+ j = -1; while (++j < m) b[i][j] = a[j][i];
+ }
+ return b;
+};
diff --git a/src/vector/vector.js b/src/vector/vector.js
new file mode 100644
index 0000000..fe140eb
--- /dev/null
+++ b/src/vector/vector.js
@@ -0,0 +1 @@
+science.vector = {};
diff --git a/test/lin/tridag-test.js b/test/lin/tridag-test.js
new file mode 100644
index 0000000..ed04362
--- /dev/null
+++ b/test/lin/tridag-test.js
@@ -0,0 +1,40 @@
+require("../../science");
+require("../../science.lin");
+
+var vows = require("vows"),
+ assert = require("assert");
+
+var suite = vows.describe("science.lin.tridag");
+
+suite.addBatch({
+ "tridag": {
+ "simple": function() {
+ var n = 31;
+ var a = [], b = [], c = [], d = [], x = [];
+ for (var i=0; i<n; i++) {
+ a[i] = 1;
+ b[i] = 2;
+ c[i] = 1;
+ d[i] = i < n / 2 ? i + 1 : n - i;
+ }
+ science.lin.tridag(a, b.slice(), c, d.slice(), x, n);
+ var matrix = [];
+ for (var i=0; i<n; i++) {
+ matrix[i] = [];
+ for (var j=0; j<n; j++) matrix[i][j] = 0;
+ }
+ for (var i=0; i<n; i++) {
+ matrix[i][i] = b[i];
+ if (i > 0) matrix[i - 1][i] = a[i];
+ if (i < n - 1) matrix[i + 1][i] = c[i];
+ }
+ var result = science.vector.multiply(matrix, x.map(function(d) { return [d]; }));
+ var epsilon = 1e-12;
+ for (var i=0; i<n; i++) {
+ assert.isTrue(result[i] - d[i] < epsilon);
+ }
+ }
+ }
+});
+
+suite.export(module);
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-javascript/science.js.git
More information about the Pkg-javascript-commits
mailing list