diff --git a/web/app/constr/solver.js b/web/app/constr/solver.js index ef74a982..e7cd8f27 100644 --- a/web/app/constr/solver.js +++ b/web/app/constr/solver.js @@ -207,6 +207,15 @@ TCAD.parametric.lock2Equals2 = function(constrs, locked) { TCAD.parametric._alg = 5; +TCAD.parametric.diagnose = function(sys) { + var jacobian = sys.makeJacobian(); + var qr = new TCAD.math.QR(jacobian); + return { + conflict : sys.constraints.length > qr.rank, + dof : sys.params.length - qr.rank + } +}; + TCAD.parametric.prepare = function(constrs, locked, aux, alg) { this.lock1(constrs, aux); @@ -234,10 +243,15 @@ TCAD.parametric.prepare = function(constrs, locked, aux, alg) { return sys.makeJacobian(); }; alg = TCAD.parametric._alg; + var _point = []; function solve(fineLevel) { if (constrs.length == 0) return; if (sys.params.length == 0) return; + if (TCAD.parametric.diagnose(sys).conflict) { + console.log("Conflicting or redundant constraints. Please fix your system."); + return; + } if (alg > 0) { switch (alg) { case 1: @@ -255,8 +269,8 @@ TCAD.parametric.prepare = function(constrs, locked, aux, alg) { break; case 5: if (optim.dog_leg(sys) !== 0) { - alg = -5; - solve(fineLevel); + //alg = -5; + //solve(fineLevel); } break; } diff --git a/web/app/math/qr.js b/web/app/math/qr.js index 9f51b555..cd940f29 100644 --- a/web/app/math/qr.js +++ b/web/app/math/qr.js @@ -1,24 +1,33 @@ -TCAD.math.qr = function(matrix) { - function arr(size) { - var out = []; - out.length = size; - for (var i = 0; i < size; ++i) { - out[i] = 0; - } - return out; +TCAD.math.vec = function(size) { + var out = []; + out.length = size; + for (var i = 0; i < size; ++i) { + out[i] = 0; } + return out; +}; + + +TCAD.math.Arrays_fill = function(a, fromIndex, toIndex,val) { + for (var i = fromIndex; i < toIndex; i++) + a[i] = val; +}; + +TCAD.math.QR = function(matrix) { + var vec = TCAD.math.vec; + this.matrix = matrix; + var nR = this.matrix.length; + var nC = this.matrix[0].length; this.qrRankingThreshold = 1e-30; //?? this.solvedCols = Math.min(nR, nC); - this.diagR = arr(nC); - this.norm = arr(nC); - this.beta = arr(nC); - this.permutation = arr(nC); + this.diagR = vec(nC); + this.norm = vec(nC); + this.beta = vec(nC); + this.permutation = vec(nC); this.rank = null; - var nR = this.matrix.length; - var nC = this.matrix[0].length; var k; var norm2; var akk; @@ -88,3 +97,131 @@ TCAD.math.qr = function(matrix) { } this.rank = this.solvedCols; }; + +TCAD.math.QR.prototype.qTy = function(y) { + var nR = this.matrix.length; + var nC = this.matrix[0].length; + + for (var k = 0; k < nC; ++k) { + var pk = this.permutation[k]; + var gamma = 0; + for (var i = k; i < nR; ++i) { + gamma += this.matrix[i][pk] * y[i]; + } + gamma *= this.beta[pk]; + for (var i = k; i < nR; ++i) { + y[i] -= gamma * this.matrix[i][pk]; + } + } +}; + +TCAD.math.QR.prototype.solve = function(qy) { + + var nR = this.matrix.length; + var nC = this.matrix[0].length; + + var vec = TCAD.math.vec; + + var diag = vec(nC); + var lmDiag = vec(nC); + var work = vec(nC); + var out = vec(nC); + + // copy R and Qty to preserve input and initialize s + // in particular, save the diagonal elements of R in lmDir + for (var j = 0; j < this.solvedCols; ++j) { + var pj = this.permutation[j]; + for (var i = j + 1; i < this.solvedCols; ++i) { + this.matrix[i][pj] = this.matrix[j][this.permutation[i]]; + } + out[j] = this.diagR[pj]; + work[j] = qy[j]; + } + + // eliminate the diagonal matrix d using a Givens rotation + for (var j = 0; j < this.solvedCols; ++j) { + + // prepare the row of d to be eliminated, locating the + // diagonal element using p from the Q.R. factorization + var pj = this.permutation[j]; + var dpj = diag[pj]; + if (dpj != 0) { + TCAD.math.Arrays_fill(lmDiag, j + 1, lmDiag.length, 0); + } + lmDiag[j] = dpj; + + // the transformations to eliminate the row of d + // modify only a single element of Qty + // beyond the first n, which is initially zero. + var qtbpj = 0; + for (var k = j; k < this.solvedCols; ++k) { + var pk = this.permutation[k]; + + // determine a Givens rotation which eliminates the + // appropriate element in the current row of d + if (lmDiag[k] != 0) { + + var sin; + var cos; + var rkk = this.matrix[k][pk]; + if (Math.abs(rkk) < Math.abs(lmDiag[k])) { + var cotan = rkk / lmDiag[k]; + sin = 1.0 / Math.sqrt(1.0 + cotan * cotan); + cos = sin * cotan; + } else { + var tan = lmDiag[k] / rkk; + cos = 1.0 / Math.sqrt(1.0 + tan * tan); + sin = cos * tan; + } + + // compute the modified diagonal element of R and + // the modified element of (Qty,0) + this.matrix[k][pk] = cos * rkk + sin * lmDiag[k]; + var temp = cos * work[k] + sin * qtbpj; + qtbpj = -sin * work[k] + cos * qtbpj; + work[k] = temp; + + // accumulate the tranformation in the row of s + for (var i = k + 1; i < this.solvedCols; ++i) { + var rik = this.matrix[i][pk]; + var temp2 = cos * rik + sin * lmDiag[i]; + lmDiag[i] = -sin * rik + cos * lmDiag[i]; + this.matrix[i][pk] = temp2; + } + } + } + + // store the diagonal element of s and restore + // the corresponding diagonal element of R + lmDiag[j] = this.matrix[j][this.permutation[j]]; + this.matrix[j][this.permutation[j]] = out[j]; + } + + // solve the triangular system for z, if the system is + // singular, then obtain a least squares solution + var nSing = this.solvedCols; + for (var j = 0; j < this.solvedCols; ++j) { + if ((lmDiag[j] == 0) && (nSing == this.solvedCols)) { + nSing = j; + } + if (nSing < this.solvedCols) { + work[j] = 0; + } + } + if (nSing > 0) { + for (var j = nSing - 1; j >= 0; --j) { + var pj = this.permutation[j]; + var sum = 0; + for (var i = j + 1; i < nSing; ++i) { + sum += this.matrix[i][pj] * work[i]; + } + work[j] = (work[j] - sum) / lmDiag[j]; + } + } + + // permute the components of z back to components of lmDir + for (var j = 0; j < out.length; ++j) { + out[this.permutation[j]] = work[j]; + } + return out; +}; \ No newline at end of file diff --git a/web/app/math/test.js b/web/app/math/test.js index d025fe89..57e90b01 100644 --- a/web/app/math/test.js +++ b/web/app/math/test.js @@ -55,13 +55,19 @@ function testCompare() { function lsolve() { var n = numeric; - A = [[1, 0, 0], [2,3,0]] - b = [10,25] + var A = [[2,3], [1,1]]; + var b = [13,5]; // var At = n.transpose(A); // var A = n.dot(At, A); // var b = n.dot(At, b); // var res = n.dot(n.dot(At, n.inv(n.dot(A, At)) ), b); - var res = optim.cg(A, [100, 100, 0], b, 1e-8, 800); +// var res = optim.cg(A, [100, 100, 0], b, 1e-8, 800); + + //var qr = new TCAD.math.QR(A); + //qr.qTy(b); + //var res = qr.solve(b); + + var res = n.solve(A, b); console.log(res); } diff --git a/web/sketcher.html b/web/sketcher.html index 791aa5a4..9717e0d7 100644 --- a/web/sketcher.html +++ b/web/sketcher.html @@ -192,6 +192,7 @@ + diff --git a/web/test.html b/web/test.html index 36caa05e..a92ead56 100644 --- a/web/test.html +++ b/web/test.html @@ -1,27 +1,44 @@ - Test Suite + TCAD + + + - - - - - - - - - + + - - - - + + + + + + - + + + + + + + +