Template:Team:Amsterdam/maarten/scripts/models

From 2012.igem.org

Revision as of 20:20, 23 July 2012 by MaartenR (Talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

// // Math.seedrandom('yipee'); Sets Math.random to a function that is // initialized using the given explicit seed. // // Math.seedrandom(); Sets Math.random to a function that is // seeded using the current time, dom state, // and other accumulated local entropy. // The generated seed string is returned. // // Math.seedrandom('yowza', true); // Seeds using the given explicit seed mixed // together with accumulated entropy. // // // Seeds using physical random bits downloaded // from random.org. // // Seeds using urandom bits from call.jsonlib.com, // which is faster than random.org. // // Examples: // // Math.seedrandom("hello"); // Use "hello" as the seed. // document.write(Math.random()); // Always 0.5463663768140734 // document.write(Math.random()); // Always 0.43973793770592234 // var rng1 = Math.random; // Remember the current prng. // // var autoseed = Math.seedrandom(); // New prng with an automatic seed. // document.write(Math.random()); // Pretty much unpredictable. // // Math.random = rng1; // Continue "hello" prng sequence. // document.write(Math.random()); // Always 0.554769432473455 // // Math.seedrandom(autoseed); // Restart at the previous seed. // document.write(Math.random()); // Repeat the 'unpredictable' value. // // Notes: // // Each time seedrandom('arg') is called, entropy from the passed seed // is accumulated in a pool to help generate future seeds for the // zero-argument form of Math.seedrandom, so entropy can be injected over // time by calling seedrandom with explicit data repeatedly. // // On speed - This javascript implementation of Math.random() is about // 3-10x slower than the built-in Math.random() because it is not native // code, but this is typically fast enough anyway. Seeding is more expensive, // especially if you use auto-seeding. Some details (timings on Chrome 4): // // Our Math.random() - avg less than 0.002 milliseconds per call // seedrandom('explicit') - avg less than 0.5 milliseconds per call // seedrandom('explicit', true) - avg less than 2 milliseconds per call // seedrandom() - avg about 38 milliseconds per call // // LICENSE (BSD): // // Copyright 2010 David Bau, all rights reserved. // // Redistribution and use in source and binary forms, with or without // modification, are permitted provided that the following conditions are met: // // 1. Redistributions of source code must retain the above copyright // notice, this list of conditions and the following disclaimer. // // 2. Redistributions in binary form must reproduce the above copyright // notice, this list of conditions and the following disclaimer in the // documentation and/or other materials provided with the distribution. // // 3. Neither the name of this module nor the names of its contributors may // be used to endorse or promote products derived from this software // without specific prior written permission. // // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // /** * All code is in an anonymous closure to keep the global namespace clean. * * @param {number=} overflow * @param {number=} startdenom */ // Patched by Seb so that seedrandom.js does not pollute the Math object. // My tests suggest that doing Math.trouble = 1 makes Math lookups about 5% // slower. numeric.seedrandom = { pow:Math.pow, random:Math.random }; (function (pool, math, width, chunks, significance, overflow, startdenom) { // // seedrandom() // This is the seedrandom function described above. // math['seedrandom'] = function seedrandom(seed, use_entropy) { var key = []; var arc4; // Flatten the seed string or build one from local entropy if needed. seed = mixkey(flatten( use_entropy ? [seed, pool] : arguments.length ? seed : [new Date().getTime(), pool, window], 3), key); // Use the seed to initialize an ARC4 generator. arc4 = new ARC4(key); // Mix the randomness into accumulated entropy. mixkey(arc4.S, pool); // Override Math.random // This function returns a random double in [0, 1) that contains // randomness in every bit of the mantissa of the IEEE 754 value. math['random'] = function random() { // Closure to return a random double: var n = arc4.g(chunks); // Start with a numerator n < 2 ^ 48 var d = startdenom; // and denominator d = 2 ^ 48. var x = 0; // and no 'extra last byte'. while (n < significance) { // Fill up all significant digits by n = (n + x) * width; // shifting numerator and d *= width; // denominator and generating a x = arc4.g(1); // new least-significant-byte. } while (n >= overflow) { // To avoid rounding up, before adding n /= 2; // last byte, shift everything d /= 2; // right using integer math until x >>>= 1; // we have exactly the desired bits. } return (n + x) / d; // Form the number within [0, 1). }; // Return the seed that was used return seed; }; // // ARC4 // // An ARC4 implementation. The constructor takes a key in the form of // an array of at most (width) integers that should be 0 <= x < (width). // // The g(count) method returns a pseudorandom integer that concatenates // the next (count) outputs from ARC4. Its return value is a number x // that is in the range 0 <= x < (width ^ count). // /** @constructor */ function ARC4(key) { var t, u, me = this, keylen = key.length; var i = 0, j = me.i = me.j = me.m = 0; me.S = []; me.c = []; // The empty key [] is treated as [0]. if (!keylen) { key = [keylen++]; } // Set up S using the standard key scheduling algorithm. while (i < width) { me.S[i] = i++; } for (i = 0; i < width; i++) { t = me.S[i]; j = lowbits(j + t + key[i % keylen]); u = me.S[j]; me.S[i] = u; me.S[j] = t; } // The "g" method returns the next (count) outputs as one number. me.g = function getnext(count) { var s = me.S; var i = lowbits(me.i + 1); var t = s[i]; var j = lowbits(me.j + t); var u = s[j]; s[i] = u; s[j] = t; var r = s[lowbits(t + u)]; while (--count) { i = lowbits(i + 1); t = s[i]; j = lowbits(j + t); u = s[j]; s[i] = u; s[j] = t; r = r * width + s[lowbits(t + u)]; } me.i = i; me.j = j; return r; }; // For robust unpredictability discard an initial batch of values. // See http://www.rsa.com/rsalabs/node.asp?id=2009 me.g(width); } // // flatten() // Converts an object tree to nested arrays of strings. // /** @param {Object=} result * @param {string=} prop * @param {string=} typ */ function flatten(obj, depth, result, prop, typ) { result = []; typ = typeof(obj); if (depth && typ == 'object') { for (prop in obj) { if (prop.indexOf('S') < 5) { // Avoid FF3 bug (local/sessionStorage) try { result.push(flatten(obj[prop], depth - 1)); } catch (e) {} } } } return (result.length ? result : obj + (typ != 'string' ? '\0' : '')); } // // mixkey() // Mixes a string seed into a key that is an array of integers, and // returns a shortened string seed that is equivalent to the result key. // /** @param {number=} smear * @param {number=} j */ function mixkey(seed, key, smear, j) { seed += ''; // Ensure the seed is a string smear = 0; for (j = 0; j < seed.length; j++) { key[lowbits(j)] = lowbits((smear ^= key[lowbits(j)] * 19) + seed.charCodeAt(j)); } seed = ''; for (j in key) { seed += String.fromCharCode(key[j]); } return seed; } // // lowbits() // A quick "n mod width" for width a power of 2. // function lowbits(n) { return n & (width - 1); } // // The following constants are related to IEEE 754 limits. // startdenom = math.pow(width, chunks); significance = math.pow(2, significance); overflow = significance * 2; // // When seedrandom.js is loaded, we immediately mix a few bits // from the built-in RNG into the entropy pool. Because we do // not want to intefere with determinstic PRNG state later, // seedrandom will not call math.random on its own again after // initialization. // mixkey(math.random(), pool); // End anonymous scope, and pass initial values. }( [], // pool: entropy pool starts empty numeric.seedrandom, // math: package containing random, pow, and seedrandom 256, // width: each RC4 output is 0 <= x < 256 6, // chunks: at least six RC4 outputs for each double 52 // significance: there are 52 significant digits in a double )); /* This file is a slightly modified version of quadprog.js from Alberto Santini. * It has been slightly modified by Sébastien Loisel to make sure that it handles * 0-based Arrays instead of 1-based Arrays. * License is in resources/LICENSE.quadprog */ (function(exports) { function base0to1(A) { if(typeof A !== "object") { return A; } var ret = [], i,n=A.length; for(i=0;i meq) { work[l] = sum; } else { work[l] = -Math.abs(sum); if (sum > 0) { for (j = 1; j <= n; j = j + 1) { amat[j][i] = -amat[j][i]; } bvec[i] = -bvec[i]; } } } for (i = 1; i <= nact; i = i + 1) { work[iwsv + iact[i]] = 0; } nvl = 0; temp = 0; for (i = 1; i <= q; i = i + 1) { if (work[iwsv + i] < temp * work[iwnbv + i]) { nvl = i; temp = work[iwsv + i] / work[iwnbv + i]; } } if (nvl === 0) { return 999; } return 0; } function fn_goto_55() { for (i = 1; i <= n; i = i + 1) { sum = 0; for (j = 1; j <= n; j = j + 1) { sum = sum + dmat[j][i] * amat[j][nvl]; } work[i] = sum; } l1 = iwzv; for (i = 1; i <= n; i = i + 1) { work[l1 + i] = 0; } for (j = nact + 1; j <= n; j = j + 1) { for (i = 1; i <= n; i = i + 1) { work[l1 + i] = work[l1 + i] + dmat[i][j] * work[j]; } } t1inf = true; for (i = nact; i >= 1; i = i - 1) { sum = work[i]; l = iwrm + (i * (i + 3)) / 2; l1 = l - i; for (j = i + 1; j <= nact; j = j + 1) { sum = sum - work[l] * work[iwrv + j]; l = l + j; } sum = sum / work[l1]; work[iwrv + i] = sum; if (iact[i] < meq) { // continue; break; } if (sum < 0) { // continue; break; } t1inf = false; it1 = i; } if (!t1inf) { t1 = work[iwuv + it1] / work[iwrv + it1]; for (i = 1; i <= nact; i = i + 1) { if (iact[i] < meq) { // continue; break; } if (work[iwrv + i] < 0) { // continue; break; } temp = work[iwuv + i] / work[iwrv + i]; if (temp < t1) { t1 = temp; it1 = i; } } } sum = 0; for (i = iwzv + 1; i <= iwzv + n; i = i + 1) { sum = sum + work[i] * work[i]; } if (Math.abs(sum) <= vsmall) { if (t1inf) { ierr[1] = 1; // GOTO 999 return 999; } else { for (i = 1; i <= nact; i = i + 1) { work[iwuv + i] = work[iwuv + i] - t1 * work[iwrv + i]; } work[iwuv + nact + 1] = work[iwuv + nact + 1] + t1; // GOTO 700 return 700; } } else { sum = 0; for (i = 1; i <= n; i = i + 1) { sum = sum + work[iwzv + i] * amat[i][nvl]; } tt = -work[iwsv + nvl] / sum; t2min = true; if (!t1inf) { if (t1 < tt) { tt = t1; t2min = false; } } for (i = 1; i <= n; i = i + 1) { sol[i] = sol[i] + tt * work[iwzv + i]; if (Math.abs(sol[i]) < vsmall) { sol[i] = 0; } } crval[1] = crval[1] + tt * sum * (tt / 2 + work[iwuv + nact + 1]); for (i = 1; i <= nact; i = i + 1) { work[iwuv + i] = work[iwuv + i] - tt * work[iwrv + i]; } work[iwuv + nact + 1] = work[iwuv + nact + 1] + tt; if (t2min) { nact = nact + 1; iact[nact] = nvl; l = iwrm + ((nact - 1) * nact) / 2 + 1; for (i = 1; i <= nact - 1; i = i + 1) { work[l] = work[i]; l = l + 1; } if (nact === n) { work[l] = work[n]; } else { for (i = n; i >= nact + 1; i = i - 1) { if (work[i] === 0) { // continue; break; } gc = Math.max(Math.abs(work[i - 1]), Math.abs(work[i])); gs = Math.min(Math.abs(work[i - 1]), Math.abs(work[i])); if (work[i - 1] >= 0) { temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); } else { temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); } gc = work[i - 1] / temp; gs = work[i] / temp; if (gc === 1) { // continue; break; } if (gc === 0) { work[i - 1] = gs * temp; for (j = 1; j <= n; j = j + 1) { temp = dmat[j][i - 1]; dmat[j][i - 1] = dmat[j][i]; dmat[j][i] = temp; } } else { work[i - 1] = temp; nu = gs / (1 + gc); for (j = 1; j <= n; j = j + 1) { temp = gc * dmat[j][i - 1] + gs * dmat[j][i]; dmat[j][i] = nu * (dmat[j][i - 1] + temp) - dmat[j][i]; dmat[j][i - 1] = temp; } } } work[l] = work[nact]; } } else { sum = -bvec[nvl]; for (j = 1; j <= n; j = j + 1) { sum = sum + sol[j] * amat[j][nvl]; } if (nvl > meq) { work[iwsv + nvl] = sum; } else { work[iwsv + nvl] = -Math.abs(sum); if (sum > 0) { for (j = 1; j <= n; j = j + 1) { amat[j][nvl] = -amat[j][nvl]; } bvec[nvl] = -bvec[nvl]; } } // GOTO 700 return 700; } } return 0; } function fn_goto_797() { l = iwrm + (it1 * (it1 + 1)) / 2 + 1; l1 = l + it1; if (work[l1] === 0) { // GOTO 798 return 798; } gc = Math.max(Math.abs(work[l1 - 1]), Math.abs(work[l1])); gs = Math.min(Math.abs(work[l1 - 1]), Math.abs(work[l1])); if (work[l1 - 1] >= 0) { temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); } else { temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); } gc = work[l1 - 1] / temp; gs = work[l1] / temp; if (gc === 1) { // GOTO 798 return 798; } if (gc === 0) { for (i = it1 + 1; i <= nact; i = i + 1) { temp = work[l1 - 1]; work[l1 - 1] = work[l1]; work[l1] = temp; l1 = l1 + i; } for (i = 1; i <= n; i = i + 1) { temp = dmat[i][it1]; dmat[i][it1] = dmat[i][it1 + 1]; dmat[i][it1 + 1] = temp; } } else { nu = gs / (1 + gc); for (i = it1 + 1; i <= nact; i = i + 1) { temp = gc * work[l1 - 1] + gs * work[l1]; work[l1] = nu * (work[l1 - 1] + temp) - work[l1]; work[l1 - 1] = temp; l1 = l1 + i; } for (i = 1; i <= n; i = i + 1) { temp = gc * dmat[i][it1] + gs * dmat[i][it1 + 1]; dmat[i][it1 + 1] = nu * (dmat[i][it1] + temp) - dmat[i][it1 + 1]; dmat[i][it1] = temp; } } return 0; } function fn_goto_798() { l1 = l - it1; for (i = 1; i <= it1; i = i + 1) { work[l1] = work[l]; l = l + 1; l1 = l1 + 1; } work[iwuv + it1] = work[iwuv + it1 + 1]; iact[it1] = iact[it1 + 1]; it1 = it1 + 1; if (it1 < nact) { // GOTO 797 return 797; } return 0; } function fn_goto_799() { work[iwuv + nact] = work[iwuv + nact + 1]; work[iwuv + nact + 1] = 0; iact[nact] = 0; nact = nact - 1; iter[2] = iter[2] + 1; return 0; } go = 0; while (true) { go = fn_goto_50(); if (go === 999) { return; } while (true) { go = fn_goto_55(); if (go === 0) { break; } if (go === 999) { return; } if (go === 700) { if (it1 === nact) { fn_goto_799(); } else { while (true) { fn_goto_797(); go = fn_goto_798(); if (go !== 797) { break; } } fn_goto_799(); } } } } } function solveQP(Dmat, dvec, Amat, bvec, meq, factorized) { Dmat = base0to1(Dmat); dvec = base0to1(dvec); Amat = base0to1(Amat); var i, n, q, nact, r, crval = [], iact = [], sol = [], work = [], iter = [], message; meq = meq || 0; factorized = factorized ? base0to1(factorized) : [undefined, 0]; bvec = bvec ? base0to1(bvec) : []; // In Fortran the array index starts from 1 n = Dmat.length - 1; q = Amat[1].length - 1; if (!bvec) { for (i = 1; i <= q; i = i + 1) { bvec[i] = 0; } } for (i = 1; i <= q; i = i + 1) { iact[i] = 0; } nact = 0; r = Math.min(n, q); for (i = 1; i <= n; i = i + 1) { sol[i] = 0; } crval[1] = 0; for (i = 1; i <= (2 * n + (r * (r + 5)) / 2 + 2 * q + 1); i = i + 1) { work[i] = 0; } for (i = 1; i <= 2; i = i + 1) { iter[i] = 0; } qpgen2(Dmat, dvec, n, n, sol, crval, Amat, bvec, n, q, meq, iact, nact, iter, work, factorized); message = ""; if (factorized[1] === 1) { message = "constraints are inconsistent, no solution!"; } if (factorized[1] === 2) { message = "matrix D in quadratic function is not positive definite!"; } return { solution: base1to0(sol), value: base1to0(crval), unconstrained_solution: base1to0(dvec), iterations: base1to0(iter), iact: base1to0(iact), message: message }; } exports.solveQP = solveQP; }(numeric)); /* Shanti Rao sent me this routine by private email. I had to modify it slightly to work on Arrays instead of using a Matrix object. It is apparently translated from http://stitchpanorama.sourceforge.net/Python/svd.py */ numeric.svd= function svd(A) { var temp; //Compute the thin SVD from G. H. Golub and C. Reinsch, Numer. Math. 14, 403-420 (1970) var prec= numeric.epsilon; //Math.pow(2,-52) // assumes double prec var tolerance= 1.e-64/prec; var itmax= 50; var c=0; var i=0; var j=0; var k=0; var l=0; var u= numeric.clone(A); var m= u.length; var n= u[0].length; if (m < n) throw "Need more rows than columns" var e = new Array(n); var q = new Array(n); for (i=0; i b) return a*Math.sqrt(1.0+(b*b/a/a)) else if (b == 0.0) return a return b*Math.sqrt(1.0+(a*a/b/b)) } //Householder's reduction to bidiagonal form var f= 0.0; var g= 0.0; var h= 0.0; var x= 0.0; var y= 0.0; var z= 0.0; var s= 0.0; for (i=0; i < n; i++) { e[i]= g; s= 0.0; l= i+1; for (j=i; j < m; j++) s += (u[j][i]*u[j][i]); if (s <= tolerance) g= 0.0; else { f= u[i][i]; g= Math.sqrt(s); if (f >= 0.0) g= -g; h= f*g-s u[i][i]=f-g; for (j=l; j < n; j++) { s= 0.0 for (k=i; k < m; k++) s += u[k][i]*u[k][j] f= s/h for (k=i; k < m; k++) u[k][j]+=f*u[k][i] } } q[i]= g s= 0.0 for (j=l; j < n; j++) s= s + u[i][j]*u[i][j] if (s <= tolerance) g= 0.0 else { f= u[i][i+1] g= Math.sqrt(s) if (f >= 0.0) g= -g h= f*g - s u[i][i+1] = f-g; for (j=l; j < n; j++) e[j]= u[i][j]/h for (j=l; j < m; j++) { s=0.0 for (k=l; k < n; k++) s += (u[j][k]*u[i][k]) for (k=l; k < n; k++) u[j][k]+=s*e[k] } } y= Math.abs(q[i])+Math.abs(e[i]) if (y>x) x=y } // accumulation of right hand gtransformations for (i=n-1; i != -1; i+= -1) { if (g != 0.0) { h= g*u[i][i+1] for (j=l; j < n; j++) v[j][i]=u[i][j]/h for (j=l; j < n; j++) { s=0.0 for (k=l; k < n; k++) s += u[i][k]*v[k][j] for (k=l; k < n; k++) v[k][j]+=(s*v[k][i]) } } for (j=l; j < n; j++) { v[i][j] = 0; v[j][i] = 0; } v[i][i] = 1; g= e[i] l= i } // accumulation of left hand transformations for (i=n-1; i != -1; i+= -1) { l= i+1 g= q[i] for (j=l; j < n; j++) u[i][j] = 0; if (g != 0.0) { h= u[i][i]*g for (j=l; j < n; j++) { s=0.0 for (k=l; k < m; k++) s += u[k][i]*u[k][j]; f= s/h for (k=i; k < m; k++) u[k][j]+=f*u[k][i]; } for (j=i; j < m; j++) u[j][i] = u[j][i]/g; } else for (j=i; j < m; j++) u[j][i] = 0; u[i][i] += 1; } // diagonalization of the bidiagonal form prec= prec*x for (k=n-1; k != -1; k+= -1) { for (var iteration=0; iteration < itmax; iteration++) { // test f splitting var test_convergence = false for (l=k; l != -1; l+= -1) { if (Math.abs(e[l]) <= prec) { test_convergence= true break } if (Math.abs(q[l-1]) <= prec) break } if (!test_convergence) { // cancellation of e[l] if l>0 c= 0.0 s= 1.0 var l1= l-1 for (i =l; i= itmax-1) throw 'Error: no convergence.' // shift from bottom 2x2 minor x= q[l] y= q[k-1] g= e[k-1] h= e[k] f= ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y) g= pythag(f,1.0) if (f < 0.0) f= ((x-z)*(x+z)+h*(y/(f-g)-h))/x else f= ((x-z)*(x+z)+h*(y/(f+g)-h))/x // next QR transformation c= 1.0 s= 1.0 for (i=l+1; i< k+1; i++) { g= e[i] y= q[i] h= s*g g= c*g z= pythag(f,h) e[i-1]= z c= f/z s= h/z f= x*c+g*s g= -x*s+g*c h= y*s y= y*c for (j=0; j < n; j++) { x= v[j][i-1] z= v[j][i] v[j][i-1] = x*c+z*s v[j][i] = -x*s+z*c } z= pythag(f,h) q[i-1]= z c= f/z s= h/z f= c*g+s*y x= -s*g+c*y for (j=0; j < m; j++) { y= u[j][i-1] z= u[j][i] u[j][i-1] = y*c+z*s u[j][i] = -y*s+z*c } } e[l]= 0.0 e[k]= f q[k]= x } } //vt= transpose(v) //return (u,q,vt) for (i=0;i= 0; j--) { if (q[j] < q[i]) { // writeln(i,'-',j) c = q[j] q[j] = q[i] q[i] = c for(k=0;k