Skip to content

Instantly share code, notes, and snippets.

@jgbos

jgbos/README.md

Last active December 15, 2015 14:58
Show Gist options
  • Save jgbos/5278157 to your computer and use it in GitHub Desktop.
Save jgbos/5278157 to your computer and use it in GitHub Desktop.
Download ZIP
Wargames
Raw
README.md

Note: Best viewed with Google Chrome; not tested on any other browser. If you view with another browser or on a smartphone the animation may not be smooth due to the computations of trajectories happening as they are drawn.

This demonstrates the power of numerical integration of nonlinear equations using Javascript in the browser. I would have never considered using Javascript for such a task if I had not seen the benchmark tests for Julia and Numeric.js.

The animation is based on an influential movie from my childhood, Wargames. I must have watched it a hundred times as a kid. Since I've been trying to learn D3, I thought I'd try to recreate a portion of the final scene (you can watch here) using javascript. Some thoughts on the scene:

  • Clearly the movie uses a Mercator projection
  • I did not realize until trying to recreate the missile trajectories myself, but the movie actually creates realistic trajectories. Pretty cool for 1983 movie, and amazing the the director let physics be portrayed accurately.
  • Calculations are done using basic rocket equations that assume the earth is spherical and not rotating, this simplifies the math but does not take away from the visualization.

I used the following tools:

My repository contains the earth model and a basic trajectory creation tool for a missile from a launch and impact point (many realistic effects are ignored like atmospheric drag and the need to boost a rocket to the desired velocity). All calculations are done in the browser. The initial point to the trajectories are created on load and the trajectories are flown during the animation.

Sorry, this animation does not support chess or tic-tac-toe :)

Raw
index.html
<!DOCTYPE html>
<meta charset="utf-8">
<style>
svg {
background: #222;
}
.graticule {
fill: none;
stroke: none;
stroke-opacity: .5;
stroke-width: .5px;
}
.land {
fill: none;
stroke: #fff;
stroke-width: 0.5px;
}
.boundary {
fill: none;
stroke: #ccc;
stroke-width: .5px;
}
.bomb {
fill: none;
stroke: #fff;
stroke-width: 2.5px;
stroke-opacity: 0.7;
}
circle.bomb {
fill: #fff;
fill-opacity: 1.0;
stroke: none;
filter: url(#glowimpact);
}
p {
position: absolute;
top: 350px;
color: #fff;
font-family: "Courier New", Courier, monospace;
font-weight: 100;
font-size: 70px;
text-align: center;
width: 960px;
}
</style>
<body>
<p id="text"></p>
<script src="http://d3js.org/d3.v3.min.js"></script>
<script src="http://d3js.org/topojson.v0.min.js"></script>
<script src="http://d3js.org/queue.v1.min.js"></script>
<script src="math.js"></script>
<script src="numeric.js"></script>
<script>
var width = 960,
height = 500;
// Scenario Data
var USSR=643,
USA=840,
russia = [[30,50],[140,70]],
usa = [[-120,25],[-67,50]],
launch_intervals = 1,
reaction_time = 200,
randloc = function (extent){
var lon = extent[0][0] + (extent[1][0]-extent[0][0])*Math.random(),
lat = extent[0][1] + (extent[1][1]-extent[0][1])*Math.random();
return [lat,lon,0];
};
// Functions for orbital calculations
var earth = math.geo.earth("spherical inertial"),
orbital = math.geo.orbital();
var title = d3.select("p");
var projection = d3.geo.mercator()
.scale((width + 1) / 2 / Math.PI)
.translate([width / 2, 400])
.precision(.1);
var path = d3.geo.path()
.projection(projection);
var graticule = d3.geo.graticule();
var svg = d3.select("body").append("svg")
.attr("width", width)
.attr("height", height);
svg.append("filter")
.attr("id", "glow")
.append("feGaussianBlur")
.attr("stdDeviation", 1);
svg.append("filter")
.attr("id", "glowimpact")
.append("feGaussianBlur")
.attr("stdDeviation", 3);
svg.append("path")
.datum(graticule)
.attr("class", "graticule")
.attr("d", path);
queue()
.defer(d3.json,"/d/4090846/world-50m.json")
.defer(scenario,russia,usa)
.defer(scenario,usa, russia)
.await(plot);
function scenario(aggressor, defender, callback){
var launches = [], impacts = [], trajectories=[];
for (var i=0;i<20;i++){
var aloc = randloc(aggressor),
dloc = randloc(defender);
launches.push(aloc);
impacts.push(dloc);
var traj = orbital(aloc,dloc);
trajectories.push({time: Math.random()*launch_intervals, state: traj, history:[], impact: false});
}
callback(null, trajectories)
}
function plot(error, world, rlocs, ulocs){
var land = topojson.object(world, world.objects.land),
countries = topojson.object(world, world.objects.countries).geometries,
borders = topojson.mesh(world, world.objects.countries, function(a, b) { return a.id !== b.id; });
var aggressors = countries.filter(function (a){ return a.id==USSR || a.id==USA;});
svg.selectAll("path.aggressors", ".graticule")
.data(aggressors).enter().append("path")
.attr("class", "aggressors")
.attr("fill","none")
.attr("stroke",function (d) { return d.id==USSR ? "red" : "blue";})
.attr("d", path);
svg.insert("path", ".graticule")
.datum(land)
.attr("class", "land")
.attr("d", path);
svg.insert("path", ".graticule")
.datum(borders)
.attr("class", "boundary")
.attr("d", path);
d3.transition()
.duration(5000)
.each("start", function() {
// title.text("USSR FIRST STRIKE");
})
.tween("rotate", function() {
var time = d3.interpolate(0,8000),
radius = d3.interpolate(1,50);
tlast = 0;
return function(k) {
var tk = time(k);
if ((tk-tlast) < 20) return;
tlast = tk;
var features = [];
var impacts = [];
rlocs.forEach(function (d,i){
propagate(d,tk);
features.push({type: "LineString", coordinates:d.history});
if (d.impact)
impacts.push(d.history[d.history.length-1]);
});
if (tk>reaction_time){
ulocs.forEach(function (d){
propagate(d,tk);
features.push({type: "LineString", coordinates:d.history});
if (d.impact)
impacts.push(d.history[d.history.length-1]);
});
}
var bombs = svg
.selectAll("path.bomb")
.data(features);
bombs.enter().append("path")
.attr("class", "bomb")
.attr("d", path);
bombs
.attr("d", path);
// var craters = svg
// .selectAll("circle.bomb")
// .data(impacts,function (d){return d[0];});
//
// craters.enter().append("circle")
// .attr("class", "bomb")
// .attr("r", 1.0)
// .attr("cx",function (d) { return projection(d)[0];})
// .attr("cy",function (d) {
// return projection(d)[1];
// });
//
// craters
// .attr("r", function (d) { return radius(k);});
};
})
.transition()
.each("end", function (){
// title.text("WINNER: NONE");
});
}
function propagate(d,tk){
if (!d.impact & tk >= d.time){
var sol = numeric.dopri(d.time,tk, d.state,function (t,x){
var g = earth(x);
return [x[3],x[4],x[5],g[0],g[1],g[2]];
},
1e-8,
1000,
function (t,y){
var gc = earth.toGeodetic(y);
return 0 - gc[2];
});
d.time = sol.x[sol.x.length-1];
d.state = sol.y[sol.x.length-1];
var gc = earth.toGeodetic(d.state);
d.impact = gc[2] < 10 && tk>100 ? true : false;
d.history.push([gc[1],gc[0]]);
}
}
</script>
Raw
math.js
math = function() {
var π = Math.PI, ε = 1e-6, math = {
version: "0.0.1"
}, math_radians = π / 180, math_degrees = 180 / π, math_document = document, math_window = window;
math.util = {};
math.util.concat = function(a, b) {
var x = [];
for (var i = 0; i < a.length; i++) x.push(a[i]);
for (var i = 0; i < b.length; i++) x.push(b[i]);
return x;
};
math.psd = {};
math.psd.collapse = function(P) {
var d = numeric.dim(P);
if (d[0] != d[1]) throw Error("P must be square");
var V = [];
for (var i = 0; i < d[0]; i++) {
for (var j = i; j < d[1]; j++) {
V.push(P[i][j]);
}
}
return V;
};
math.psd.dilate = function(V) {
var d = numeric.dim(V);
if (d.length != 1) throw Error("Cannot dilate a matrix");
var n = -.5 + Math.sqrt(.5 * .5 + 2 * d[0]);
var P = [];
k = 0;
for (var i = 0; i < n; i++) P[i] = [];
for (var i = 0; i < n; i++) {
for (var j = i; j < n; j++) {
P[i][j] = V[k];
if (i != j) P[j][i] = V[k];
k++;
}
}
return P;
};
math.quaternion = {};
math.quaternion.q = function q(w, x, y, z) {
this.w = w;
this.x = x;
this.y = y;
this.z = z;
};
math.quaternion.q.prototype.add = function(_) {
return new math.quaternion.q(this.w + _.w, this.x + _.x, this.y + _.y, this.z + _.z);
};
math.quaternion.q.prototype.sub = function(_) {
return new math.quaternion.q(this.w - _.w, this.x - _.x, this.y - _.y, this.z - _.z);
};
math.quaternion.q.prototype.dot = function(_) {
return this.w * _.w + this.x * _.x + this.y * _.y + this.z * _.z;
};
math.quaternion.q.prototype.mul = function(_) {
var a1 = this.w, b1 = this.x, c1 = this.y, d1 = this.z;
if (typeof _ == "object") {
var a2 = _.w, b2 = _.x, c2 = _.y, d2 = _.z;
return new math.quaternion.q(a1 * a2 - b1 * b2 - c1 * c2 - d1 * d2, a1 * b2 + b1 * a2 + c1 * d2 - d1 * c2, a1 * c2 - b1 * d2 + c1 * a2 + d1 * b2, a1 * d2 + b1 * c2 - c1 * b2 + d1 * a2);
} else if (typeof _ == "number") return new math.quaternion.q(a1 * _, b1 * _, c1 * _, d1 * _);
};
math.quaternion.q.prototype.conj = function() {
return new math.quaternion.q(this.w, -this.x, -this.y, -this.z);
};
math.quaternion.q.prototype.norm = function() {
return numeric.norm2([ this.w, this.x, this.y, this.z ]);
};
math.quaternion.q.prototype.yaw = function() {
return Math.atan2(2 * (this.w * this.z + this.x * this.y), 1 - 2 * (this.y * this.y + this.z * this.z));
};
math.quaternion.q.prototype.pitch = function() {
return Math.asin(2 * (this.w * this.y - this.z * this.x));
};
math.quaternion.q.prototype.roll = function() {
return Math.atan2(2 * (this.w * this.x + this.y * this.z), 1 - 2 * (this.x * this.x + this.y * this.y));
};
math.quaternion.q.prototype.dcm = function() {
var w = this.w, x = this.x, y = this.y, z = this.z;
return [ [ w * w + x * x - y * y - z * z, 2 * (x * y - w * z), 2 * (x * z + w * y) ], [ 2 * (y * x + w * z), w * w - x * x + y * y - z * z, 2 * (y * z - w * x) ], [ 2 * (z * x - w * y), 2 * (w * x + y * z), w * w - x * x - y * y + z * z ] ];
};
math.quaternion.q.prototype.vec = function() {
return [ this.w, this.x, this.y, this.z ];
};
math.quaternion.identity = function() {
return new math.quaternion.q(1, 0, 0, 0);
};
math.quaternion.rotateX = function(theta) {
return new math.quaternion.q(Math.cos(theta / 2), Math.sin(theta / 2), 0, 0);
};
math.quaternion.rotateY = function(theta) {
return new math.quaternion.q(Math.cos(theta / 2), 0, Math.sin(theta / 2), 0);
};
math.quaternion.rotateZ = function(theta) {
return new math.quaternion.q(Math.cos(theta / 2), 0, 0, Math.sin(theta / 2));
};
math.quaternion.rotate = function(v, theta) {
if (v.length != 3) throw new Error("math.quaternion.rotate defined for R3");
var n = Math.sin(theta / 2) / numeric.norm2(v);
return new math.quaternion.q(Math.cos(theta / 2), v[0] * n, v[1] * n, v[2] * n);
};
math.quaternion.euler = function(yaw, pitch, roll) {
temp1 = math.quaternion.rotateZ(yaw);
temp2 = temp1.mul(math.quaternion.rotateY(pitch));
return temp2.mul(math.quaternion.rotateX(roll));
};
math.quaternion.permuteZXY = function() {
return new math.quaternion.q(.5, .5, .5, .5);
};
math.quaternion.permuteYZX = function() {
return new math.quaternion.q(.5, -.5, -.5, -.5);
};
math.quaternion.swapXYflipZ = function() {
return new math.quaternion.q(0, Math.sqrt(2) / 2, Math.sqrt(2) / 2, 0);
};
math.quaternion.swapXZflipY = function() {
return new math.quaternion.q(0, Math.sqrt(2) / 2, 0, Math.sqrt(2) / 2);
};
math.quaternion.swapYZflipX = function() {
return new math.quaternion.q(0, 0, Math.sqrt(2) / 2, Math.sqrt(2) / 2);
};
math.quaternion.slerp = function(q1, q2, u) {
var theta = Math.acos(q1.dot(q2));
var st = Math.sin(theta);
var su = Math.sin(u * theta);
var suu = Math.sin((1 - u) * theta);
var a1 = q1.mul(suu / st);
var a2 = q2.mul(su / st);
var b = a1.add(a2);
return b.mul(1 / b.norm());
};
math.models = {};
math.models.poly = function(order, n) {
var J, ns = (order + 1) * n;
function derivative(t, x) {
calcJ(t, x);
return numeric.dot(J, x);
}
function calcJ(t, x) {
if (x.length != ns) throw Error("Wrong number of states");
J = numeric.rep([ ns, ns ], 0);
var I = numeric.identity(ns - 1);
numeric.setBlock(J, [ 0, n ], [ ns - n - 1, ns - 1 ], I);
}
derivative.jacobian = function(t, x) {
if (arguments.length == 2) calcJ(t, x);
return J;
};
derivative.numStates = function(_) {
if (!arguments.length) return ns;
ns = _;
return derivative;
};
return derivative;
};
math.models.meas = {};
math.models.meas.simple = function(m, ns) {
var J;
function to(t, x) {
calcJ(t, x);
return numeric.dot(J, x);
}
function calcJ(t, x) {
if (!ns) ns = x.length;
J = numeric.rep([ m, ns ], 0);
numeric.setBlock(J, [ 0, 0 ], [ m - 1, m - 1 ], numeric.identity(m));
}
to.jacobian = function(t, x) {
if (arguments.length == 2) calcJ(t, x);
return J;
};
return to;
};
math.state = function(t, x, P) {
var extra, flattened;
s = {};
s.t = t;
s.x = x;
s.P = P;
s.flatten = function() {
var xout = [];
for (var i in x) xout.push(x[i]);
var pvec = math.psd.collapse(P);
for (var i in pvec) xout.push(pvec[i]);
s.flattened = xout;
return xout;
};
s.meta = function(_) {
if (!arguments.length) return extra;
extra = _;
return s;
};
return s;
};
math.state.dilate = function(t, vec) {
if (arguments.length == 1) {
vec = t;
t = 0;
}
var n = 1;
switch (vec.length) {
case 2:
n = 1;
break;
case 5:
n = 2;
break;
case 9:
n = 3;
break;
}
var x = numeric.getBlock(vec, [ 0 ], [ n - 1 ]);
var P = math.psd.dilate(numeric.getBlock(vec, [ n ], [ vec.length - 1 ]));
return new math.state(t, x, P);
};
math.kalman = {};
math.kalman.dynamics = function(f, q) {
return function(t, xP) {
var state = math.state.dilate(t, xP);
var x = state.x;
var P = state.P;
var xd = f(t, x);
var F = f.jacobian();
var Ft = numeric.transpose(F);
var Pd = numeric.add(numeric.dot(F, P), numeric.dot(P, Ft));
if (q) var Pd = numeric.add(Pd, q(t, x));
var p = math.psd.collapse(Pd);
for (var i = 0; i < p.length; i++) xd.push(p[i]);
return xd;
};
};
math.kalman.predict = function(f, q) {
return function(state, time_to) {
if (state.t == time_to.t) return state;
var vec = state.flatten(), fdyn = math.kalman.dynamics(f, q);
var sol = numeric.dopri(state.t, time_to, vec, fdyn);
var n = sol.x.length - 1;
return math.state.dilate(time_to, sol.y[n]);
};
};
math.kalman.update = function(h) {
return function(state, meas) {
var t = state.t, x = state.x, P = state.P, z = meas.z, R = meas.R;
if (!R.length) R = [ [ R ] ];
var zp = h(t, x), H = h.jacobian(), Ht = numeric.transpose(H), Rp = numeric.dot(H, numeric.dot(P, Ht));
var S = numeric.add(Rp, R);
var K = numeric.dot(numeric.dot(P, Ht), numeric.inv(S));
var innovation = numeric.sub(z, zp);
var xE = numeric.add(x, numeric.dot(K, innovation));
var I = numeric.identity(x.length), kh = numeric.dot(K, H), ikh = numeric.sub(I, kh), ikht = numeric.transpose(ikh), krk = numeric.dot(K, numeric.dot(R, numeric.transpose(K)));
var Pe = numeric.add(numeric.dot(ikh, numeric.dot(P, ikht)), krk);
var s = new math.state(t, xE, Pe);
return s.meta({
Gain: K,
innovation: innovation,
S: S
});
};
};
math.geo = {};
math.geo.earth = math_geo_earth;
function math_geo_earth(earth_type) {
var Re = 6378137, Rp = 6356752.3142, μ = 3986004418e5, μA = 35e7, ω = 72921151467e-15, f = 1 / 298.257223563, ε = .08181919084262149, J2 = .0010826299890519;
if (arguments.length == 1) {
type = earth_type.split(" ");
type.forEach(function(t) {
switch (t) {
case "wgs84":
break;
case "spherical":
Re = 6371e3;
f = 0;
ε = 0;
J2 = 0;
break;
case "inertial":
ω = 0;
break;
}
});
}
var model = function(state) {
var Ω = [ [ 0, -ω, 0 ], [ ω, 0, 0 ], [ 0, 0, 0 ] ];
var ΩΩ = [ [ -ω * ω, 0, 0 ], [ 0, -ω * ω, 0 ], [ 0, 0, 0 ] ];
var r = state.slice(0, 3);
var v = state.slice(3, 6);
var R = Math.sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]), RSq = R * R, ReSq = Re * Re, ur = [ r[0] / R, r[1] / R, r[2] / R ], up = [ 0, 0, 1 ], Cc = r[2] / R, CcSq = Cc * Cc;
var g = [ 0, 0, 0 ];
for (var i = 0; i < 3; i++) {
g[i] = -μ / RSq * (ur[i] - 1.5 * (J2 * ReSq / RSq) * ((5 * CcSq - 1) * ur[i] - 2 * Cc * up[i]));
}
var a = [ 0, 0, 0 ];
var coriolis = [ 0, 0, 0 ];
var centrifugal = [ 0, 0, 0 ];
for (var i = 0; i < 3; i++) {
if (v.length == 3) {
for (var j = 0; j < 3; j++) {
coriolis[i] += -2 * Ω[i][j] * v[j];
centrifugal[i] += -ΩΩ[i][j] * r[j];
}
}
a[i] = g[i] + coriolis[i] + centrifugal[i];
}
return a;
};
model.toGeodetic = function(position) {
var x = position[0], y = position[1], z = position[2], xSq = x * x, ySq = y * y;
var eSq = ε * ε, r = Math.sqrt(xSq + ySq), b = Re * Math.sqrt(1 - eSq);
var lon = Math.atan2(y, x);
if (ε == 0) {
lat = Math.atan2(z, r);
H = Math.sqrt(x * x + y * y + z * z) - Re;
} else {
var l = eSq / 2, lSq = l * l, m = Math.pow(r / Re, 2), n = Math.pow((1 - eSq) * z / b, 2), i = -(2 * lSq + m + n) / 2, k = lSq * (lSq - m - n), q = Math.pow(m + n - 4 * lSq, 3) / 216 + m * n * lSq, D = Math.sqrt((2 * q - m * n * lSq) * m * n * lSq), beta = i / 3 - Math.pow(q + D, 1 / 3) - Math.pow(q - D, 1 / 3);
var sign = 0;
if (m - n > 0) {
sign = 1;
} else if (m - n < 0) {
sign = -1;
}
var t = Math.sqrt(Math.sqrt(beta * beta - k) - (beta + i) / 2) - sign * Math.sqrt((beta - i) / 2);
var r0 = r / (t + l);
var z0 = (1 - eSq) * z / (t - l);
var lat = Math.atan(z0 / ((1 - eSq) * r0));
sign = 0;
if (t - 1 + l > 0) {
sign = 1;
} else if (t - 1 + l < 0) {
sign = -1;
}
var H = sign * Math.sqrt((r - r0) * (r - r0) + (z - z0) * (z - z0));
}
return [ lat * math_degrees, lon * math_degrees, H ];
};
model.toCartesian = function(gc) {
var lat = gc[0] * math_radians, lon = gc[1] * math_radians, alt = gc[2], slat = Math.sin(lat), clat = Math.cos(lat);
var R = Re / Math.sqrt(1 - ε * ε * slat * slat);
var x = (R + alt) * clat * Math.cos(lon);
var y = (R + alt) * clat * Math.sin(lon);
var z = (R + alt - ε * ε * R) * slat;
return [ x, y, z ];
};
model.semimajor_axis = Re;
model.semiminor_axis = Rp;
model.eccentricity = ε;
model.flattening = f;
model.rotation_rate = ω;
model.GM = μ;
model.GM_atm = μA;
model.J2 = J2;
return model;
}
math.geo.orbital = math_geo_orbital;
function math_geo_orbital() {
var earth = math.geo.earth("spherical inertial"), Re = earth.semimajor_axis, μ = earth.GM;
var γ = 45 * math_radians;
function model(launch, impact) {
var R0 = earth.toCartesian(launch), r0 = numeric.norm2(R0), ur0 = numeric.mul(R0, 1 / r0);
var Rf = earth.toCartesian(impact), rf = numeric.norm2(Rf), urf = numeric.mul(Rf, 1 / rf);
var φ = Math.acos(numeric.dot(ur0, urf)), V = Math.sqrt(μ * (1 - Math.cos(φ)) / (r0 * Math.cos(γ) * (r0 * Math.cos(γ) / Re - Math.cos(φ + γ))));
var tempK = cross(ur0, urf), K = numeric.mul(tempK, 1 / numeric.norm2(tempK));
var q = math.quaternion.rotate(K, Math.PI / 2 - γ), qv = new math.quaternion.q(0, ur0[0], ur0[1], ur0[2]), rotR = q.mul(qv.mul(q.conj())).vec().slice(1, 4), urotR = numeric.mul(rotR, 1 / numeric.norm2(rotR)), velocity = numeric.mul(urotR, V);
return [ 0, 0, 0, 0, 0, 0 ].map(function(d, i) {
if (i < 3) return R0[i]; else return velocity[i - 3];
});
}
model.flightTime = function(r, φ, V) {
var cg = Math.cos(γ), sg = Math.sin(γ), tg = Math.tan(γ), cp = Math.cos(φ), sp = Math.sin(φ), cgp = Math.cos(γ + φ), cot = 1 / Math.tan(φ / 2), λ = r * V * V / μ;
var first_term = (tg * (1 - cp) + (1 - λ) * sp) / ((2 - λ) * ((1 - cp) / (λ * cg * cg) + cgp / cg)), sec_term = 2 * cg / (λ * Math.pow(2 / λ - 1, 1.5)) * Math.atan2(Math.sqrt(2 / λ - 1), cg * cot - sg);
return r / (V * cg) * (first_term + sec_term);
};
function cross(a, b) {
return [ b[2] * a[1] - b[1] * a[2], b[0] * a[2] - b[2] * a[0], b[1] * a[0] - b[0] * a[1] ];
}
model.velocity = function(_) {
if (!arguments.length) return V;
V = _;
return model;
};
model.angle = function(_) {
if (!arguments.length) return γ;
γ = _;
return model;
};
return model;
}
math.radar = {};
math.radar.header = math_radar_header;
function math_radar_header() {
var time = 0, nf = -30, f0 = 2424e5, bw = 1335e5, location = [ 0, 0 ], numberSamples = Math.pow(2, 15), numberChannels = 1269;
var header = {};
header.origin = function() {
var elements = [], l = 3e8 / f0;
for (var i = 0; i < numberChannels; i++) elements.push({
x: (i - (numberChannels - 1) / 2) * (l / 2) + location[0],
y: location[1]
});
return elements;
};
header.time = function(_) {
if (!arguments.length) return time;
time = _;
return header;
};
header.f0 = function(_) {
if (!arguments.length) return f0;
f0 = _;
return header;
};
header.bw = function(_) {
if (!arguments.length) return bw;
bw = _;
return header;
};
header.nf = function(_) {
if (!arguments.length) return nf;
nf = _;
return header;
};
header.loc = function(_) {
if (!arguments.length) return location;
location = _;
return header;
};
header.numberSamples = function(_) {
if (!arguments.length) return numberSamples;
numberSamples = _;
return header;
};
return header;
}
math.radar.pb = math_radar_pulse_builder;
function math_radar_pulse_builder() {
var pulse, header, inFrequencyDomain = true;
function pb() {
var z = numeric.rep([ header.numberSamples() ], 0);
pulse = numeric.t(z, z);
}
pb.toTimeDomain = function() {
if (!inFrequencyDomain) return pulse;
inFrequencyDomain = false;
return pulse.ifft();
};
pb.addNoise = function(noiseFloor) {
var sigf = Math.sqrt(.5 * pulse.x.length) * Math.pow(10, noiseFloor / 20);
pulse.x = pulse.x.map(function(d) {
return d + jStat.randn() * sigf;
});
pulse.y = pulse.y.map(function(d) {
return d + jStat.randn() * sigf;
});
};
pb.addPointResponse = function(state, power) {
var f = numeric.linspace(header.f0() - header.bw() / 2, header.f0() + header.bw() / 2, pulse.x.length);
var k = numeric.mul(4 * Math.PI / 3e8, f);
var loc = header.loc();
var r = numeric.norm2(numeric.sub(state, loc));
var A = Math.pow(10, power / 20);
for (var i = 0; i < pulse.x.length; i++) {
pulse.x[i] += A * Math.cos(-k[i] * r);
pulse.y[i] += A * Math.sin(-k[i] * r);
}
};
pb.pulse = function(_) {
if (!arguments.length) return pulse;
pulse = _;
return pb;
};
pb.header = function(_) {
if (!arguments.length) return header;
header = _;
return pb;
};
pb.header = function(_) {
if (!arguments.length) return header;
header = _;
return pb;
};
return pb;
}
math.plot = {};
math.plot.axis = math_plot_axis;
function math_plot_axis() {
var axis = d3.svg.axis(), scale = d3.scale.linear(), ticks = null, grid, gridStroke = function(d) {
return d ? "#ccc" : "#666";
}, interp, label, labelMargin = 60, container;
axis.scale(scale).orient("bottom").tickFormat(function(d) {
return d;
});
function chart(selection) {
selection.each(function(data) {
container = d3.select(this);
var g = container.selectAll("g.math.axis").data([ data ]);
var g1 = g.enter().append("g").attr("class", "math axis");
g1.append("g").attr("class", "the axis");
if (ticks !== null) axis.ticks(ticks);
chart.draw();
});
}
function draw_label() {
var g = container.select("g.the.axis").selectAll("text.label").data([ 0 ]);
var scale = axis.scale();
switch (axis.orient()) {
case "bottom":
var w = scale.range()[1];
g.enter().append("text").attr("class", "label").attr("x", w / 2).attr("y", labelMargin).attr("dx", "1em").style("text-anchor", "middle").text(label);
break;
case "left":
g.enter().append("text").attr("class", "label").attr("transform", "rotate(-90)").attr("y", 0 - labelMargin).attr("x", 0 - scale.range()[0] / 2).attr("dy", "1em").style("text-anchor", "middle").text(label);
break;
}
}
chart.draw = function() {
var g = container.select("g.the.axis");
if (interp) {
g.transition().ease("linear").call(axis);
} else g.call(axis);
if (grid) {
g.selectAll("g").each(function() {
d3.select(this).append("svg:line").attr("class", "axis line").style("stroke", gridStroke).attr("x1", grid.x1).attr("x2", grid.x2).attr("y1", grid.y1).attr("y2", grid.y2);
});
}
draw_label();
};
chart.axis = axis;
d3.rebind(chart, axis, "orient", "tickValues", "tickSubdivide", "tickSize", "tickPadding", "tickFormat");
chart.margin = function(_) {
if (!arguments.length) return labelMargin;
labelMargin = _;
return chart;
};
chart.scale = function(_) {
if (!arguments.length) return scale;
scale = _;
axis.scale(scale);
return chart;
};
chart.ticks = function(_) {
if (!arguments.length) return ticks;
ticks = _;
return chart;
};
chart.grid = function(_) {
if (!arguments.length) return grid;
grid = _;
return chart;
};
chart.label = function(_) {
if (!arguments.length) return label;
label = _;
return chart;
};
chart.interp = function(_) {
if (!arguments.length) return interp;
interp = _;
return chart;
};
return chart;
}
math.plot.zoom = math_plot_zoom;
function math_plot_zoom() {
var container;
function chart(selection, f) {
selection.each(function() {
var xAxis = f.xAxis();
var yAxis = f.yAxis();
var width = f.width();
var height = f.height();
var margin = f.margin();
var zoomable = f.zoom();
var xyzoom = d3.behavior.zoom().x(xAxis.scale()).y(yAxis.scale()).on("zoom", zoomable ? f.draw : null);
var xzoom = d3.behavior.zoom().x(xAxis.scale()).on("zoom", zoomable ? f.draw : null);
var yzoom = d3.behavior.zoom().y(yAxis.scale()).on("zoom", zoomable ? f.draw : null);
container = d3.select(this).append("g").attr("class", "zoom").attr("transform", "translate(" + margin.left + "," + margin.top + ")");
container.append("svg:rect").attr("class", "zoom xy box").attr("width", width - margin.left - margin.right).attr("height", height - margin.top - margin.bottom).style("visibility", "hidden").attr("pointer-events", "all").call(xyzoom);
container.append("svg:rect").attr("class", "zoom x box").attr("width", width).attr("height", margin.bottom).attr("transform", "translate(" + 0 + "," + (height - margin.top - margin.bottom) + ")").style("visibility", "hidden").attr("pointer-events", "all").call(xzoom);
container.append("svg:rect").attr("class", "zoom y box").attr("width", margin.left).attr("height", height).attr("transform", "translate(" + -margin.left + "," + 0 + ")").style("visibility", "hidden").attr("pointer-events", "all").call(yzoom);
});
}
chart.update = function(f) {
var xAxis = f.xAxis();
var yAxis = f.yAxis();
var zoomable = f.zoom();
var xyzoom = d3.behavior.zoom().x(xAxis.scale()).y(yAxis.scale()).on("zoom", zoomable ? f.draw : null);
var xzoom = d3.behavior.zoom().x(xAxis.scale()).on("zoom", zoomable ? f.draw : null);
var yzoom = d3.behavior.zoom().y(yAxis.scale()).on("zoom", zoomable ? f.draw : null);
container.select("rect.zoom.xy.box").call(xyzoom);
container.select("rect.zoom.x.box").call(xzoom);
container.select("rect.zoom.y.box").call(yzoom);
};
return chart;
}
math.plot.scatter = math_plot_scatter;
function math_plot_scatter() {
var margin = {
top: 60,
bottom: 80,
left: 60,
right: 0
}, width = 940, height = 600, xValue = function(d) {
return d.x;
}, yValue = function(d) {
return d.y;
}, cValue = function() {
return "steelblue";
}, fValue = null, xAxis = math.plot.axis(), yAxis = math.plot.axis(), cScale = d3.scale.linear().domain([ 0, 1 ]).range([ "hsl(250, 100%, 50%)", "hsl(0, 100%, 50%)" ]).interpolate(d3.interpolateHsl), markerSize = function() {
return 1;
}, dfilter = function() {
return true;
}, zoom = math.plot.zoom(), zoomable = false, container;
function chart(selection) {
selection.each(function(data) {
data = data.map(function(d, i) {
if (!fValue) return [ xValue.call(data, d, i), yValue.call(data, d, i), cValue.call(data, d, i) ]; else return [ xValue.call(data, d, i), yValue.call(data, d, i), cValue.call(data, d, i), fValue.call(data, d, i) ];
});
container = d3.select(this).selectAll("svg").data([ data ]);
container.enter().append("svg");
var g = container.append("g").attr("transform", "translate(" + margin.left + "," + margin.top + ")");
g.append("defs").append("clipPath").attr("id", "clip").append("rect").attr("width", width - margin.left - margin.right).attr("height", height - margin.top - margin.bottom);
g.append("svg:rect").attr("class", "border").attr("width", width - margin.left - margin.right).attr("height", height - margin.top - margin.bottom).style("stroke", "black").style("fill", "none");
g.append("g").attr("class", "x axis").attr("transform", "translate(" + 0 + "," + (height - margin.top - margin.bottom) + ")");
g.append("g").attr("class", "y axis");
g.append("g").attr("class", "scatter").attr("clip-path", "url(#clip)");
g.append("g").attr("class", "legend");
if (cValue) cScale.domain(d3.extent(data, function(d) {
return d[2];
}));
xAxis.grid({
x1: 0,
x2: 0,
y1: -(height - margin.top - margin.bottom),
y2: 0
}).orient("bottom").tickPadding(10);
yAxis.grid({
x1: width - margin.left - margin.right,
x2: 0,
y1: 0,
y2: 0
}).orient("left").tickPadding(10);
xAxis.scale().domain(d3.extent(data, function(d) {
return d[0];
})).range([ 0, width - margin.left - margin.right ]);
container.select("g.x.axis").call(xAxis);
yAxis.scale().domain(d3.extent(data, function(d) {
return d[1];
})).range([ height - margin.top - margin.bottom, 0 ]);
container.select("g.y.axis").call(yAxis);
container.call(zoom, chart);
chart.draw();
});
return chart;
}
function update() {
var gs = container.select("g.scatter");
var circle = gs.selectAll("circle").data(function(d) {
return d;
});
circle.style("fill", function(d) {
return dfilter(d) ? C(d) : "none";
}).attr("cx", function(d) {
return X(d);
}).attr("cy", function(d) {
return Y(d);
}).attr("r", markerSize);
circle.enter().append("svg:circle").attr("class", "points").style("fill", function(d) {
return C(d);
}).attr("cx", function(d) {
return X(d);
}).attr("cy", function(d) {
return Y(d);
}).attr("r", markerSize);
circle.exit().remove();
}
chart.draw = function() {
container.select("g.x.axis").call(xAxis);
container.select("g.y.axis").call(yAxis);
update();
zoom.update(chart);
};
function X(d) {
var xScale = xAxis.scale();
return xScale(d[0]);
}
function Y(d) {
var yScale = yAxis.scale();
return yScale(d[1]);
}
function C(d) {
return cScale(d[2]);
}
chart.xaxis = xAxis;
d3.rebind(chart, xAxis, "orient", "tickValues", "tickSubdivide", "tickSize", "tickPadding", "tickFormat");
chart.yaxis = yAxis;
d3.rebind(chart, yAxis, "orient", "tickValues", "tickSubdivide", "tickSize", "tickPadding", "tickFormat");
chart.margin = function(_) {
if (!arguments.length) return margin;
margin.top = typeof _.top != "undefined" ? _.top : margin.top;
margin.right = typeof _.right != "undefined" ? _.right : margin.right;
margin.bottom = typeof _.bottom != "undefined" ? _.bottom : margin.bottom;
margin.left = typeof _.left != "undefined" ? _.left : margin.left;
return chart;
};
chart.width = function(_) {
if (!arguments.length) return width;
width = _;
return chart;
};
chart.height = function(_) {
if (!arguments.length) return height;
height = _;
return chart;
};
chart.markersize = function(_) {
if (!arguments.length) return markerSize;
markerSize = _;
return chart;
};
chart.xAxis = function(_) {
if (!arguments.length) return xAxis;
xAxis = _;
return chart;
};
chart.yAxis = function(_) {
if (!arguments.length) return yAxis;
yAxis = _;
return chart;
};
chart.caxis = function(_) {
if (!arguments.length) return cScale.domain();
cScale.domain(_);
return chart;
};
chart.x = function(_) {
if (!arguments.length) return xValue;
xValue = _;
return chart;
};
chart.y = function(_) {
if (!arguments.length) return yValue;
yValue = _;
return chart;
};
chart.c = function(_) {
if (!arguments.length) return cValue;
cValue = _;
return chart;
};
chart.filter = function(_, __) {
if (!arguments.length) return [ fValue, dfilter ];
fValue = _;
dfilter = __;
return chart;
};
chart.zoom = function(_) {
if (!arguments.length) return zoomable;
zoomable = _;
return chart;
};
return chart;
}
math.plot.line = math_plot_line;
function math_plot_line() {
var margin = {
top: 60,
bottom: 80,
left: 60,
right: 0
}, width = 940, height = 600, zoom = math.plot.zoom(), zoomable = false, xValue = function(d) {
return d.x;
}, yValue = function(d) {
return d.y;
}, groupValue = function(d) {
return d.id;
}, xAxis = math.plot.axis(), yAxis = math.plot.axis(), color = d3.scale.category20c(), lineWidth = function() {
return "2pt";
}, interp = false, ylim, container;
function chart(selection) {
selection.each(function(data) {
var flattened_data = data.map(function(d, i) {
return [ xValue.call(data, d, i), yValue.call(data, d, i), groupValue.call(data, d, i) ];
});
data = d3.nest().key(function(d) {
return I(d);
}).entries(flattened_data);
container = d3.select(this).selectAll("svg").data([ data ]);
container.enter().append("svg");
var g = container.append("g").attr("transform", "translate(" + margin.left + "," + margin.top + ")");
g.append("defs").append("clipPath").attr("id", "clip").append("rect").attr("width", width - margin.left - margin.right).attr("height", height - margin.top - margin.bottom);
g.append("svg:rect").attr("class", "border").attr("width", width - margin.left - margin.right).attr("height", height - margin.top - margin.bottom).style("stroke", "black").style("fill", "none");
g.append("g").attr("class", "x axis").attr("transform", "translate(" + 0 + "," + (height - margin.top - margin.bottom) + ")");
g.append("g").attr("class", "y axis");
g.append("g").attr("class", "lineChart").attr("clip-path", "url(#clip)");
g.append("g").attr("class", "legend");
xAxis.grid({
x1: 0,
x2: 0,
y1: -(height - margin.top - margin.bottom),
y2: 0
}).orient("bottom").tickPadding(10);
xAxis.scale().domain(d3.extent(flattened_data, function(d) {
return d[0];
})).range([ 0, width - margin.left - margin.right ]);
container.select("g.x.axis").call(xAxis);
yAxis.grid({
x1: width - margin.left - margin.right,
x2: 0,
y1: 0,
y2: 0
}).orient("left").tickPadding(10);
if (ylim) {
yAxis.scale().domain(ylim).range([ height - margin.top - margin.bottom, 0 ]);
} else {
yAxis.scale().domain(d3.extent(flattened_data, function(d) {
return d[1];
})).range([ height - margin.top - margin.bottom, 0 ]);
}
container.select("g.y.axis").call(yAxis);
container.call(zoom, chart);
chart.draw();
});
}
chart.update = function() {
var line = d3.svg.line().x(function(d) {
return X(d);
}).y(function(d) {
return Y(d);
}).interpolate("linear");
var gl = container.select("g.lineChart");
var l = gl.selectAll("path.lines").data(function(d) {
return d;
});
if (interp) {
l.transition().ease("linear");
}
l.attr("class", "lines").style("stroke", function(d) {
return color(+d.key);
}).style("stroke-width", lineWidth).style("fill", "none").attr("d", function(d) {
return line(d.values);
});
l.enter().append("svg:path").attr("class", "lines").style("stroke", function(d) {
return color(+d.key);
}).style("stroke-width", lineWidth).style("fill", "none").attr("d", function(d) {
return line(d.values);
});
if (interp) {
l.exit().transition().ease("linear").remove();
} else {
l.exit().remove();
}
};
chart.draw = function() {
container.select("g.x.axis").call(xAxis);
container.select("g.y.axis").call(yAxis);
chart.update();
zoom.update(chart);
};
function X(d) {
var xScale = xAxis.scale();
return xScale(d[0]);
}
function Y(d) {
var yScale = yAxis.scale();
return yScale(d[1]);
}
function I(d) {
return d[2];
}
chart.xaxis = xAxis;
d3.rebind(chart, xAxis, "orient", "tickValues", "tickSubdivide", "tickSize", "tickPadding", "tickFormat");
chart.yaxis = yAxis;
d3.rebind(chart, yAxis, "orient", "tickValues", "tickSubdivide", "tickSize", "tickPadding", "tickFormat");
chart.margin = function(_) {
if (!arguments.length) return margin;
margin.top = typeof _.top != "undefined" ? _.top : margin.top;
margin.right = typeof _.right != "undefined" ? _.right : margin.right;
margin.bottom = typeof _.bottom != "undefined" ? _.bottom : margin.bottom;
margin.left = typeof _.left != "undefined" ? _.left : margin.left;
return chart;
};
chart.width = function(_) {
if (!arguments.length) return width;
width = _;
return chart;
};
chart.height = function(_) {
if (!arguments.length) return height;
height = _;
return chart;
};
chart.xAxis = function(_) {
if (!arguments.length) return xAxis;
xAxis = _;
return chart;
};
chart.yAxis = function(_) {
if (!arguments.length) return yAxis;
yAxis = _;
return chart;
};
chart.x = function(_) {
if (!arguments.length) return xValue;
xValue = _;
return chart;
};
chart.y = function(_) {
if (!arguments.length) return yValue;
yValue = _;
return chart;
};
chart.group = function(_) {
if (!arguments.length) return groupValue;
groupValue = _;
return chart;
};
chart.linewidth = function(_) {
if (!arguments.length) return lineWidth;
lineWidth = _;
return chart;
};
chart.colorscale = function(_) {
if (!arguments.length) return color;
color = _;
return chart;
};
chart.zoom = function(_) {
if (!arguments.length) return zoomable;
zoomable = _;
return chart;
};
chart.interp = function(_) {
if (!arguments.length) return interp;
interp = _;
return chart;
};
chart.ylim = function(_) {
if (!arguments.length) return ylim;
ylim = _;
return chart;
};
return chart;
}
return math;
}();
Raw
numeric.js
"use strict";
var numeric = (typeof exports === "undefined")?(function numeric() {}):(exports);
if(typeof global !== "undefined") { global.numeric = numeric; }
numeric.version = "1.2.6";
// 1. Utility functions
numeric.bench = function bench (f,interval) {
var t1,t2,n,i;
if(typeof interval === "undefined") { interval = 15; }
n = 0.5;
t1 = new Date();
while(1) {
n*=2;
for(i=n;i>3;i-=4) { f(); f(); f(); f(); }
while(i>0) { f(); i--; }
t2 = new Date();
if(t2-t1 > interval) break;
}
for(i=n;i>3;i-=4) { f(); f(); f(); f(); }
while(i>0) { f(); i--; }
t2 = new Date();
return 1000*(3*n-1)/(t2-t1);
}
numeric._myIndexOf = (function _myIndexOf(w) {
var n = this.length,k;
for(k=0;k<n;++k) if(this[k]===w) return k;
return -1;
});
numeric.myIndexOf = (Array.prototype.indexOf)?Array.prototype.indexOf:numeric._myIndexOf;
numeric.Function = Function;
numeric.precision = 4;
numeric.largeArray = 50;
numeric.prettyPrint = function prettyPrint(x) {
function fmtnum(x) {
if(x === 0) { return '0'; }
if(isNaN(x)) { return 'NaN'; }
if(x<0) { return '-'+fmtnum(-x); }
if(isFinite(x)) {
var scale = Math.floor(Math.log(x) / Math.log(10));
var normalized = x / Math.pow(10,scale);
var basic = normalized.toPrecision(numeric.precision);
if(parseFloat(basic) === 10) { scale++; normalized = 1; basic = normalized.toPrecision(numeric.precision); }
return parseFloat(basic).toString()+'e'+scale.toString();
}
return 'Infinity';
}
var ret = [];
function foo(x) {
var k;
if(typeof x === "undefined") { ret.push(Array(numeric.precision+8).join(' ')); return false; }
if(typeof x === "string") { ret.push('"'+x+'"'); return false; }
if(typeof x === "boolean") { ret.push(x.toString()); return false; }
if(typeof x === "number") {
var a = fmtnum(x);
var b = x.toPrecision(numeric.precision);
var c = parseFloat(x.toString()).toString();
var d = [a,b,c,parseFloat(b).toString(),parseFloat(c).toString()];
for(k=1;k<d.length;k++) { if(d[k].length < a.length) a = d[k]; }
ret.push(Array(numeric.precision+8-a.length).join(' ')+a);
return false;
}
if(x === null) { ret.push("null"); return false; }
if(typeof x === "function") {
ret.push(x.toString());
var flag = false;
for(k in x) { if(x.hasOwnProperty(k)) {
if(flag) ret.push(',\n');
else ret.push('\n{');
flag = true;
ret.push(k);
ret.push(': \n');
foo(x[k]);
} }
if(flag) ret.push('}\n');
return true;
}
if(x instanceof Array) {
if(x.length > numeric.largeArray) { ret.push('...Large Array...'); return true; }
var flag = false;
ret.push('[');
for(k=0;k<x.length;k++) { if(k>0) { ret.push(','); if(flag) ret.push('\n '); } flag = foo(x[k]); }
ret.push(']');
return true;
}
ret.push('{');
var flag = false;
for(k in x) { if(x.hasOwnProperty(k)) { if(flag) ret.push(',\n'); flag = true; ret.push(k); ret.push(': \n'); foo(x[k]); } }
ret.push('}');
return true;
}
foo(x);
return ret.join('');
}
numeric.parseDate = function parseDate(d) {
function foo(d) {
if(typeof d === 'string') { return Date.parse(d.replace(/-/g,'/')); }
if(!(d instanceof Array)) { throw new Error("parseDate: parameter must be arrays of strings"); }
var ret = [],k;
for(k=0;k<d.length;k++) { ret[k] = foo(d[k]); }
return ret;
}
return foo(d);
}
numeric.parseFloat = function parseFloat_(d) {
function foo(d) {
if(typeof d === 'string') { return parseFloat(d); }
if(!(d instanceof Array)) { throw new Error("parseFloat: parameter must be arrays of strings"); }
var ret = [],k;
for(k=0;k<d.length;k++) { ret[k] = foo(d[k]); }
return ret;
}
return foo(d);
}
numeric.parseCSV = function parseCSV(t) {
var foo = t.split('\n');
var j,k;
var ret = [];
var pat = /(([^'",]*)|('[^']*')|("[^"]*")),/g;
var patnum = /^\s*(([+-]?[0-9]+(\.[0-9]*)?(e[+-]?[0-9]+)?)|([+-]?[0-9]*(\.[0-9]+)?(e[+-]?[0-9]+)?))\s*$/;
var stripper = function(n) { return n.substr(0,n.length-1); }
var count = 0;
for(k=0;k<foo.length;k++) {
var bar = (foo[k]+",").match(pat),baz;
if(bar.length>0) {
ret[count] = [];
for(j=0;j<bar.length;j++) {
baz = stripper(bar[j]);
if(patnum.test(baz)) { ret[count][j] = parseFloat(baz); }
else ret[count][j] = baz;
}
count++;
}
}
return ret;
}
numeric.toCSV = function toCSV(A) {
var s = numeric.dim(A);
var i,j,m,n,row,ret;
m = s[0];
n = s[1];
ret = [];
for(i=0;i<m;i++) {
row = [];
for(j=0;j<m;j++) { row[j] = A[i][j].toString(); }
ret[i] = row.join(', ');
}
return ret.join('\n')+'\n';
}
numeric.getURL = function getURL(url) {
var client = new XMLHttpRequest();
client.open("GET",url,false);
client.send();
return client;
}
numeric.imageURL = function imageURL(img) {
function base64(A) {
var n = A.length, i,x,y,z,p,q,r,s;
var key = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/=";
var ret = "";
for(i=0;i<n;i+=3) {
x = A[i];
y = A[i+1];
z = A[i+2];
p = x >> 2;
q = ((x & 3) << 4) + (y >> 4);
r = ((y & 15) << 2) + (z >> 6);
s = z & 63;
if(i+1>=n) { r = s = 64; }
else if(i+2>=n) { s = 64; }
ret += key.charAt(p) + key.charAt(q) + key.charAt(r) + key.charAt(s);
}
return ret;
}
function crc32Array (a,from,to) {
if(typeof from === "undefined") { from = 0; }
if(typeof to === "undefined") { to = a.length; }
var table = [0x00000000, 0x77073096, 0xEE0E612C, 0x990951BA, 0x076DC419, 0x706AF48F, 0xE963A535, 0x9E6495A3,
0x0EDB8832, 0x79DCB8A4, 0xE0D5E91E, 0x97D2D988, 0x09B64C2B, 0x7EB17CBD, 0xE7B82D07, 0x90BF1D91,
0x1DB71064, 0x6AB020F2, 0xF3B97148, 0x84BE41DE, 0x1ADAD47D, 0x6DDDE4EB, 0xF4D4B551, 0x83D385C7,
0x136C9856, 0x646BA8C0, 0xFD62F97A, 0x8A65C9EC, 0x14015C4F, 0x63066CD9, 0xFA0F3D63, 0x8D080DF5,
0x3B6E20C8, 0x4C69105E, 0xD56041E4, 0xA2677172, 0x3C03E4D1, 0x4B04D447, 0xD20D85FD, 0xA50AB56B,
0x35B5A8FA, 0x42B2986C, 0xDBBBC9D6, 0xACBCF940, 0x32D86CE3, 0x45DF5C75, 0xDCD60DCF, 0xABD13D59,
0x26D930AC, 0x51DE003A, 0xC8D75180, 0xBFD06116, 0x21B4F4B5, 0x56B3C423, 0xCFBA9599, 0xB8BDA50F,
0x2802B89E, 0x5F058808, 0xC60CD9B2, 0xB10BE924, 0x2F6F7C87, 0x58684C11, 0xC1611DAB, 0xB6662D3D,
0x76DC4190, 0x01DB7106, 0x98D220BC, 0xEFD5102A, 0x71B18589, 0x06B6B51F, 0x9FBFE4A5, 0xE8B8D433,
0x7807C9A2, 0x0F00F934, 0x9609A88E, 0xE10E9818, 0x7F6A0DBB, 0x086D3D2D, 0x91646C97, 0xE6635C01,
0x6B6B51F4, 0x1C6C6162, 0x856530D8, 0xF262004E, 0x6C0695ED, 0x1B01A57B, 0x8208F4C1, 0xF50FC457,
0x65B0D9C6, 0x12B7E950, 0x8BBEB8EA, 0xFCB9887C, 0x62DD1DDF, 0x15DA2D49, 0x8CD37CF3, 0xFBD44C65,
0x4DB26158, 0x3AB551CE, 0xA3BC0074, 0xD4BB30E2, 0x4ADFA541, 0x3DD895D7, 0xA4D1C46D, 0xD3D6F4FB,
0x4369E96A, 0x346ED9FC, 0xAD678846, 0xDA60B8D0, 0x44042D73, 0x33031DE5, 0xAA0A4C5F, 0xDD0D7CC9,
0x5005713C, 0x270241AA, 0xBE0B1010, 0xC90C2086, 0x5768B525, 0x206F85B3, 0xB966D409, 0xCE61E49F,
0x5EDEF90E, 0x29D9C998, 0xB0D09822, 0xC7D7A8B4, 0x59B33D17, 0x2EB40D81, 0xB7BD5C3B, 0xC0BA6CAD,
0xEDB88320, 0x9ABFB3B6, 0x03B6E20C, 0x74B1D29A, 0xEAD54739, 0x9DD277AF, 0x04DB2615, 0x73DC1683,
0xE3630B12, 0x94643B84, 0x0D6D6A3E, 0x7A6A5AA8, 0xE40ECF0B, 0x9309FF9D, 0x0A00AE27, 0x7D079EB1,
0xF00F9344, 0x8708A3D2, 0x1E01F268, 0x6906C2FE, 0xF762575D, 0x806567CB, 0x196C3671, 0x6E6B06E7,
0xFED41B76, 0x89D32BE0, 0x10DA7A5A, 0x67DD4ACC, 0xF9B9DF6F, 0x8EBEEFF9, 0x17B7BE43, 0x60B08ED5,
0xD6D6A3E8, 0xA1D1937E, 0x38D8C2C4, 0x4FDFF252, 0xD1BB67F1, 0xA6BC5767, 0x3FB506DD, 0x48B2364B,
0xD80D2BDA, 0xAF0A1B4C, 0x36034AF6, 0x41047A60, 0xDF60EFC3, 0xA867DF55, 0x316E8EEF, 0x4669BE79,
0xCB61B38C, 0xBC66831A, 0x256FD2A0, 0x5268E236, 0xCC0C7795, 0xBB0B4703, 0x220216B9, 0x5505262F,
0xC5BA3BBE, 0xB2BD0B28, 0x2BB45A92, 0x5CB36A04, 0xC2D7FFA7, 0xB5D0CF31, 0x2CD99E8B, 0x5BDEAE1D,
0x9B64C2B0, 0xEC63F226, 0x756AA39C, 0x026D930A, 0x9C0906A9, 0xEB0E363F, 0x72076785, 0x05005713,
0x95BF4A82, 0xE2B87A14, 0x7BB12BAE, 0x0CB61B38, 0x92D28E9B, 0xE5D5BE0D, 0x7CDCEFB7, 0x0BDBDF21,
0x86D3D2D4, 0xF1D4E242, 0x68DDB3F8, 0x1FDA836E, 0x81BE16CD, 0xF6B9265B, 0x6FB077E1, 0x18B74777,
0x88085AE6, 0xFF0F6A70, 0x66063BCA, 0x11010B5C, 0x8F659EFF, 0xF862AE69, 0x616BFFD3, 0x166CCF45,
0xA00AE278, 0xD70DD2EE, 0x4E048354, 0x3903B3C2, 0xA7672661, 0xD06016F7, 0x4969474D, 0x3E6E77DB,
0xAED16A4A, 0xD9D65ADC, 0x40DF0B66, 0x37D83BF0, 0xA9BCAE53, 0xDEBB9EC5, 0x47B2CF7F, 0x30B5FFE9,
0xBDBDF21C, 0xCABAC28A, 0x53B39330, 0x24B4A3A6, 0xBAD03605, 0xCDD70693, 0x54DE5729, 0x23D967BF,
0xB3667A2E, 0xC4614AB8, 0x5D681B02, 0x2A6F2B94, 0xB40BBE37, 0xC30C8EA1, 0x5A05DF1B, 0x2D02EF8D];
var crc = -1, y = 0, n = a.length,i;
for (i = from; i < to; i++) {
y = (crc ^ a[i]) & 0xFF;
crc = (crc >>> 8) ^ table[y];
}
return crc ^ (-1);
}
var h = img[0].length, w = img[0][0].length, s1, s2, next,k,length,a,b,i,j,adler32,crc32;
var stream = [
137, 80, 78, 71, 13, 10, 26, 10, // 0: PNG signature
0,0,0,13, // 8: IHDR Chunk length
73, 72, 68, 82, // 12: "IHDR"
(w >> 24) & 255, (w >> 16) & 255, (w >> 8) & 255, w&255, // 16: Width
(h >> 24) & 255, (h >> 16) & 255, (h >> 8) & 255, h&255, // 20: Height
8, // 24: bit depth
2, // 25: RGB
0, // 26: deflate
0, // 27: no filter
0, // 28: no interlace
-1,-2,-3,-4, // 29: CRC
-5,-6,-7,-8, // 33: IDAT Chunk length
73, 68, 65, 84, // 37: "IDAT"
// RFC 1950 header starts here
8, // 41: RFC1950 CMF
29 // 42: RFC1950 FLG
];
crc32 = crc32Array(stream,12,29);
stream[29] = (crc32>>24)&255;
stream[30] = (crc32>>16)&255;
stream[31] = (crc32>>8)&255;
stream[32] = (crc32)&255;
s1 = 1;
s2 = 0;
for(i=0;i<h;i++) {
if(i<h-1) { stream.push(0); }
else { stream.push(1); }
a = (3*w+1+(i===0))&255; b = ((3*w+1+(i===0))>>8)&255;
stream.push(a); stream.push(b);
stream.push((~a)&255); stream.push((~b)&255);
if(i===0) stream.push(0);
for(j=0;j<w;j++) {
for(k=0;k<3;k++) {
a = img[k][i][j];
if(a>255) a = 255;
else if(a<0) a=0;
else a = Math.round(a);
s1 = (s1 + a )%65521;
s2 = (s2 + s1)%65521;
stream.push(a);
}
}
stream.push(0);
}
adler32 = (s2<<16)+s1;
stream.push((adler32>>24)&255);
stream.push((adler32>>16)&255);
stream.push((adler32>>8)&255);
stream.push((adler32)&255);
length = stream.length - 41;
stream[33] = (length>>24)&255;
stream[34] = (length>>16)&255;
stream[35] = (length>>8)&255;
stream[36] = (length)&255;
crc32 = crc32Array(stream,37);
stream.push((crc32>>24)&255);
stream.push((crc32>>16)&255);
stream.push((crc32>>8)&255);
stream.push((crc32)&255);
stream.push(0);
stream.push(0);
stream.push(0);
stream.push(0);
// a = stream.length;
stream.push(73); // I
stream.push(69); // E
stream.push(78); // N
stream.push(68); // D
stream.push(174); // CRC1
stream.push(66); // CRC2
stream.push(96); // CRC3
stream.push(130); // CRC4
return 'data:image/png;base64,'+base64(stream);
}
// 2. Linear algebra with Arrays.
numeric._dim = function _dim(x) {
var ret = [];
while(typeof x === "object") { ret.push(x.length); x = x[0]; }
return ret;
}
numeric.dim = function dim(x) {
var y,z;
if(typeof x === "object") {
y = x[0];
if(typeof y === "object") {
z = y[0];
if(typeof z === "object") {
return numeric._dim(x);
}
return [x.length,y.length];
}
return [x.length];
}
return [];
}
numeric.mapreduce = function mapreduce(body,init) {
return Function('x','accum','_s','_k',
'if(typeof accum === "undefined") accum = '+init+';\n'+
'if(typeof x === "number") { var xi = x; '+body+'; return accum; }\n'+
'if(typeof _s === "undefined") _s = numeric.dim(x);\n'+
'if(typeof _k === "undefined") _k = 0;\n'+
'var _n = _s[_k];\n'+
'var i,xi;\n'+
'if(_k < _s.length-1) {\n'+
' for(i=_n-1;i>=0;i--) {\n'+
' accum = arguments.callee(x[i],accum,_s,_k+1);\n'+
' }'+
' return accum;\n'+
'}\n'+
'for(i=_n-1;i>=1;i-=2) { \n'+
' xi = x[i];\n'+
' '+body+';\n'+
' xi = x[i-1];\n'+
' '+body+';\n'+
'}\n'+
'if(i === 0) {\n'+
' xi = x[i];\n'+
' '+body+'\n'+
'}\n'+
'return accum;'
);
}
numeric.mapreduce2 = function mapreduce2(body,setup) {
return Function('x',
'var n = x.length;\n'+
'var i,xi;\n'+setup+';\n'+
'for(i=n-1;i!==-1;--i) { \n'+
' xi = x[i];\n'+
' '+body+';\n'+
'}\n'+
'return accum;'
);
}
numeric.same = function same(x,y) {
var i,n;
if(!(x instanceof Array) || !(y instanceof Array)) { return false; }
n = x.length;
if(n !== y.length) { return false; }
for(i=0;i<n;i++) {
if(x[i] === y[i]) { continue; }
if(typeof x[i] === "object") { if(!same(x[i],y[i])) return false; }
else { return false; }
}
return true;
}
numeric.rep = function rep(s,v,k) {
if(typeof k === "undefined") { k=0; }
var n = s[k], ret = Array(n), i;
if(k === s.length-1) {
for(i=n-2;i>=0;i-=2) { ret[i+1] = v; ret[i] = v; }
if(i===-1) { ret[0] = v; }
return ret;
}
for(i=n-1;i>=0;i--) { ret[i] = numeric.rep(s,v,k+1); }
return ret;
}
numeric.dotMMsmall = function dotMMsmall(x,y) {
var i,j,k,p,q,r,ret,foo,bar,woo,i0,k0,p0,r0;
p = x.length; q = y.length; r = y[0].length;
ret = Array(p);
for(i=p-1;i>=0;i--) {
foo = Array(r);
bar = x[i];
for(k=r-1;k>=0;k--) {
woo = bar[q-1]*y[q-1][k];
for(j=q-2;j>=1;j-=2) {
i0 = j-1;
woo += bar[j]*y[j][k] + bar[i0]*y[i0][k];
}
if(j===0) { woo += bar[0]*y[0][k]; }
foo[k] = woo;
}
ret[i] = foo;
}
return ret;
}
numeric._getCol = function _getCol(A,j,x) {
var n = A.length, i;
for(i=n-1;i>0;--i) {
x[i] = A[i][j];
--i;
x[i] = A[i][j];
}
if(i===0) x[0] = A[0][j];
}
numeric.dotMMbig = function dotMMbig(x,y){
var gc = numeric._getCol, p = y.length, v = Array(p);
var m = x.length, n = y[0].length, A = new Array(m), xj;
var VV = numeric.dotVV;
var i,j,k,z;
--p;
--m;
for(i=m;i!==-1;--i) A[i] = Array(n);
--n;
for(i=n;i!==-1;--i) {
gc(y,i,v);
for(j=m;j!==-1;--j) {
z=0;
xj = x[j];
A[j][i] = VV(xj,v);
}
}
return A;
}
numeric.dotMV = function dotMV(x,y) {
var p = x.length, q = y.length,i;
var ret = Array(p), dotVV = numeric.dotVV;
for(i=p-1;i>=0;i--) { ret[i] = dotVV(x[i],y); }
return ret;
}
numeric.dotVM = function dotVM(x,y) {
var i,j,k,p,q,r,ret,foo,bar,woo,i0,k0,p0,r0,s1,s2,s3,baz,accum;
p = x.length; q = y[0].length;
ret = Array(q);
for(k=q-1;k>=0;k--) {
woo = x[p-1]*y[p-1][k];
for(j=p-2;j>=1;j-=2) {
i0 = j-1;
woo += x[j]*y[j][k] + x[i0]*y[i0][k];
}
if(j===0) { woo += x[0]*y[0][k]; }
ret[k] = woo;
}
return ret;
}
numeric.dotVV = function dotVV(x,y) {
var i,n=x.length,i1,ret = x[n-1]*y[n-1];
for(i=n-2;i>=1;i-=2) {
i1 = i-1;
ret += x[i]*y[i] + x[i1]*y[i1];
}
if(i===0) { ret += x[0]*y[0]; }
return ret;
}
numeric.dot = function dot(x,y) {
var d = numeric.dim;
switch(d(x).length*1000+d(y).length) {
case 2002:
if(y.length < 10) return numeric.dotMMsmall(x,y);
else return numeric.dotMMbig(x,y);
case 2001: return numeric.dotMV(x,y);
case 1002: return numeric.dotVM(x,y);
case 1001: return numeric.dotVV(x,y);
case 1000: return numeric.mulVS(x,y);
case 1: return numeric.mulSV(x,y);
case 0: return x*y;
default: throw new Error('numeric.dot only works on vectors and matrices');
}
}
numeric.diag = function diag(d) {
var i,i1,j,n = d.length, A = Array(n), Ai;
for(i=n-1;i>=0;i--) {
Ai = Array(n);
i1 = i+2;
for(j=n-1;j>=i1;j-=2) {
Ai[j] = 0;
Ai[j-1] = 0;
}
if(j>i) { Ai[j] = 0; }
Ai[i] = d[i];
for(j=i-1;j>=1;j-=2) {
Ai[j] = 0;
Ai[j-1] = 0;
}
if(j===0) { Ai[0] = 0; }
A[i] = Ai;
}
return A;
}
numeric.getDiag = function(A) {
var n = Math.min(A.length,A[0].length),i,ret = Array(n);
for(i=n-1;i>=1;--i) {
ret[i] = A[i][i];
--i;
ret[i] = A[i][i];
}
if(i===0) {
ret[0] = A[0][0];
}
return ret;
}
numeric.identity = function identity(n) { return numeric.diag(numeric.rep([n],1)); }
numeric.pointwise = function pointwise(params,body,setup) {
if(typeof setup === "undefined") { setup = ""; }
var fun = [];
var k;
var avec = /\[i\]$/,p,thevec = '';
var haveret = false;
for(k=0;k<params.length;k++) {
if(avec.test(params[k])) {
p = params[k].substring(0,params[k].length-3);
thevec = p;
} else { p = params[k]; }
if(p==='ret') haveret = true;
fun.push(p);
}
fun[params.length] = '_s';
fun[params.length+1] = '_k';
fun[params.length+2] = (
'if(typeof _s === "undefined") _s = numeric.dim('+thevec+');\n'+
'if(typeof _k === "undefined") _k = 0;\n'+
'var _n = _s[_k];\n'+
'var i'+(haveret?'':', ret = Array(_n)')+';\n'+
'if(_k < _s.length-1) {\n'+
' for(i=_n-1;i>=0;i--) ret[i] = arguments.callee('+params.join(',')+',_s,_k+1);\n'+
' return ret;\n'+
'}\n'+
setup+'\n'+
'for(i=_n-1;i!==-1;--i) {\n'+
' '+body+'\n'+
'}\n'+
'return ret;'
);
return Function.apply(null,fun);
}
numeric.pointwise2 = function pointwise2(params,body,setup) {
if(typeof setup === "undefined") { setup = ""; }
var fun = [];
var k;
var avec = /\[i\]$/,p,thevec = '';
var haveret = false;
for(k=0;k<params.length;k++) {
if(avec.test(params[k])) {
p = params[k].substring(0,params[k].length-3);
thevec = p;
} else { p = params[k]; }
if(p==='ret') haveret = true;
fun.push(p);
}
fun[params.length] = (
'var _n = '+thevec+'.length;\n'+
'var i'+(haveret?'':', ret = Array(_n)')+';\n'+
setup+'\n'+
'for(i=_n-1;i!==-1;--i) {\n'+
body+'\n'+
'}\n'+
'return ret;'
);
return Function.apply(null,fun);
}
numeric._biforeach = (function _biforeach(x,y,s,k,f) {
if(k === s.length-1) { f(x,y); return; }
var i,n=s[k];
for(i=n-1;i>=0;i--) { _biforeach(typeof x==="object"?x[i]:x,typeof y==="object"?y[i]:y,s,k+1,f); }
});
numeric._biforeach2 = (function _biforeach2(x,y,s,k,f) {
if(k === s.length-1) { return f(x,y); }
var i,n=s[k],ret = Array(n);
for(i=n-1;i>=0;--i) { ret[i] = _biforeach2(typeof x==="object"?x[i]:x,typeof y==="object"?y[i]:y,s,k+1,f); }
return ret;
});
numeric._foreach = (function _foreach(x,s,k,f) {
if(k === s.length-1) { f(x); return; }
var i,n=s[k];
for(i=n-1;i>=0;i--) { _foreach(x[i],s,k+1,f); }
});
numeric._foreach2 = (function _foreach2(x,s,k,f) {
if(k === s.length-1) { return f(x); }
var i,n=s[k], ret = Array(n);
for(i=n-1;i>=0;i--) { ret[i] = _foreach2(x[i],s,k+1,f); }
return ret;
});
/*numeric.anyV = numeric.mapreduce('if(xi) return true;','false');
numeric.allV = numeric.mapreduce('if(!xi) return false;','true');
numeric.any = function(x) { if(typeof x.length === "undefined") return x; return numeric.anyV(x); }
numeric.all = function(x) { if(typeof x.length === "undefined") return x; return numeric.allV(x); }*/
numeric.ops2 = {
add: '+',
sub: '-',
mul: '*',
div: '/',
mod: '%',
and: '&&',
or: '||',
eq: '===',
neq: '!==',
lt: '<',
gt: '>',
leq: '<=',
geq: '>=',
band: '&',
bor: '|',
bxor: '^',
lshift: '<<',
rshift: '>>',
rrshift: '>>>'
};
numeric.opseq = {
addeq: '+=',
subeq: '-=',
muleq: '*=',
diveq: '/=',
modeq: '%=',
lshifteq: '<<=',
rshifteq: '>>=',
rrshifteq: '>>>=',
bandeq: '&=',
boreq: '|=',
bxoreq: '^='
};
numeric.mathfuns = ['abs','acos','asin','atan','ceil','cos',
'exp','floor','log','round','sin','sqrt','tan',
'isNaN','isFinite'];
numeric.mathfuns2 = ['atan2','pow','max','min'];
numeric.ops1 = {
neg: '-',
not: '!',
bnot: '~',
clone: ''
};
numeric.mapreducers = {
any: ['if(xi) return true;','var accum = false;'],
all: ['if(!xi) return false;','var accum = true;'],
sum: ['accum += xi;','var accum = 0;'],
prod: ['accum *= xi;','var accum = 1;'],
norm2Squared: ['accum += xi*xi;','var accum = 0;'],
norminf: ['accum = max(accum,abs(xi));','var accum = 0, max = Math.max, abs = Math.abs;'],
norm1: ['accum += abs(xi)','var accum = 0, abs = Math.abs;'],
sup: ['accum = max(accum,xi);','var accum = -Infinity, max = Math.max;'],
inf: ['accum = min(accum,xi);','var accum = Infinity, min = Math.min;']
};
(function () {
var i,o;
for(i=0;i<numeric.mathfuns2.length;++i) {
o = numeric.mathfuns2[i];
numeric.ops2[o] = o;
}
for(i in numeric.ops2) {
if(numeric.ops2.hasOwnProperty(i)) {
o = numeric.ops2[i];
var code, codeeq, setup = '';
if(numeric.myIndexOf.call(numeric.mathfuns2,i)!==-1) {
setup = 'var '+o+' = Math.'+o+';\n';
code = function(r,x,y) { return r+' = '+o+'('+x+','+y+')'; };
codeeq = function(x,y) { return x+' = '+o+'('+x+','+y+')'; };
} else {
code = function(r,x,y) { return r+' = '+x+' '+o+' '+y; };
if(numeric.opseq.hasOwnProperty(i+'eq')) {
codeeq = function(x,y) { return x+' '+o+'= '+y; };
} else {
codeeq = function(x,y) { return x+' = '+x+' '+o+' '+y; };
}
}
numeric[i+'VV'] = numeric.pointwise2(['x[i]','y[i]'],code('ret[i]','x[i]','y[i]'),setup);
numeric[i+'SV'] = numeric.pointwise2(['x','y[i]'],code('ret[i]','x','y[i]'),setup);
numeric[i+'VS'] = numeric.pointwise2(['x[i]','y'],code('ret[i]','x[i]','y'),setup);
numeric[i] = Function(
'var n = arguments.length, i, x = arguments[0], y;\n'+
'var VV = numeric.'+i+'VV, VS = numeric.'+i+'VS, SV = numeric.'+i+'SV;\n'+
'var dim = numeric.dim;\n'+
'for(i=1;i!==n;++i) { \n'+
' y = arguments[i];\n'+
' if(typeof x === "object") {\n'+
' if(typeof y === "object") x = numeric._biforeach2(x,y,dim(x),0,VV);\n'+
' else x = numeric._biforeach2(x,y,dim(x),0,VS);\n'+
' } else if(typeof y === "object") x = numeric._biforeach2(x,y,dim(y),0,SV);\n'+
' else '+codeeq('x','y')+'\n'+
'}\nreturn x;\n');
numeric[o] = numeric[i];
numeric[i+'eqV'] = numeric.pointwise2(['ret[i]','x[i]'], codeeq('ret[i]','x[i]'),setup);
numeric[i+'eqS'] = numeric.pointwise2(['ret[i]','x'], codeeq('ret[i]','x'),setup);
numeric[i+'eq'] = Function(
'var n = arguments.length, i, x = arguments[0], y;\n'+
'var V = numeric.'+i+'eqV, S = numeric.'+i+'eqS\n'+
'var s = numeric.dim(x);\n'+
'for(i=1;i!==n;++i) { \n'+
' y = arguments[i];\n'+
' if(typeof y === "object") numeric._biforeach(x,y,s,0,V);\n'+
' else numeric._biforeach(x,y,s,0,S);\n'+
'}\nreturn x;\n');
}
}
for(i=0;i<numeric.mathfuns2.length;++i) {
o = numeric.mathfuns2[i];
delete numeric.ops2[o];
}
for(i=0;i<numeric.mathfuns.length;++i) {
o = numeric.mathfuns[i];
numeric.ops1[o] = o;
}
for(i in numeric.ops1) {
if(numeric.ops1.hasOwnProperty(i)) {
setup = '';
o = numeric.ops1[i];
if(numeric.myIndexOf.call(numeric.mathfuns,i)!==-1) {
if(Math.hasOwnProperty(o)) setup = 'var '+o+' = Math.'+o+';\n';
}
numeric[i+'eqV'] = numeric.pointwise2(['ret[i]'],'ret[i] = '+o+'(ret[i]);',setup);
numeric[i+'eq'] = Function('x',
'if(typeof x !== "object") return '+o+'x\n'+
'var i;\n'+
'var V = numeric.'+i+'eqV;\n'+
'var s = numeric.dim(x);\n'+
'numeric._foreach(x,s,0,V);\n'+
'return x;\n');
numeric[i+'V'] = numeric.pointwise2(['x[i]'],'ret[i] = '+o+'(x[i]);',setup);
numeric[i] = Function('x',
'if(typeof x !== "object") return '+o+'(x)\n'+
'var i;\n'+
'var V = numeric.'+i+'V;\n'+
'var s = numeric.dim(x);\n'+
'return numeric._foreach2(x,s,0,V);\n');
}
}
for(i=0;i<numeric.mathfuns.length;++i) {
o = numeric.mathfuns[i];
delete numeric.ops1[o];
}
for(i in numeric.mapreducers) {
if(numeric.mapreducers.hasOwnProperty(i)) {
o = numeric.mapreducers[i];
numeric[i+'V'] = numeric.mapreduce2(o[0],o[1]);
numeric[i] = Function('x','s','k',
o[1]+
'if(typeof x !== "object") {'+
' xi = x;\n'+
o[0]+';\n'+
' return accum;\n'+
'}'+
'if(typeof s === "undefined") s = numeric.dim(x);\n'+
'if(typeof k === "undefined") k = 0;\n'+
'if(k === s.length-1) return numeric.'+i+'V(x);\n'+
'var xi;\n'+
'var n = x.length, i;\n'+
'for(i=n-1;i!==-1;--i) {\n'+
' xi = arguments.callee(x[i]);\n'+
o[0]+';\n'+
'}\n'+
'return accum;\n');
}
}
}());
numeric.truncVV = numeric.pointwise(['x[i]','y[i]'],'ret[i] = round(x[i]/y[i])*y[i];','var round = Math.round;');
numeric.truncVS = numeric.pointwise(['x[i]','y'],'ret[i] = round(x[i]/y)*y;','var round = Math.round;');
numeric.truncSV = numeric.pointwise(['x','y[i]'],'ret[i] = round(x/y[i])*y[i];','var round = Math.round;');
numeric.trunc = function trunc(x,y) {
if(typeof x === "object") {
if(typeof y === "object") return numeric.truncVV(x,y);
return numeric.truncVS(x,y);
}
if (typeof y === "object") return numeric.truncSV(x,y);
return Math.round(x/y)*y;
}
numeric.inv = function inv(x) {
var s = numeric.dim(x), abs = Math.abs, m = s[0], n = s[1];
var A = numeric.clone(x), Ai, Aj;
var I = numeric.identity(m), Ii, Ij;
var i,j,k,x;
for(j=0;j<n;++j) {
var i0 = -1;
var v0 = -1;
for(i=j;i!==m;++i) { k = abs(A[i][j]); if(k>v0) { i0 = i; v0 = k; } }
Aj = A[i0]; A[i0] = A[j]; A[j] = Aj;
Ij = I[i0]; I[i0] = I[j]; I[j] = Ij;
x = Aj[j];
for(k=j;k!==n;++k) Aj[k] /= x;
for(k=n-1;k!==-1;--k) Ij[k] /= x;
for(i=m-1;i!==-1;--i) {
if(i!==j) {
Ai = A[i];
Ii = I[i];
x = Ai[j];
for(k=j+1;k!==n;++k) Ai[k] -= Aj[k]*x;
for(k=n-1;k>0;--k) { Ii[k] -= Ij[k]*x; --k; Ii[k] -= Ij[k]*x; }
if(k===0) Ii[0] -= Ij[0]*x;
}
}
}
return I;
}
numeric.det = function det(x) {
var s = numeric.dim(x);
if(s.length !== 2 || s[0] !== s[1]) { throw new Error('numeric: det() only works on square matrices'); }
var n = s[0], ret = 1,i,j,k,A = numeric.clone(x),Aj,Ai,alpha,temp,k1,k2,k3;
for(j=0;j<n-1;j++) {
k=j;
for(i=j+1;i<n;i++) { if(Math.abs(A[i][j]) > Math.abs(A[k][j])) { k = i; } }
if(k !== j) {
temp = A[k]; A[k] = A[j]; A[j] = temp;
ret *= -1;
}
Aj = A[j];
for(i=j+1;i<n;i++) {
Ai = A[i];
alpha = Ai[j]/Aj[j];
for(k=j+1;k<n-1;k+=2) {
k1 = k+1;
Ai[k] -= Aj[k]*alpha;
Ai[k1] -= Aj[k1]*alpha;
}
if(k!==n) { Ai[k] -= Aj[k]*alpha; }
}
if(Aj[j] === 0) { return 0; }
ret *= Aj[j];
}
return ret*A[j][j];
}
numeric.transpose = function transpose(x) {
var i,j,m = x.length,n = x[0].length, ret=Array(n),A0,A1,Bj;
for(j=0;j<n;j++) ret[j] = Array(m);
for(i=m-1;i>=1;i-=2) {
A1 = x[i];
A0 = x[i-1];
for(j=n-1;j>=1;--j) {
Bj = ret[j]; Bj[i] = A1[j]; Bj[i-1] = A0[j];
--j;
Bj = ret[j]; Bj[i] = A1[j]; Bj[i-1] = A0[j];
}
if(j===0) {
Bj = ret[0]; Bj[i] = A1[0]; Bj[i-1] = A0[0];
}
}
if(i===0) {
A0 = x[0];
for(j=n-1;j>=1;--j) {
ret[j][0] = A0[j];
--j;
ret[j][0] = A0[j];
}
if(j===0) { ret[0][0] = A0[0]; }
}
return ret;
}
numeric.negtranspose = function negtranspose(x) {
var i,j,m = x.length,n = x[0].length, ret=Array(n),A0,A1,Bj;
for(j=0;j<n;j++) ret[j] = Array(m);
for(i=m-1;i>=1;i-=2) {
A1 = x[i];
A0 = x[i-1];
for(j=n-1;j>=1;--j) {
Bj = ret[j]; Bj[i] = -A1[j]; Bj[i-1] = -A0[j];
--j;
Bj = ret[j]; Bj[i] = -A1[j]; Bj[i-1] = -A0[j];
}
if(j===0) {
Bj = ret[0]; Bj[i] = -A1[0]; Bj[i-1] = -A0[0];
}
}
if(i===0) {
A0 = x[0];
for(j=n-1;j>=1;--j) {
ret[j][0] = -A0[j];
--j;
ret[j][0] = -A0[j];
}
if(j===0) { ret[0][0] = -A0[0]; }
}
return ret;
}
numeric._random = function _random(s,k) {
var i,n=s[k],ret=Array(n), rnd;
if(k === s.length-1) {
rnd = Math.random;
for(i=n-1;i>=1;i-=2) {
ret[i] = rnd();
ret[i-1] = rnd();
}
if(i===0) { ret[0] = rnd(); }
return ret;
}
for(i=n-1;i>=0;i--) ret[i] = _random(s,k+1);
return ret;
}
numeric.random = function random(s) { return numeric._random(s,0); }
numeric.norm2 = function norm2(x) { return Math.sqrt(numeric.norm2Squared(x)); }
numeric.linspace = function linspace(a,b,n) {
if(typeof n === "undefined") n = Math.max(Math.round(b-a)+1,1);
if(n<2) { return n===1?[a]:[]; }
var i,ret = Array(n);
n--;
for(i=n;i>=0;i--) { ret[i] = (i*b+(n-i)*a)/n; }
return ret;
}
numeric.getBlock = function getBlock(x,from,to) {
var s = numeric.dim(x);
function foo(x,k) {
var i,a = from[k], n = to[k]-a, ret = Array(n);
if(k === s.length-1) {
for(i=n;i>=0;i--) { ret[i] = x[i+a]; }
return ret;
}
for(i=n;i>=0;i--) { ret[i] = foo(x[i+a],k+1); }
return ret;
}
return foo(x,0);
}
numeric.setBlock = function setBlock(x,from,to,B) {
var s = numeric.dim(x);
function foo(x,y,k) {
var i,a = from[k], n = to[k]-a;
if(k === s.length-1) { for(i=n;i>=0;i--) { x[i+a] = y[i]; } }
for(i=n;i>=0;i--) { foo(x[i+a],y[i],k+1); }
}
foo(x,B,0);
return x;
}
numeric.getRange = function getRange(A,I,J) {
var m = I.length, n = J.length;
var i,j;
var B = Array(m), Bi, AI;
for(i=m-1;i!==-1;--i) {
B[i] = Array(n);
Bi = B[i];
AI = A[I[i]];
for(j=n-1;j!==-1;--j) Bi[j] = AI[J[j]];
}
return B;
}
numeric.blockMatrix = function blockMatrix(X) {
var s = numeric.dim(X);
if(s.length<4) return numeric.blockMatrix([X]);
var m=s[0],n=s[1],M,N,i,j,Xij;
M = 0; N = 0;
for(i=0;i<m;++i) M+=X[i][0].length;
for(j=0;j<n;++j) N+=X[0][j][0].length;
var Z = Array(M);
for(i=0;i<M;++i) Z[i] = Array(N);
var I=0,J,ZI,k,l,Xijk;
for(i=0;i<m;++i) {
J=N;
for(j=n-1;j!==-1;--j) {
Xij = X[i][j];
J -= Xij[0].length;
for(k=Xij.length-1;k!==-1;--k) {
Xijk = Xij[k];
ZI = Z[I+k];
for(l = Xijk.length-1;l!==-1;--l) ZI[J+l] = Xijk[l];
}
}
I += X[i][0].length;
}
return Z;
}
numeric.tensor = function tensor(x,y) {
if(typeof x === "number" || typeof y === "number") return numeric.mul(x,y);
var s1 = numeric.dim(x), s2 = numeric.dim(y);
if(s1.length !== 1 || s2.length !== 1) {
throw new Error('numeric: tensor product is only defined for vectors');
}
var m = s1[0], n = s2[0], A = Array(m), Ai, i,j,xi;
for(i=m-1;i>=0;i--) {
Ai = Array(n);
xi = x[i];
for(j=n-1;j>=3;--j) {
Ai[j] = xi * y[j];
--j;
Ai[j] = xi * y[j];
--j;
Ai[j] = xi * y[j];
--j;
Ai[j] = xi * y[j];
}
while(j>=0) { Ai[j] = xi * y[j]; --j; }
A[i] = Ai;
}
return A;
}
// 3. The Tensor type T
numeric.T = function T(x,y) { this.x = x; this.y = y; }
numeric.t = function t(x,y) { return new numeric.T(x,y); }
numeric.Tbinop = function Tbinop(rr,rc,cr,cc,setup) {
var io = numeric.indexOf;
if(typeof setup !== "string") {
var k;
setup = '';
for(k in numeric) {
if(numeric.hasOwnProperty(k) && (rr.indexOf(k)>=0 || rc.indexOf(k)>=0 || cr.indexOf(k)>=0 || cc.indexOf(k)>=0) && k.length>1) {
setup += 'var '+k+' = numeric.'+k+';\n';
}
}
}
return Function(['y'],
'var x = this;\n'+
'if(!(y instanceof numeric.T)) { y = new numeric.T(y); }\n'+
setup+'\n'+
'if(x.y) {'+
' if(y.y) {'+
' return new numeric.T('+cc+');\n'+
' }\n'+
' return new numeric.T('+cr+');\n'+
'}\n'+
'if(y.y) {\n'+
' return new numeric.T('+rc+');\n'+
'}\n'+
'return new numeric.T('+rr+');\n'
);
}
numeric.T.prototype.add = numeric.Tbinop(
'add(x.x,y.x)',
'add(x.x,y.x),y.y',
'add(x.x,y.x),x.y',
'add(x.x,y.x),add(x.y,y.y)');
numeric.T.prototype.sub = numeric.Tbinop(
'sub(x.x,y.x)',
'sub(x.x,y.x),neg(y.y)',
'sub(x.x,y.x),x.y',
'sub(x.x,y.x),sub(x.y,y.y)');
numeric.T.prototype.mul = numeric.Tbinop(
'mul(x.x,y.x)',
'mul(x.x,y.x),mul(x.x,y.y)',
'mul(x.x,y.x),mul(x.y,y.x)',
'sub(mul(x.x,y.x),mul(x.y,y.y)),add(mul(x.x,y.y),mul(x.y,y.x))');
numeric.T.prototype.reciprocal = function reciprocal() {
var mul = numeric.mul, div = numeric.div;
if(this.y) {
var d = numeric.add(mul(this.x,this.x),mul(this.y,this.y));
return new numeric.T(div(this.x,d),div(numeric.neg(this.y),d));
}
return new T(div(1,this.x));
}
numeric.T.prototype.div = function div(y) {
if(!(y instanceof numeric.T)) y = new numeric.T(y);
if(y.y) { return this.mul(y.reciprocal()); }
var div = numeric.div;
if(this.y) { return new numeric.T(div(this.x,y.x),div(this.y,y.x)); }
return new numeric.T(div(this.x,y.x));
}
numeric.T.prototype.dot = numeric.Tbinop(
'dot(x.x,y.x)',
'dot(x.x,y.x),dot(x.x,y.y)',
'dot(x.x,y.x),dot(x.y,y.x)',
'sub(dot(x.x,y.x),dot(x.y,y.y)),add(dot(x.x,y.y),dot(x.y,y.x))'
);
numeric.T.prototype.transpose = function transpose() {
var t = numeric.transpose, x = this.x, y = this.y;
if(y) { return new numeric.T(t(x),t(y)); }
return new numeric.T(t(x));
}
numeric.T.prototype.transjugate = function transjugate() {
var t = numeric.transpose, x = this.x, y = this.y;
if(y) { return new numeric.T(t(x),numeric.negtranspose(y)); }
return new numeric.T(t(x));
}
numeric.Tunop = function Tunop(r,c,s) {
if(typeof s !== "string") { s = ''; }
return Function(
'var x = this;\n'+
s+'\n'+
'if(x.y) {'+
' '+c+';\n'+
'}\n'+
r+';\n'
);
}
numeric.T.prototype.exp = numeric.Tunop(
'return new numeric.T(ex)',
'return new numeric.T(mul(cos(x.y),ex),mul(sin(x.y),ex))',
'var ex = numeric.exp(x.x), cos = numeric.cos, sin = numeric.sin, mul = numeric.mul;');
numeric.T.prototype.conj = numeric.Tunop(
'return new numeric.T(x.x);',
'return new numeric.T(x.x,numeric.neg(x.y));');
numeric.T.prototype.neg = numeric.Tunop(
'return new numeric.T(neg(x.x));',
'return new numeric.T(neg(x.x),neg(x.y));',
'var neg = numeric.neg;');
numeric.T.prototype.sin = numeric.Tunop(
'return new numeric.T(numeric.sin(x.x))',
'return x.exp().sub(x.neg().exp()).div(new numeric.T(0,2));');
numeric.T.prototype.cos = numeric.Tunop(
'return new numeric.T(numeric.cos(x.x))',
'return x.exp().add(x.neg().exp()).div(2);');
numeric.T.prototype.abs = numeric.Tunop(
'return new numeric.T(numeric.abs(x.x));',
'return new numeric.T(numeric.sqrt(numeric.add(mul(x.x,x.x),mul(x.y,x.y))));',
'var mul = numeric.mul;');
numeric.T.prototype.log = numeric.Tunop(
'return new numeric.T(numeric.log(x.x));',
'var theta = new numeric.T(numeric.atan2(x.y,x.x)), r = x.abs();\n'+
'return new numeric.T(numeric.log(r.x),theta.x);');
numeric.T.prototype.norm2 = numeric.Tunop(
'return numeric.norm2(x.x);',
'var f = numeric.norm2Squared;\n'+
'return Math.sqrt(f(x.x)+f(x.y));');
numeric.T.prototype.inv = function inv() {
var A = this;
if(typeof A.y === "undefined") { return new numeric.T(numeric.inv(A.x)); }
var n = A.x.length, i, j, k;
var Rx = numeric.identity(n),Ry = numeric.rep([n,n],0);
var Ax = numeric.clone(A.x), Ay = numeric.clone(A.y);
var Aix, Aiy, Ajx, Ajy, Rix, Riy, Rjx, Rjy;
var i,j,k,d,d1,ax,ay,bx,by,temp;
for(i=0;i<n;i++) {
ax = Ax[i][i]; ay = Ay[i][i];
d = ax*ax+ay*ay;
k = i;
for(j=i+1;j<n;j++) {
ax = Ax[j][i]; ay = Ay[j][i];
d1 = ax*ax+ay*ay;
if(d1 > d) { k=j; d = d1; }
}
if(k!==i) {
temp = Ax[i]; Ax[i] = Ax[k]; Ax[k] = temp;
temp = Ay[i]; Ay[i] = Ay[k]; Ay[k] = temp;
temp = Rx[i]; Rx[i] = Rx[k]; Rx[k] = temp;
temp = Ry[i]; Ry[i] = Ry[k]; Ry[k] = temp;
}
Aix = Ax[i]; Aiy = Ay[i];
Rix = Rx[i]; Riy = Ry[i];
ax = Aix[i]; ay = Aiy[i];
for(j=i+1;j<n;j++) {
bx = Aix[j]; by = Aiy[j];
Aix[j] = (bx*ax+by*ay)/d;
Aiy[j] = (by*ax-bx*ay)/d;
}
for(j=0;j<n;j++) {
bx = Rix[j]; by = Riy[j];
Rix[j] = (bx*ax+by*ay)/d;
Riy[j] = (by*ax-bx*ay)/d;
}
for(j=i+1;j<n;j++) {
Ajx = Ax[j]; Ajy = Ay[j];
Rjx = Rx[j]; Rjy = Ry[j];
ax = Ajx[i]; ay = Ajy[i];
for(k=i+1;k<n;k++) {
bx = Aix[k]; by = Aiy[k];
Ajx[k] -= bx*ax-by*ay;
Ajy[k] -= by*ax+bx*ay;
}
for(k=0;k<n;k++) {
bx = Rix[k]; by = Riy[k];
Rjx[k] -= bx*ax-by*ay;
Rjy[k] -= by*ax+bx*ay;
}
}
}
for(i=n-1;i>0;i--) {
Rix = Rx[i]; Riy = Ry[i];
for(j=i-1;j>=0;j--) {
Rjx = Rx[j]; Rjy = Ry[j];
ax = Ax[j][i]; ay = Ay[j][i];
for(k=n-1;k>=0;k--) {
bx = Rix[k]; by = Riy[k];
Rjx[k] -= ax*bx - ay*by;
Rjy[k] -= ax*by + ay*bx;
}
}
}
return new numeric.T(Rx,Ry);
}
numeric.T.prototype.get = function get(i) {
var x = this.x, y = this.y, k = 0, ik, n = i.length;
if(y) {
while(k<n) {
ik = i[k];
x = x[ik];
y = y[ik];
k++;
}
return new numeric.T(x,y);
}
while(k<n) {
ik = i[k];
x = x[ik];
k++;
}
return new numeric.T(x);
}
numeric.T.prototype.set = function set(i,v) {
var x = this.x, y = this.y, k = 0, ik, n = i.length, vx = v.x, vy = v.y;
if(n===0) {
if(vy) { this.y = vy; }
else if(y) { this.y = undefined; }
this.x = x;
return this;
}
if(vy) {
if(y) { /* ok */ }
else {
y = numeric.rep(numeric.dim(x),0);
this.y = y;
}
while(k<n-1) {
ik = i[k];
x = x[ik];
y = y[ik];
k++;
}
ik = i[k];
x[ik] = vx;
y[ik] = vy;
return this;
}
if(y) {
while(k<n-1) {
ik = i[k];
x = x[ik];
y = y[ik];
k++;
}
ik = i[k];
x[ik] = vx;
if(vx instanceof Array) y[ik] = numeric.rep(numeric.dim(vx),0);
else y[ik] = 0;
return this;
}
while(k<n-1) {
ik = i[k];
x = x[ik];
k++;
}
ik = i[k];
x[ik] = vx;
return this;
}
numeric.T.prototype.getRows = function getRows(i0,i1) {
var n = i1-i0+1, j;
var rx = Array(n), ry, x = this.x, y = this.y;
for(j=i0;j<=i1;j++) { rx[j-i0] = x[j]; }
if(y) {
ry = Array(n);
for(j=i0;j<=i1;j++) { ry[j-i0] = y[j]; }
return new numeric.T(rx,ry);
}
return new numeric.T(rx);
}
numeric.T.prototype.setRows = function setRows(i0,i1,A) {
var j;
var rx = this.x, ry = this.y, x = A.x, y = A.y;
for(j=i0;j<=i1;j++) { rx[j] = x[j-i0]; }
if(y) {
if(!ry) { ry = numeric.rep(numeric.dim(rx),0); this.y = ry; }
for(j=i0;j<=i1;j++) { ry[j] = y[j-i0]; }
} else if(ry) {
for(j=i0;j<=i1;j++) { ry[j] = numeric.rep([x[j-i0].length],0); }
}
return this;
}
numeric.T.prototype.getRow = function getRow(k) {
var x = this.x, y = this.y;
if(y) { return new numeric.T(x[k],y[k]); }
return new numeric.T(x[k]);
}
numeric.T.prototype.setRow = function setRow(i,v) {
var rx = this.x, ry = this.y, x = v.x, y = v.y;
rx[i] = x;
if(y) {
if(!ry) { ry = numeric.rep(numeric.dim(rx),0); this.y = ry; }
ry[i] = y;
} else if(ry) {
ry = numeric.rep([x.length],0);
}
return this;
}
numeric.T.prototype.getBlock = function getBlock(from,to) {
var x = this.x, y = this.y, b = numeric.getBlock;
if(y) { return new numeric.T(b(x,from,to),b(y,from,to)); }
return new numeric.T(b(x,from,to));
}
numeric.T.prototype.setBlock = function setBlock(from,to,A) {
if(!(A instanceof numeric.T)) A = new numeric.T(A);
var x = this.x, y = this.y, b = numeric.setBlock, Ax = A.x, Ay = A.y;
if(Ay) {
if(!y) { this.y = numeric.rep(numeric.dim(this),0); y = this.y; }
b(x,from,to,Ax);
b(y,from,to,Ay);
return this;
}
b(x,from,to,Ax);
if(y) b(y,from,to,numeric.rep(numeric.dim(Ax),0));
}
numeric.T.rep = function rep(s,v) {
var T = numeric.T;
if(!(v instanceof T)) v = new T(v);
var x = v.x, y = v.y, r = numeric.rep;
if(y) return new T(r(s,x),r(s,y));
return new T(r(s,x));
}
numeric.T.diag = function diag(d) {
if(!(d instanceof numeric.T)) d = new numeric.T(d);
var x = d.x, y = d.y, diag = numeric.diag;
if(y) return new numeric.T(diag(x),diag(y));
return new numeric.T(diag(x));
}
numeric.T.eig = function eig() {
if(this.y) { throw new Error('eig: not implemented for complex matrices.'); }
return numeric.eig(this.x);
}
numeric.T.identity = function identity(n) { return new numeric.T(numeric.identity(n)); }
numeric.T.prototype.getDiag = function getDiag() {
var n = numeric;
var x = this.x, y = this.y;
if(y) { return new n.T(n.getDiag(x),n.getDiag(y)); }
return new n.T(n.getDiag(x));
}
// 4. Eigenvalues of real matrices
numeric.house = function house(x) {
var v = numeric.clone(x);
var s = x[0] >= 0 ? 1 : -1;
var alpha = s*numeric.norm2(x);
v[0] += alpha;
var foo = numeric.norm2(v);
if(foo === 0) { /* this should not happen */ throw new Error('eig: internal error'); }
return numeric.div(v,foo);
}
numeric.toUpperHessenberg = function toUpperHessenberg(me) {
var s = numeric.dim(me);
if(s.length !== 2 || s[0] !== s[1]) { throw new Error('numeric: toUpperHessenberg() only works on square matrices'); }
var m = s[0], i,j,k,x,v,A = numeric.clone(me),B,C,Ai,Ci,Q = numeric.identity(m),Qi;
for(j=0;j<m-2;j++) {
x = Array(m-j-1);
for(i=j+1;i<m;i++) { x[i-j-1] = A[i][j]; }
if(numeric.norm2(x)>0) {
v = numeric.house(x);
B = numeric.getBlock(A,[j+1,j],[m-1,m-1]);
C = numeric.tensor(v,numeric.dot(v,B));
for(i=j+1;i<m;i++) { Ai = A[i]; Ci = C[i-j-1]; for(k=j;k<m;k++) Ai[k] -= 2*Ci[k-j]; }
B = numeric.getBlock(A,[0,j+1],[m-1,m-1]);
C = numeric.tensor(numeric.dot(B,v),v);
for(i=0;i<m;i++) { Ai = A[i]; Ci = C[i]; for(k=j+1;k<m;k++) Ai[k] -= 2*Ci[k-j-1]; }
B = Array(m-j-1);
for(i=j+1;i<m;i++) B[i-j-1] = Q[i];
C = numeric.tensor(v,numeric.dot(v,B));
for(i=j+1;i<m;i++) { Qi = Q[i]; Ci = C[i-j-1]; for(k=0;k<m;k++) Qi[k] -= 2*Ci[k]; }
}
}
return {H:A, Q:Q};
}
numeric.epsilon = 2.220446049250313e-16;
numeric.QRFrancis = function(H,maxiter) {
if(typeof maxiter === "undefined") { maxiter = 10000; }
H = numeric.clone(H);
var H0 = numeric.clone(H);
var s = numeric.dim(H),m=s[0],x,v,a,b,c,d,det,tr, Hloc, Q = numeric.identity(m), Qi, Hi, B, C, Ci,i,j,k,iter;
if(m<3) { return {Q:Q, B:[ [0,m-1] ]}; }
var epsilon = numeric.epsilon;
for(iter=0;iter<maxiter;iter++) {
for(j=0;j<m-1;j++) {
if(Math.abs(H[j+1][j]) < epsilon*(Math.abs(H[j][j])+Math.abs(H[j+1][j+1]))) {
var QH1 = numeric.QRFrancis(numeric.getBlock(H,[0,0],[j,j]),maxiter);
var QH2 = numeric.QRFrancis(numeric.getBlock(H,[j+1,j+1],[m-1,m-1]),maxiter);
B = Array(j+1);
for(i=0;i<=j;i++) { B[i] = Q[i]; }
C = numeric.dot(QH1.Q,B);
for(i=0;i<=j;i++) { Q[i] = C[i]; }
B = Array(m-j-1);
for(i=j+1;i<m;i++) { B[i-j-1] = Q[i]; }
C = numeric.dot(QH2.Q,B);
for(i=j+1;i<m;i++) { Q[i] = C[i-j-1]; }
return {Q:Q,B:QH1.B.concat(numeric.add(QH2.B,j+1))};
}
}
a = H[m-2][m-2]; b = H[m-2][m-1];
c = H[m-1][m-2]; d = H[m-1][m-1];
tr = a+d;
det = (a*d-b*c);
Hloc = numeric.getBlock(H, [0,0], [2,2]);
if(tr*tr>=4*det) {
var s1,s2;
s1 = 0.5*(tr+Math.sqrt(tr*tr-4*det));
s2 = 0.5*(tr-Math.sqrt(tr*tr-4*det));
Hloc = numeric.add(numeric.sub(numeric.dot(Hloc,Hloc),
numeric.mul(Hloc,s1+s2)),
numeric.diag(numeric.rep([3],s1*s2)));
} else {
Hloc = numeric.add(numeric.sub(numeric.dot(Hloc,Hloc),
numeric.mul(Hloc,tr)),
numeric.diag(numeric.rep([3],det)));
}
x = [Hloc[0][0],Hloc[1][0],Hloc[2][0]];
v = numeric.house(x);
B = [H[0],H[1],H[2]];
C = numeric.tensor(v,numeric.dot(v,B));
for(i=0;i<3;i++) { Hi = H[i]; Ci = C[i]; for(k=0;k<m;k++) Hi[k] -= 2*Ci[k]; }
B = numeric.getBlock(H, [0,0],[m-1,2]);
C = numeric.tensor(numeric.dot(B,v),v);
for(i=0;i<m;i++) { Hi = H[i]; Ci = C[i]; for(k=0;k<3;k++) Hi[k] -= 2*Ci[k]; }
B = [Q[0],Q[1],Q[2]];
C = numeric.tensor(v,numeric.dot(v,B));
for(i=0;i<3;i++) { Qi = Q[i]; Ci = C[i]; for(k=0;k<m;k++) Qi[k] -= 2*Ci[k]; }
var J;
for(j=0;j<m-2;j++) {
for(k=j;k<=j+1;k++) {
if(Math.abs(H[k+1][k]) < epsilon*(Math.abs(H[k][k])+Math.abs(H[k+1][k+1]))) {
var QH1 = numeric.QRFrancis(numeric.getBlock(H,[0,0],[k,k]),maxiter);
var QH2 = numeric.QRFrancis(numeric.getBlock(H,[k+1,k+1],[m-1,m-1]),maxiter);
B = Array(k+1);
for(i=0;i<=k;i++) { B[i] = Q[i]; }
C = numeric.dot(QH1.Q,B);
for(i=0;i<=k;i++) { Q[i] = C[i]; }
B = Array(m-k-1);
for(i=k+1;i<m;i++) { B[i-k-1] = Q[i]; }
C = numeric.dot(QH2.Q,B);
for(i=k+1;i<m;i++) { Q[i] = C[i-k-1]; }
return {Q:Q,B:QH1.B.concat(numeric.add(QH2.B,k+1))};
}
}
J = Math.min(m-1,j+3);
x = Array(J-j);
for(i=j+1;i<=J;i++) { x[i-j-1] = H[i][j]; }
v = numeric.house(x);
B = numeric.getBlock(H, [j+1,j],[J,m-1]);
C = numeric.tensor(v,numeric.dot(v,B));
for(i=j+1;i<=J;i++) { Hi = H[i]; Ci = C[i-j-1]; for(k=j;k<m;k++) Hi[k] -= 2*Ci[k-j]; }
B = numeric.getBlock(H, [0,j+1],[m-1,J]);
C = numeric.tensor(numeric.dot(B,v),v);
for(i=0;i<m;i++) { Hi = H[i]; Ci = C[i]; for(k=j+1;k<=J;k++) Hi[k] -= 2*Ci[k-j-1]; }
B = Array(J-j);
for(i=j+1;i<=J;i++) B[i-j-1] = Q[i];
C = numeric.tensor(v,numeric.dot(v,B));
for(i=j+1;i<=J;i++) { Qi = Q[i]; Ci = C[i-j-1]; for(k=0;k<m;k++) Qi[k] -= 2*Ci[k]; }
}
}
throw new Error('numeric: eigenvalue iteration does not converge -- increase maxiter?');
}
numeric.eig = function eig(A,maxiter) {
var QH = numeric.toUpperHessenberg(A);
var QB = numeric.QRFrancis(QH.H,maxiter);
var T = numeric.T;
var n = A.length,i,k,flag = false,B = QB.B,H = numeric.dot(QB.Q,numeric.dot(QH.H,numeric.transpose(QB.Q)));
var Q = new T(numeric.dot(QB.Q,QH.Q)),Q0;
var m = B.length,j;
var a,b,c,d,p1,p2,disc,x,y,p,q,n1,n2;
var sqrt = Math.sqrt;
for(k=0;k<m;k++) {
i = B[k][0];
if(i === B[k][1]) {
// nothing
} else {
j = i+1;
a = H[i][i];
b = H[i][j];
c = H[j][i];
d = H[j][j];
if(b === 0 && c === 0) continue;
p1 = -a-d;
p2 = a*d-b*c;
disc = p1*p1-4*p2;
if(disc>=0) {
if(p1<0) x = -0.5*(p1-sqrt(disc));
else x = -0.5*(p1+sqrt(disc));
n1 = (a-x)*(a-x)+b*b;
n2 = c*c+(d-x)*(d-x);
if(n1>n2) {
n1 = sqrt(n1);
p = (a-x)/n1;
q = b/n1;
} else {
n2 = sqrt(n2);
p = c/n2;
q = (d-x)/n2;
}
Q0 = new T([[q,-p],[p,q]]);
Q.setRows(i,j,Q0.dot(Q.getRows(i,j)));
} else {
x = -0.5*p1;
y = 0.5*sqrt(-disc);
n1 = (a-x)*(a-x)+b*b;
n2 = c*c+(d-x)*(d-x);
if(n1>n2) {
n1 = sqrt(n1+y*y);
p = (a-x)/n1;
q = b/n1;
x = 0;
y /= n1;
} else {
n2 = sqrt(n2+y*y);
p = c/n2;
q = (d-x)/n2;
x = y/n2;
y = 0;
}
Q0 = new T([[q,-p],[p,q]],[[x,y],[y,-x]]);
Q.setRows(i,j,Q0.dot(Q.getRows(i,j)));
}
}
}
var R = Q.dot(A).dot(Q.transjugate()), n = A.length, E = numeric.T.identity(n);
for(j=0;j<n;j++) {
if(j>0) {
for(k=j-1;k>=0;k--) {
var Rk = R.get([k,k]), Rj = R.get([j,j]);
if(numeric.neq(Rk.x,Rj.x) || numeric.neq(Rk.y,Rj.y)) {
x = R.getRow(k).getBlock([k],[j-1]);
y = E.getRow(j).getBlock([k],[j-1]);
E.set([j,k],(R.get([k,j]).neg().sub(x.dot(y))).div(Rk.sub(Rj)));
} else {
E.setRow(j,E.getRow(k));
continue;
}
}
}
}
for(j=0;j<n;j++) {
x = E.getRow(j);
E.setRow(j,x.div(x.norm2()));
}
E = E.transpose();
E = Q.transjugate().dot(E);
return { lambda:R.getDiag(), E:E };
};
// 5. Compressed Column Storage matrices
numeric.ccsSparse = function ccsSparse(A) {
var m = A.length,n,foo, i,j, counts = [];
for(i=m-1;i!==-1;--i) {
foo = A[i];
for(j in foo) {
j = parseInt(j);
while(j>=counts.length) counts[counts.length] = 0;
if(foo[j]!==0) counts[j]++;
}
}
var n = counts.length;
var Ai = Array(n+1);
Ai[0] = 0;
for(i=0;i<n;++i) Ai[i+1] = Ai[i] + counts[i];
var Aj = Array(Ai[n]), Av = Array(Ai[n]);
for(i=m-1;i!==-1;--i) {
foo = A[i];
for(j in foo) {
if(foo[j]!==0) {
counts[j]--;
Aj[Ai[j]+counts[j]] = i;
Av[Ai[j]+counts[j]] = foo[j];
}
}
}
return [Ai,Aj,Av];
}
numeric.ccsFull = function ccsFull(A) {
var Ai = A[0], Aj = A[1], Av = A[2], s = numeric.ccsDim(A), m = s[0], n = s[1], i,j,j0,j1,k;
var B = numeric.rep([m,n],0);
for(i=0;i<n;i++) {
j0 = Ai[i];
j1 = Ai[i+1];
for(j=j0;j<j1;++j) { B[Aj[j]][i] = Av[j]; }
}
return B;
}
numeric.ccsTSolve = function ccsTSolve(A,b,x,bj,xj) {
var Ai = A[0], Aj = A[1], Av = A[2],m = Ai.length-1, max = Math.max,n=0;
if(typeof bj === "undefined") x = numeric.rep([m],0);
if(typeof bj === "undefined") bj = numeric.linspace(0,x.length-1);
if(typeof xj === "undefined") xj = [];
function dfs(j) {
var k;
if(x[j] !== 0) return;
x[j] = 1;
for(k=Ai[j];k<Ai[j+1];++k) dfs(Aj[k]);
xj[n] = j;
++n;
}
var i,j,j0,j1,k,l,l0,l1,a;
for(i=bj.length-1;i!==-1;--i) { dfs(bj[i]); }
xj.length = n;
for(i=xj.length-1;i!==-1;--i) { x[xj[i]] = 0; }
for(i=bj.length-1;i!==-1;--i) { j = bj[i]; x[j] = b[j]; }
for(i=xj.length-1;i!==-1;--i) {
j = xj[i];
j0 = Ai[j];
j1 = max(Ai[j+1],j0);
for(k=j0;k!==j1;++k) { if(Aj[k] === j) { x[j] /= Av[k]; break; } }
a = x[j];
for(k=j0;k!==j1;++k) {
l = Aj[k];
if(l !== j) x[l] -= a*Av[k];
}
}
return x;
}
numeric.ccsDFS = function ccsDFS(n) {
this.k = Array(n);
this.k1 = Array(n);
this.j = Array(n);
}
numeric.ccsDFS.prototype.dfs = function dfs(J,Ai,Aj,x,xj,Pinv) {
var m = 0,foo,n=xj.length;
var k = this.k, k1 = this.k1, j = this.j,km,k11;
if(x[J]!==0) return;
x[J] = 1;
j[0] = J;
k[0] = km = Ai[J];
k1[0] = k11 = Ai[J+1];
while(1) {
if(km >= k11) {
xj[n] = j[m];
if(m===0) return;
++n;
--m;
km = k[m];
k11 = k1[m];
} else {
foo = Pinv[Aj[km]];
if(x[foo] === 0) {
x[foo] = 1;
k[m] = km;
++m;
j[m] = foo;
km = Ai[foo];
k1[m] = k11 = Ai[foo+1];
} else ++km;
}
}
}
numeric.ccsLPSolve = function ccsLPSolve(A,B,x,xj,I,Pinv,dfs) {
var Ai = A[0], Aj = A[1], Av = A[2],m = Ai.length-1, n=0;
var Bi = B[0], Bj = B[1], Bv = B[2];
var i,i0,i1,j,J,j0,j1,k,l,l0,l1,a;
i0 = Bi[I];
i1 = Bi[I+1];
xj.length = 0;
for(i=i0;i<i1;++i) { dfs.dfs(Pinv[Bj[i]],Ai,Aj,x,xj,Pinv); }
for(i=xj.length-1;i!==-1;--i) { x[xj[i]] = 0; }
for(i=i0;i!==i1;++i) { j = Pinv[Bj[i]]; x[j] = Bv[i]; }
for(i=xj.length-1;i!==-1;--i) {
j = xj[i];
j0 = Ai[j];
j1 = Ai[j+1];
for(k=j0;k<j1;++k) { if(Pinv[Aj[k]] === j) { x[j] /= Av[k]; break; } }
a = x[j];
for(k=j0;k<j1;++k) {
l = Pinv[Aj[k]];
if(l !== j) x[l] -= a*Av[k];
}
}
return x;
}
numeric.ccsLUP1 = function ccsLUP1(A,threshold) {
var m = A[0].length-1;
var L = [numeric.rep([m+1],0),[],[]], U = [numeric.rep([m+1], 0),[],[]];
var Li = L[0], Lj = L[1], Lv = L[2], Ui = U[0], Uj = U[1], Uv = U[2];
var x = numeric.rep([m],0), xj = numeric.rep([m],0);
var i,j,k,j0,j1,a,e,c,d,K;
var sol = numeric.ccsLPSolve, max = Math.max, abs = Math.abs;
var P = numeric.linspace(0,m-1),Pinv = numeric.linspace(0,m-1);
var dfs = new numeric.ccsDFS(m);
if(typeof threshold === "undefined") { threshold = 1; }
for(i=0;i<m;++i) {
sol(L,A,x,xj,i,Pinv,dfs);
a = -1;
e = -1;
for(j=xj.length-1;j!==-1;--j) {
k = xj[j];
if(k <= i) continue;
c = abs(x[k]);
if(c > a) { e = k; a = c; }
}
if(abs(x[i])<threshold*a) {
j = P[i];
a = P[e];
P[i] = a; Pinv[a] = i;
P[e] = j; Pinv[j] = e;
a = x[i]; x[i] = x[e]; x[e] = a;
}
a = Li[i];
e = Ui[i];
d = x[i];
Lj[a] = P[i];
Lv[a] = 1;
++a;
for(j=xj.length-1;j!==-1;--j) {
k = xj[j];
c = x[k];
xj[j] = 0;
x[k] = 0;
if(k<=i) { Uj[e] = k; Uv[e] = c; ++e; }
else { Lj[a] = P[k]; Lv[a] = c/d; ++a; }
}
Li[i+1] = a;
Ui[i+1] = e;
}
for(j=Lj.length-1;j!==-1;--j) { Lj[j] = Pinv[Lj[j]]; }
return {L:L, U:U, P:P, Pinv:Pinv};
}
numeric.ccsDFS0 = function ccsDFS0(n) {
this.k = Array(n);
this.k1 = Array(n);
this.j = Array(n);
}
numeric.ccsDFS0.prototype.dfs = function dfs(J,Ai,Aj,x,xj,Pinv,P) {
var m = 0,foo,n=xj.length;
var k = this.k, k1 = this.k1, j = this.j,km,k11;
if(x[J]!==0) return;
x[J] = 1;
j[0] = J;
k[0] = km = Ai[Pinv[J]];
k1[0] = k11 = Ai[Pinv[J]+1];
while(1) {
if(isNaN(km)) throw new Error("Ow!");
if(km >= k11) {
xj[n] = Pinv[j[m]];
if(m===0) return;
++n;
--m;
km = k[m];
k11 = k1[m];
} else {
foo = Aj[km];
if(x[foo] === 0) {
x[foo] = 1;
k[m] = km;
++m;
j[m] = foo;
foo = Pinv[foo];
km = Ai[foo];
k1[m] = k11 = Ai[foo+1];
} else ++km;
}
}
}
numeric.ccsLPSolve0 = function ccsLPSolve0(A,B,y,xj,I,Pinv,P,dfs) {
var Ai = A[0], Aj = A[1], Av = A[2],m = Ai.length-1, n=0;
var Bi = B[0], Bj = B[1], Bv = B[2];
var i,i0,i1,j,J,j0,j1,k,l,l0,l1,a;
i0 = Bi[I];
i1 = Bi[I+1];
xj.length = 0;
for(i=i0;i<i1;++i) { dfs.dfs(Bj[i],Ai,Aj,y,xj,Pinv,P); }
for(i=xj.length-1;i!==-1;--i) { j = xj[i]; y[P[j]] = 0; }
for(i=i0;i!==i1;++i) { j = Bj[i]; y[j] = Bv[i]; }
for(i=xj.length-1;i!==-1;--i) {
j = xj[i];
l = P[j];
j0 = Ai[j];
j1 = Ai[j+1];
for(k=j0;k<j1;++k) { if(Aj[k] === l) { y[l] /= Av[k]; break; } }
a = y[l];
for(k=j0;k<j1;++k) y[Aj[k]] -= a*Av[k];
y[l] = a;
}
}
numeric.ccsLUP0 = function ccsLUP0(A,threshold) {
var m = A[0].length-1;
var L = [numeric.rep([m+1],0),[],[]], U = [numeric.rep([m+1], 0),[],[]];
var Li = L[0], Lj = L[1], Lv = L[2], Ui = U[0], Uj = U[1], Uv = U[2];
var y = numeric.rep([m],0), xj = numeric.rep([m],0);
var i,j,k,j0,j1,a,e,c,d,K;
var sol = numeric.ccsLPSolve0, max = Math.max, abs = Math.abs;
var P = numeric.linspace(0,m-1),Pinv = numeric.linspace(0,m-1);
var dfs = new numeric.ccsDFS0(m);
if(typeof threshold === "undefined") { threshold = 1; }
for(i=0;i<m;++i) {
sol(L,A,y,xj,i,Pinv,P,dfs);
a = -1;
e = -1;
for(j=xj.length-1;j!==-1;--j) {
k = xj[j];
if(k <= i) continue;
c = abs(y[P[k]]);
if(c > a) { e = k; a = c; }
}
if(abs(y[P[i]])<threshold*a) {
j = P[i];
a = P[e];
P[i] = a; Pinv[a] = i;
P[e] = j; Pinv[j] = e;
}
a = Li[i];
e = Ui[i];
d = y[P[i]];
Lj[a] = P[i];
Lv[a] = 1;
++a;
for(j=xj.length-1;j!==-1;--j) {
k = xj[j];
c = y[P[k]];
xj[j] = 0;
y[P[k]] = 0;
if(k<=i) { Uj[e] = k; Uv[e] = c; ++e; }
else { Lj[a] = P[k]; Lv[a] = c/d; ++a; }
}
Li[i+1] = a;
Ui[i+1] = e;
}
for(j=Lj.length-1;j!==-1;--j) { Lj[j] = Pinv[Lj[j]]; }
return {L:L, U:U, P:P, Pinv:Pinv};
}
numeric.ccsLUP = numeric.ccsLUP0;
numeric.ccsDim = function ccsDim(A) { return [numeric.sup(A[1])+1,A[0].length-1]; }
numeric.ccsGetBlock = function ccsGetBlock(A,i,j) {
var s = numeric.ccsDim(A),m=s[0],n=s[1];
if(typeof i === "undefined") { i = numeric.linspace(0,m-1); }
else if(typeof i === "number") { i = [i]; }
if(typeof j === "undefined") { j = numeric.linspace(0,n-1); }
else if(typeof j === "number") { j = [j]; }
var p,p0,p1,P = i.length,q,Q = j.length,r,jq,ip;
var Bi = numeric.rep([n],0), Bj=[], Bv=[], B = [Bi,Bj,Bv];
var Ai = A[0], Aj = A[1], Av = A[2];
var x = numeric.rep([m],0),count=0,flags = numeric.rep([m],0);
for(q=0;q<Q;++q) {
jq = j[q];
var q0 = Ai[jq];
var q1 = Ai[jq+1];
for(p=q0;p<q1;++p) {
r = Aj[p];
flags[r] = 1;
x[r] = Av[p];
}
for(p=0;p<P;++p) {
ip = i[p];
if(flags[ip]) {
Bj[count] = p;
Bv[count] = x[i[p]];
++count;
}
}
for(p=q0;p<q1;++p) {
r = Aj[p];
flags[r] = 0;
}
Bi[q+1] = count;
}
return B;
}
numeric.ccsDot = function ccsDot(A,B) {
var Ai = A[0], Aj = A[1], Av = A[2];
var Bi = B[0], Bj = B[1], Bv = B[2];
var sA = numeric.ccsDim(A), sB = numeric.ccsDim(B);
var m = sA[0], n = sA[1], o = sB[1];
var x = numeric.rep([m],0), flags = numeric.rep([m],0), xj = Array(m);
var Ci = numeric.rep([o],0), Cj = [], Cv = [], C = [Ci,Cj,Cv];
var i,j,k,j0,j1,i0,i1,l,p,a,b;
for(k=0;k!==o;++k) {
j0 = Bi[k];
j1 = Bi[k+1];
p = 0;
for(j=j0;j<j1;++j) {
a = Bj[j];
b = Bv[j];
i0 = Ai[a];
i1 = Ai[a+1];
for(i=i0;i<i1;++i) {
l = Aj[i];
if(flags[l]===0) {
xj[p] = l;
flags[l] = 1;
p = p+1;
}
x[l] = x[l] + Av[i]*b;
}
}
j0 = Ci[k];
j1 = j0+p;
Ci[k+1] = j1;
for(j=p-1;j!==-1;--j) {
b = j0+j;
i = xj[j];
Cj[b] = i;
Cv[b] = x[i];
flags[i] = 0;
x[i] = 0;
}
Ci[k+1] = Ci[k]+p;
}
return C;
}
numeric.ccsLUPSolve = function ccsLUPSolve(LUP,B) {
var L = LUP.L, U = LUP.U, P = LUP.P;
var Bi = B[0];
var flag = false;
if(typeof Bi !== "object") { B = [[0,B.length],numeric.linspace(0,B.length-1),B]; Bi = B[0]; flag = true; }
var Bj = B[1], Bv = B[2];
var n = L[0].length-1, m = Bi.length-1;
var x = numeric.rep([n],0), xj = Array(n);
var b = numeric.rep([n],0), bj = Array(n);
var Xi = numeric.rep([m+1],0), Xj = [], Xv = [];
var sol = numeric.ccsTSolve;
var i,j,j0,j1,k,J,N=0;
for(i=0;i<m;++i) {
k = 0;
j0 = Bi[i];
j1 = Bi[i+1];
for(j=j0;j<j1;++j) {
J = LUP.Pinv[Bj[j]];
bj[k] = J;
b[J] = Bv[j];
++k;
}
bj.length = k;
sol(L,b,x,bj,xj);
for(j=bj.length-1;j!==-1;--j) b[bj[j]] = 0;
sol(U,x,b,xj,bj);
if(flag) return b;
for(j=xj.length-1;j!==-1;--j) x[xj[j]] = 0;
for(j=bj.length-1;j!==-1;--j) {
J = bj[j];
Xj[N] = J;
Xv[N] = b[J];
b[J] = 0;
++N;
}
Xi[i+1] = N;
}
return [Xi,Xj,Xv];
}
numeric.ccsbinop = function ccsbinop(body,setup) {
if(typeof setup === "undefined") setup='';
return Function('X','Y',
'var Xi = X[0], Xj = X[1], Xv = X[2];\n'+
'var Yi = Y[0], Yj = Y[1], Yv = Y[2];\n'+
'var n = Xi.length-1,m = Math.max(numeric.sup(Xj),numeric.sup(Yj))+1;\n'+
'var Zi = numeric.rep([n+1],0), Zj = [], Zv = [];\n'+
'var x = numeric.rep([m],0),y = numeric.rep([m],0);\n'+
'var xk,yk,zk;\n'+
'var i,j,j0,j1,k,p=0;\n'+
setup+
'for(i=0;i<n;++i) {\n'+
' j0 = Xi[i]; j1 = Xi[i+1];\n'+
' for(j=j0;j!==j1;++j) {\n'+
' k = Xj[j];\n'+
' x[k] = 1;\n'+
' Zj[p] = k;\n'+
' ++p;\n'+
' }\n'+
' j0 = Yi[i]; j1 = Yi[i+1];\n'+
' for(j=j0;j!==j1;++j) {\n'+
' k = Yj[j];\n'+
' y[k] = Yv[j];\n'+
' if(x[k] === 0) {\n'+
' Zj[p] = k;\n'+
' ++p;\n'+
' }\n'+
' }\n'+
' Zi[i+1] = p;\n'+
' j0 = Xi[i]; j1 = Xi[i+1];\n'+
' for(j=j0;j!==j1;++j) x[Xj[j]] = Xv[j];\n'+
' j0 = Zi[i]; j1 = Zi[i+1];\n'+
' for(j=j0;j!==j1;++j) {\n'+
' k = Zj[j];\n'+
' xk = x[k];\n'+
' yk = y[k];\n'+
body+'\n'+
' Zv[j] = zk;\n'+
' }\n'+
' j0 = Xi[i]; j1 = Xi[i+1];\n'+
' for(j=j0;j!==j1;++j) x[Xj[j]] = 0;\n'+
' j0 = Yi[i]; j1 = Yi[i+1];\n'+
' for(j=j0;j!==j1;++j) y[Yj[j]] = 0;\n'+
'}\n'+
'return [Zi,Zj,Zv];'
);
};
(function() {
var k,A,B,C;
for(k in numeric.ops2) {
if(isFinite(eval('1'+numeric.ops2[k]+'0'))) A = '[Y[0],Y[1],numeric.'+k+'(X,Y[2])]';
else A = 'NaN';
if(isFinite(eval('0'+numeric.ops2[k]+'1'))) B = '[X[0],X[1],numeric.'+k+'(X[2],Y)]';
else B = 'NaN';
if(isFinite(eval('1'+numeric.ops2[k]+'0')) && isFinite(eval('0'+numeric.ops2[k]+'1'))) C = 'numeric.ccs'+k+'MM(X,Y)';
else C = 'NaN';
numeric['ccs'+k+'MM'] = numeric.ccsbinop('zk = xk '+numeric.ops2[k]+'yk;');
numeric['ccs'+k] = Function('X','Y',
'if(typeof X === "number") return '+A+';\n'+
'if(typeof Y === "number") return '+B+';\n'+
'return '+C+';\n'
);
}
}());
numeric.ccsScatter = function ccsScatter(A) {
var Ai = A[0], Aj = A[1], Av = A[2];
var n = numeric.sup(Aj)+1,m=Ai.length;
var Ri = numeric.rep([n],0),Rj=Array(m), Rv = Array(m);
var counts = numeric.rep([n],0),i;
for(i=0;i<m;++i) counts[Aj[i]]++;
for(i=0;i<n;++i) Ri[i+1] = Ri[i] + counts[i];
var ptr = Ri.slice(0),k,Aii;
for(i=0;i<m;++i) {
Aii = Aj[i];
k = ptr[Aii];
Rj[k] = Ai[i];
Rv[k] = Av[i];
ptr[Aii]=ptr[Aii]+1;
}
return [Ri,Rj,Rv];
}
numeric.ccsGather = function ccsGather(A) {
var Ai = A[0], Aj = A[1], Av = A[2];
var n = Ai.length-1,m = Aj.length;
var Ri = Array(m), Rj = Array(m), Rv = Array(m);
var i,j,j0,j1,p;
p=0;
for(i=0;i<n;++i) {
j0 = Ai[i];
j1 = Ai[i+1];
for(j=j0;j!==j1;++j) {
Rj[p] = i;
Ri[p] = Aj[j];
Rv[p] = Av[j];
++p;
}
}
return [Ri,Rj,Rv];
}
// The following sparse linear algebra routines are deprecated.
numeric.sdim = function dim(A,ret,k) {
if(typeof ret === "undefined") { ret = []; }
if(typeof A !== "object") return ret;
if(typeof k === "undefined") { k=0; }
if(!(k in ret)) { ret[k] = 0; }
if(A.length > ret[k]) ret[k] = A.length;
var i;
for(i in A) {
if(A.hasOwnProperty(i)) dim(A[i],ret,k+1);
}
return ret;
};
numeric.sclone = function clone(A,k,n) {
if(typeof k === "undefined") { k=0; }
if(typeof n === "undefined") { n = numeric.sdim(A).length; }
var i,ret = Array(A.length);
if(k === n-1) {
for(i in A) { if(A.hasOwnProperty(i)) ret[i] = A[i]; }
return ret;
}
for(i in A) {
if(A.hasOwnProperty(i)) ret[i] = clone(A[i],k+1,n);
}
return ret;
}
numeric.sdiag = function diag(d) {
var n = d.length,i,ret = Array(n),i1,i2,i3;
for(i=n-1;i>=1;i-=2) {
i1 = i-1;
ret[i] = []; ret[i][i] = d[i];
ret[i1] = []; ret[i1][i1] = d[i1];
}
if(i===0) { ret[0] = []; ret[0][0] = d[i]; }
return ret;
}
numeric.sidentity = function identity(n) { return numeric.sdiag(numeric.rep([n],1)); }
numeric.stranspose = function transpose(A) {
var ret = [], n = A.length, i,j,Ai;
for(i in A) {
if(!(A.hasOwnProperty(i))) continue;
Ai = A[i];
for(j in Ai) {
if(!(Ai.hasOwnProperty(j))) continue;
if(typeof ret[j] !== "object") { ret[j] = []; }
ret[j][i] = Ai[j];
}
}
return ret;
}
numeric.sLUP = function LUP(A,tol) {
throw new Error("The function numeric.sLUP had a bug in it and has been removed. Please use the new numeric.ccsLUP function instead.");
};
numeric.sdotMM = function dotMM(A,B) {
var p = A.length, q = B.length, BT = numeric.stranspose(B), r = BT.length, Ai, BTk;
var i,j,k,accum;
var ret = Array(p),reti;
for(i=p-1;i>=0;i--) {
reti = [];
Ai = A[i];
for(k=r-1;k>=0;k--) {
accum = 0;
BTk = BT[k];
for(j in Ai) {
if(!(Ai.hasOwnProperty(j))) continue;
if(j in BTk) { accum += Ai[j]*BTk[j]; }
}
if(accum) reti[k] = accum;
}
ret[i] = reti;
}
return ret;
}
numeric.sdotMV = function dotMV(A,x) {
var p = A.length, Ai, i,j;
var ret = Array(p), accum;
for(i=p-1;i>=0;i--) {
Ai = A[i];
accum = 0;
for(j in Ai) {
if(!(Ai.hasOwnProperty(j))) continue;
if(x[j]) accum += Ai[j]*x[j];
}
if(accum) ret[i] = accum;
}
return ret;
}
numeric.sdotVM = function dotMV(x,A) {
var i,j,Ai,alpha;
var ret = [], accum;
for(i in x) {
if(!x.hasOwnProperty(i)) continue;
Ai = A[i];
alpha = x[i];
for(j in Ai) {
if(!Ai.hasOwnProperty(j)) continue;
if(!ret[j]) { ret[j] = 0; }
ret[j] += alpha*Ai[j];
}
}
return ret;
}
numeric.sdotVV = function dotVV(x,y) {
var i,ret=0;
for(i in x) { if(x[i] && y[i]) ret+= x[i]*y[i]; }
return ret;
}
numeric.sdot = function dot(A,B) {
var m = numeric.sdim(A).length, n = numeric.sdim(B).length;
var k = m*1000+n;
switch(k) {
case 0: return A*B;
case 1001: return numeric.sdotVV(A,B);
case 2001: return numeric.sdotMV(A,B);
case 1002: return numeric.sdotVM(A,B);
case 2002: return numeric.sdotMM(A,B);
default: throw new Error('numeric.sdot not implemented for tensors of order '+m+' and '+n);
}
}
numeric.sscatter = function scatter(V) {
var n = V[0].length, Vij, i, j, m = V.length, A = [], Aj;
for(i=n-1;i>=0;--i) {
if(!V[m-1][i]) continue;
Aj = A;
for(j=0;j<m-2;j++) {
Vij = V[j][i];
if(!Aj[Vij]) Aj[Vij] = [];
Aj = Aj[Vij];
}
Aj[V[j][i]] = V[j+1][i];
}
return A;
}
numeric.sgather = function gather(A,ret,k) {
if(typeof ret === "undefined") ret = [];
if(typeof k === "undefined") k = [];
var n,i,Ai;
n = k.length;
for(i in A) {
if(A.hasOwnProperty(i)) {
k[n] = parseInt(i);
Ai = A[i];
if(typeof Ai === "number") {
if(Ai) {
if(ret.length === 0) {
for(i=n+1;i>=0;--i) ret[i] = [];
}
for(i=n;i>=0;--i) ret[i].push(k[i]);
ret[n+1].push(Ai);
}
} else gather(Ai,ret,k);
}
}
if(k.length>n) k.pop();
return ret;
}
// 6. Coordinate matrices
numeric.cLU = function LU(A) {
var I = A[0], J = A[1], V = A[2];
var p = I.length, m=0, i,j,k,a,b,c;
for(i=0;i<p;i++) if(I[i]>m) m=I[i];
m++;
var L = Array(m), U = Array(m), left = numeric.rep([m],Infinity), right = numeric.rep([m],-Infinity);
var Ui, Uj,alpha;
for(k=0;k<p;k++) {
i = I[k];
j = J[k];
if(j<left[i]) left[i] = j;
if(j>right[i]) right[i] = j;
}
for(i=0;i<m-1;i++) { if(right[i] > right[i+1]) right[i+1] = right[i]; }
for(i=m-1;i>=1;i--) { if(left[i]<left[i-1]) left[i-1] = left[i]; }
var countL = 0, countU = 0;
for(i=0;i<m;i++) {
U[i] = numeric.rep([right[i]-left[i]+1],0);
L[i] = numeric.rep([i-left[i]],0);
countL += i-left[i]+1;
countU += right[i]-i+1;
}
for(k=0;k<p;k++) { i = I[k]; U[i][J[k]-left[i]] = V[k]; }
for(i=0;i<m-1;i++) {
a = i-left[i];
Ui = U[i];
for(j=i+1;left[j]<=i && j<m;j++) {
b = i-left[j];
c = right[i]-i;
Uj = U[j];
alpha = Uj[b]/Ui[a];
if(alpha) {
for(k=1;k<=c;k++) { Uj[k+b] -= alpha*Ui[k+a]; }
L[j][i-left[j]] = alpha;
}
}
}
var Ui = [], Uj = [], Uv = [], Li = [], Lj = [], Lv = [];
var p,q,foo;
p=0; q=0;
for(i=0;i<m;i++) {
a = left[i];
b = right[i];
foo = U[i];
for(j=i;j<=b;j++) {
if(foo[j-a]) {
Ui[p] = i;
Uj[p] = j;
Uv[p] = foo[j-a];
p++;
}
}
foo = L[i];
for(j=a;j<i;j++) {
if(foo[j-a]) {
Li[q] = i;
Lj[q] = j;
Lv[q] = foo[j-a];
q++;
}
}
Li[q] = i;
Lj[q] = i;
Lv[q] = 1;
q++;
}
return {U:[Ui,Uj,Uv], L:[Li,Lj,Lv]};
};
numeric.cLUsolve = function LUsolve(lu,b) {
var L = lu.L, U = lu.U, ret = numeric.clone(b);
var Li = L[0], Lj = L[1], Lv = L[2];
var Ui = U[0], Uj = U[1], Uv = U[2];
var p = Ui.length, q = Li.length;
var m = ret.length,i,j,k;
k = 0;
for(i=0;i<m;i++) {
while(Lj[k] < i) {
ret[i] -= Lv[k]*ret[Lj[k]];
k++;
}
k++;
}
k = p-1;
for(i=m-1;i>=0;i--) {
while(Uj[k] > i) {
ret[i] -= Uv[k]*ret[Uj[k]];
k--;
}
ret[i] /= Uv[k];
k--;
}
return ret;
};
numeric.cgrid = function grid(n,shape) {
if(typeof n === "number") n = [n,n];
var ret = numeric.rep(n,-1);
var i,j,count;
if(typeof shape !== "function") {
switch(shape) {
case 'L':
shape = function(i,j) { return (i>=n[0]/2 || j<n[1]/2); }
break;
default:
shape = function(i,j) { return true; };
break;
}
}
count=0;
for(i=1;i<n[0]-1;i++) for(j=1;j<n[1]-1;j++)
if(shape(i,j)) {
ret[i][j] = count;
count++;
}
return ret;
}
numeric.cdelsq = function delsq(g) {
var dir = [[-1,0],[0,-1],[0,1],[1,0]];
var s = numeric.dim(g), m = s[0], n = s[1], i,j,k,p,q;
var Li = [], Lj = [], Lv = [];
for(i=1;i<m-1;i++) for(j=1;j<n-1;j++) {
if(g[i][j]<0) continue;
for(k=0;k<4;k++) {
p = i+dir[k][0];
q = j+dir[k][1];
if(g[p][q]<0) continue;
Li.push(g[i][j]);
Lj.push(g[p][q]);
Lv.push(-1);
}
Li.push(g[i][j]);
Lj.push(g[i][j]);
Lv.push(4);
}
return [Li,Lj,Lv];
}
numeric.cdotMV = function dotMV(A,x) {
var ret, Ai = A[0], Aj = A[1], Av = A[2],k,p=Ai.length,N;
N=0;
for(k=0;k<p;k++) { if(Ai[k]>N) N = Ai[k]; }
N++;
ret = numeric.rep([N],0);
for(k=0;k<p;k++) { ret[Ai[k]]+=Av[k]*x[Aj[k]]; }
return ret;
}
// 7. Splines
numeric.Spline = function Spline(x,yl,yr,kl,kr) { this.x = x; this.yl = yl; this.yr = yr; this.kl = kl; this.kr = kr; }
numeric.Spline.prototype._at = function _at(x1,p) {
var x = this.x;
var yl = this.yl;
var yr = this.yr;
var kl = this.kl;
var kr = this.kr;
var x1,a,b,t;
var add = numeric.add, sub = numeric.sub, mul = numeric.mul;
a = sub(mul(kl[p],x[p+1]-x[p]),sub(yr[p+1],yl[p]));
b = add(mul(kr[p+1],x[p]-x[p+1]),sub(yr[p+1],yl[p]));
t = (x1-x[p])/(x[p+1]-x[p]);
var s = t*(1-t);
return add(add(add(mul(1-t,yl[p]),mul(t,yr[p+1])),mul(a,s*(1-t))),mul(b,s*t));
}
numeric.Spline.prototype.at = function at(x0) {
if(typeof x0 === "number") {
var x = this.x;
var n = x.length;
var p,q,mid,floor = Math.floor,a,b,t;
p = 0;
q = n-1;
while(q-p>1) {
mid = floor((p+q)/2);
if(x[mid] <= x0) p = mid;
else q = mid;
}
return this._at(x0,p);
}
var n = x0.length, i, ret = Array(n);
for(i=n-1;i!==-1;--i) ret[i] = this.at(x0[i]);
return ret;
}
numeric.Spline.prototype.diff = function diff() {
var x = this.x;
var yl = this.yl;
var yr = this.yr;
var kl = this.kl;
var kr = this.kr;
var n = yl.length;
var i,dx,dy;
var zl = kl, zr = kr, pl = Array(n), pr = Array(n);
var add = numeric.add, mul = numeric.mul, div = numeric.div, sub = numeric.sub;
for(i=n-1;i!==-1;--i) {
dx = x[i+1]-x[i];
dy = sub(yr[i+1],yl[i]);
pl[i] = div(add(mul(dy, 6),mul(kl[i],-4*dx),mul(kr[i+1],-2*dx)),dx*dx);
pr[i+1] = div(add(mul(dy,-6),mul(kl[i], 2*dx),mul(kr[i+1], 4*dx)),dx*dx);
}
return new numeric.Spline(x,zl,zr,pl,pr);
}
numeric.Spline.prototype.roots = function roots() {
function sqr(x) { return x*x; }
function heval(y0,y1,k0,k1,x) {
var A = k0*2-(y1-y0);
var B = -k1*2+(y1-y0);
var t = (x+1)*0.5;
var s = t*(1-t);
return (1-t)*y0+t*y1+A*s*(1-t)+B*s*t;
}
var ret = [];
var x = this.x, yl = this.yl, yr = this.yr, kl = this.kl, kr = this.kr;
if(typeof yl[0] === "number") {
yl = [yl];
yr = [yr];
kl = [kl];
kr = [kr];
}
var m = yl.length,n=x.length-1,i,j,k,y,s,t;
var ai,bi,ci,di, ret = Array(m),ri,k0,k1,y0,y1,A,B,D,dx,cx,stops,z0,z1,zm,t0,t1,tm;
var sqrt = Math.sqrt;
for(i=0;i!==m;++i) {
ai = yl[i];
bi = yr[i];
ci = kl[i];
di = kr[i];
ri = [];
for(j=0;j!==n;j++) {
if(j>0 && bi[j]*ai[j]<0) ri.push(x[j]);
dx = (x[j+1]-x[j]);
cx = x[j];
y0 = ai[j];
y1 = bi[j+1];
k0 = ci[j]/dx;
k1 = di[j+1]/dx;
D = sqr(k0-k1+3*(y0-y1)) + 12*k1*y0;
A = k1+3*y0+2*k0-3*y1;
B = 3*(k1+k0+2*(y0-y1));
if(D<=0) {
z0 = A/B;
if(z0>x[j] && z0<x[j+1]) stops = [x[j],z0,x[j+1]];
else stops = [x[j],x[j+1]];
} else {
z0 = (A-sqrt(D))/B;
z1 = (A+sqrt(D))/B;
stops = [x[j]];
if(z0>x[j] && z0<x[j+1]) stops.push(z0);
if(z1>x[j] && z1<x[j+1]) stops.push(z1);
stops.push(x[j+1]);
}
t0 = stops[0];
z0 = this._at(t0,j);
for(k=0;k<stops.length-1;k++) {
t1 = stops[k+1];
z1 = this._at(t1,j);
if(z0 === 0) {
ri.push(t0);
t0 = t1;
z0 = z1;
continue;
}
if(z1 === 0 || z0*z1>0) {
t0 = t1;
z0 = z1;
continue;
}
var side = 0;
while(1) {
tm = (z0*t1-z1*t0)/(z0-z1);
if(tm <= t0 || tm >= t1) { break; }
zm = this._at(tm,j);
if(zm*z1>0) {
t1 = tm;
z1 = zm;
if(side === -1) z0*=0.5;
side = -1;
} else if(zm*z0>0) {
t0 = tm;
z0 = zm;
if(side === 1) z1*=0.5;
side = 1;
} else break;
}
ri.push(tm);
t0 = stops[k+1];
z0 = this._at(t0, j);
}
if(z1 === 0) ri.push(t1);
}
ret[i] = ri;
}
if(typeof this.yl[0] === "number") return ret[0];
return ret;
}
numeric.spline = function spline(x,y,k1,kn) {
var n = x.length, b = [], dx = [], dy = [];
var i;
var sub = numeric.sub,mul = numeric.mul,add = numeric.add;
for(i=n-2;i>=0;i--) { dx[i] = x[i+1]-x[i]; dy[i] = sub(y[i+1],y[i]); }
if(typeof k1 === "string" || typeof kn === "string") {
k1 = kn = "periodic";
}
// Build sparse tridiagonal system
var T = [[],[],[]];
switch(typeof k1) {
case "undefined":
b[0] = mul(3/(dx[0]*dx[0]),dy[0]);
T[0].push(0,0);
T[1].push(0,1);
T[2].push(2/dx[0],1/dx[0]);
break;
case "string":
b[0] = add(mul(3/(dx[n-2]*dx[n-2]),dy[n-2]),mul(3/(dx[0]*dx[0]),dy[0]));
T[0].push(0,0,0);
T[1].push(n-2,0,1);
T[2].push(1/dx[n-2],2/dx[n-2]+2/dx[0],1/dx[0]);
break;
default:
b[0] = k1;
T[0].push(0);
T[1].push(0);
T[2].push(1);
break;
}
for(i=1;i<n-1;i++) {
b[i] = add(mul(3/(dx[i-1]*dx[i-1]),dy[i-1]),mul(3/(dx[i]*dx[i]),dy[i]));
T[0].push(i,i,i);
T[1].push(i-1,i,i+1);
T[2].push(1/dx[i-1],2/dx[i-1]+2/dx[i],1/dx[i]);
}
switch(typeof kn) {
case "undefined":
b[n-1] = mul(3/(dx[n-2]*dx[n-2]),dy[n-2]);
T[0].push(n-1,n-1);
T[1].push(n-2,n-1);
T[2].push(1/dx[n-2],2/dx[n-2]);
break;
case "string":
T[1][T[1].length-1] = 0;
break;
default:
b[n-1] = kn;
T[0].push(n-1);
T[1].push(n-1);
T[2].push(1);
break;
}
if(typeof b[0] !== "number") b = numeric.transpose(b);
else b = [b];
var k = Array(b.length);
if(typeof k1 === "string") {
for(i=k.length-1;i!==-1;--i) {
k[i] = numeric.ccsLUPSolve(numeric.ccsLUP(numeric.ccsScatter(T)),b[i]);
k[i][n-1] = k[i][0];
}
} else {
for(i=k.length-1;i!==-1;--i) {
k[i] = numeric.cLUsolve(numeric.cLU(T),b[i]);
}
}
if(typeof y[0] === "number") k = k[0];
else k = numeric.transpose(k);
return new numeric.Spline(x,y,y,k,k);
}
// 8. FFT
numeric.fftpow2 = function fftpow2(x,y) {
var n = x.length;
if(n === 1) return;
var cos = Math.cos, sin = Math.sin, i,j;
var xe = Array(n/2), ye = Array(n/2), xo = Array(n/2), yo = Array(n/2);
j = n/2;
for(i=n-1;i!==-1;--i) {
--j;
xo[j] = x[i];
yo[j] = y[i];
--i;
xe[j] = x[i];
ye[j] = y[i];
}
fftpow2(xe,ye);
fftpow2(xo,yo);
j = n/2;
var t,k = (-6.2831853071795864769252867665590057683943387987502116419/n),ci,si;
for(i=n-1;i!==-1;--i) {
--j;
if(j === -1) j = n/2-1;
t = k*i;
ci = cos(t);
si = sin(t);
x[i] = xe[j] + ci*xo[j] - si*yo[j];
y[i] = ye[j] + ci*yo[j] + si*xo[j];
}
}
numeric._ifftpow2 = function _ifftpow2(x,y) {
var n = x.length;
if(n === 1) return;
var cos = Math.cos, sin = Math.sin, i,j;
var xe = Array(n/2), ye = Array(n/2), xo = Array(n/2), yo = Array(n/2);
j = n/2;
for(i=n-1;i!==-1;--i) {
--j;
xo[j] = x[i];
yo[j] = y[i];
--i;
xe[j] = x[i];
ye[j] = y[i];
}
_ifftpow2(xe,ye);
_ifftpow2(xo,yo);
j = n/2;
var t,k = (6.2831853071795864769252867665590057683943387987502116419/n),ci,si;
for(i=n-1;i!==-1;--i) {
--j;
if(j === -1) j = n/2-1;
t = k*i;
ci = cos(t);
si = sin(t);
x[i] = xe[j] + ci*xo[j] - si*yo[j];
y[i] = ye[j] + ci*yo[j] + si*xo[j];
}
}
numeric.ifftpow2 = function ifftpow2(x,y) {
numeric._ifftpow2(x,y);
numeric.diveq(x,x.length);
numeric.diveq(y,y.length);
}
numeric.convpow2 = function convpow2(ax,ay,bx,by) {
numeric.fftpow2(ax,ay);
numeric.fftpow2(bx,by);
var i,n = ax.length,axi,bxi,ayi,byi;
for(i=n-1;i!==-1;--i) {
axi = ax[i]; ayi = ay[i]; bxi = bx[i]; byi = by[i];
ax[i] = axi*bxi-ayi*byi;
ay[i] = axi*byi+ayi*bxi;
}
numeric.ifftpow2(ax,ay);
}
numeric.T.prototype.fft = function fft() {
var x = this.x, y = this.y;
var n = x.length, log = Math.log, log2 = log(2),
p = Math.ceil(log(2*n-1)/log2), m = Math.pow(2,p);
var cx = numeric.rep([m],0), cy = numeric.rep([m],0), cos = Math.cos, sin = Math.sin;
var k, c = (-3.141592653589793238462643383279502884197169399375105820/n),t;
var a = numeric.rep([m],0), b = numeric.rep([m],0),nhalf = Math.floor(n/2);
for(k=0;k<n;k++) a[k] = x[k];
if(typeof y !== "undefined") for(k=0;k<n;k++) b[k] = y[k];
cx[0] = 1;
for(k=1;k<=m/2;k++) {
t = c*k*k;
cx[k] = cos(t);
cy[k] = sin(t);
cx[m-k] = cos(t);
cy[m-k] = sin(t)
}
var X = new numeric.T(a,b), Y = new numeric.T(cx,cy);
X = X.mul(Y);
numeric.convpow2(X.x,X.y,numeric.clone(Y.x),numeric.neg(Y.y));
X = X.mul(Y);
X.x.length = n;
X.y.length = n;
return X;
}
numeric.T.prototype.ifft = function ifft() {
var x = this.x, y = this.y;
var n = x.length, log = Math.log, log2 = log(2),
p = Math.ceil(log(2*n-1)/log2), m = Math.pow(2,p);
var cx = numeric.rep([m],0), cy = numeric.rep([m],0), cos = Math.cos, sin = Math.sin;
var k, c = (3.141592653589793238462643383279502884197169399375105820/n),t;
var a = numeric.rep([m],0), b = numeric.rep([m],0),nhalf = Math.floor(n/2);
for(k=0;k<n;k++) a[k] = x[k];
if(typeof y !== "undefined") for(k=0;k<n;k++) b[k] = y[k];
cx[0] = 1;
for(k=1;k<=m/2;k++) {
t = c*k*k;
cx[k] = cos(t);
cy[k] = sin(t);
cx[m-k] = cos(t);
cy[m-k] = sin(t)
}
var X = new numeric.T(a,b), Y = new numeric.T(cx,cy);
X = X.mul(Y);
numeric.convpow2(X.x,X.y,numeric.clone(Y.x),numeric.neg(Y.y));
X = X.mul(Y);
X.x.length = n;
X.y.length = n;
return X.div(n);
}
//9. Unconstrained optimization
numeric.gradient = function gradient(f,x) {
var n = x.length;
var f0 = f(x);
if(isNaN(f0)) throw new Error('gradient: f(x) is a NaN!');
var max = Math.max;
var i,x0 = numeric.clone(x),f1,f2, J = Array(n);
var div = numeric.div, sub = numeric.sub,errest,roundoff,max = Math.max,eps = 1e-3,abs = Math.abs, min = Math.min;
var t0,t1,t2,it=0,d1,d2,N;
for(i=0;i<n;i++) {
var h = max(1e-6*f0,1e-8);
while(1) {
++it;
if(it>20) { throw new Error("Numerical gradient fails"); }
x0[i] = x[i]+h;
f1 = f(x0);
x0[i] = x[i]-h;
f2 = f(x0);
x0[i] = x[i];
if(isNaN(f1) || isNaN(f2)) { h/=16; continue; }
J[i] = (f1-f2)/(2*h);
t0 = x[i]-h;
t1 = x[i];
t2 = x[i]+h;
d1 = (f1-f0)/h;
d2 = (f0-f2)/h;
N = max(abs(J[i]),abs(f0),abs(f1),abs(f2),abs(t0),abs(t1),abs(t2),1e-8);
errest = min(max(abs(d1-J[i]),abs(d2-J[i]),abs(d1-d2))/N,h/N);
if(errest>eps) { h/=16; }
else break;
}
}
return J;
}
numeric.uncmin = function uncmin(f,x0,tol,gradient,maxit,callback,options) {
var grad = numeric.gradient;
if(typeof options === "undefined") { options = {}; }
if(typeof tol === "undefined") { tol = 1e-8; }
if(typeof gradient === "undefined") { gradient = function(x) { return grad(f,x); }; }
if(typeof maxit === "undefined") maxit = 1000;
x0 = numeric.clone(x0);
var n = x0.length;
var f0 = f(x0),f1,df0;
if(isNaN(f0)) throw new Error('uncmin: f(x0) is a NaN!');
var max = Math.max, norm2 = numeric.norm2;
tol = max(tol,numeric.epsilon);
var step,g0,g1,H1 = options.Hinv || numeric.identity(n);
var dot = numeric.dot, inv = numeric.inv, sub = numeric.sub, add = numeric.add, ten = numeric.tensor, div = numeric.div, mul = numeric.mul;
var all = numeric.all, isfinite = numeric.isFinite, neg = numeric.neg;
var it=0,i,s,x1,y,Hy,Hs,ys,i0,t,nstep,t1,t2;
var msg = "";
g0 = gradient(x0);
while(it<maxit) {
if(typeof callback === "function") { if(callback(it,x0,f0,g0,H1)) { msg = "Callback returned true"; break; } }
if(!all(isfinite(g0))) { msg = "Gradient has Infinity or NaN"; break; }
step = neg(dot(H1,g0));
if(!all(isfinite(step))) { msg = "Search direction has Infinity or NaN"; break; }
nstep = norm2(step);
if(nstep < tol) { msg="Newton step smaller than tol"; break; }
t = 1;
df0 = dot(g0,step);
// line search
x1 = x0;
while(it < maxit) {
if(t*nstep < tol) { break; }
s = mul(step,t);
x1 = add(x0,s);
f1 = f(x1);
if(f1-f0 >= 0.1*t*df0 || isNaN(f1)) {
t *= 0.5;
++it;
continue;
}
break;
}
if(t*nstep < tol) { msg = "Line search step size smaller than tol"; break; }
if(it === maxit) { msg = "maxit reached during line search"; break; }
g1 = gradient(x1);
y = sub(g1,g0);
ys = dot(y,s);
Hy = dot(H1,y);
H1 = sub(add(H1,
mul(
(ys+dot(y,Hy))/(ys*ys),
ten(s,s) )),
div(add(ten(Hy,s),ten(s,Hy)),ys));
x0 = x1;
f0 = f1;
g0 = g1;
++it;
}
return {solution: x0, f: f0, gradient: g0, invHessian: H1, iterations:it, message: msg};
}
// 10. Ode solver (Dormand-Prince)
numeric.Dopri = function Dopri(x,y,f,ymid,iterations,msg,events) {
this.x = x;
this.y = y;
this.f = f;
this.ymid = ymid;
this.iterations = iterations;
this.events = events;
this.message = msg;
}
numeric.Dopri.prototype._at = function _at(xi,j) {
function sqr(x) { return x*x; }
var sol = this;
var xs = sol.x;
var ys = sol.y;
var k1 = sol.f;
var ymid = sol.ymid;
var n = xs.length;
var x0,x1,xh,y0,y1,yh,xi;
var floor = Math.floor,h;
var c = 0.5;
var add = numeric.add, mul = numeric.mul,sub = numeric.sub, p,q,w;
x0 = xs[j];
x1 = xs[j+1];
y0 = ys[j];
y1 = ys[j+1];
h = x1-x0;
xh = x0+c*h;
yh = ymid[j];
p = sub(k1[j ],mul(y0,1/(x0-xh)+2/(x0-x1)));
q = sub(k1[j+1],mul(y1,1/(x1-xh)+2/(x1-x0)));
w = [sqr(xi - x1) * (xi - xh) / sqr(x0 - x1) / (x0 - xh),
sqr(xi - x0) * sqr(xi - x1) / sqr(x0 - xh) / sqr(x1 - xh),
sqr(xi - x0) * (xi - xh) / sqr(x1 - x0) / (x1 - xh),
(xi - x0) * sqr(xi - x1) * (xi - xh) / sqr(x0-x1) / (x0 - xh),
(xi - x1) * sqr(xi - x0) * (xi - xh) / sqr(x0-x1) / (x1 - xh)];
return add(add(add(add(mul(y0,w[0]),
mul(yh,w[1])),
mul(y1,w[2])),
mul( p,w[3])),
mul( q,w[4]));
}
numeric.Dopri.prototype.at = function at(x) {
var i,j,k,floor = Math.floor;
if(typeof x !== "number") {
var n = x.length, ret = Array(n);
for(i=n-1;i!==-1;--i) {
ret[i] = this.at(x[i]);
}
return ret;
}
var x0 = this.x;
i = 0; j = x0.length-1;
while(j-i>1) {
k = floor(0.5*(i+j));
if(x0[k] <= x) i = k;
else j = k;
}
return this._at(x,i);
}
numeric.dopri = function dopri(x0,x1,y0,f,tol,maxit,event) {
if(typeof tol === "undefined") { tol = 1e-6; }
if(typeof maxit === "undefined") { maxit = 1000; }
var xs = [x0], ys = [y0], k1 = [f(x0,y0)], k2,k3,k4,k5,k6,k7, ymid = [];
var A2 = 1/5;
var A3 = [3/40,9/40];
var A4 = [44/45,-56/15,32/9];
var A5 = [19372/6561,-25360/2187,64448/6561,-212/729];
var A6 = [9017/3168,-355/33,46732/5247,49/176,-5103/18656];
var b = [35/384,0,500/1113,125/192,-2187/6784,11/84];
var bm = [0.5*6025192743/30085553152,
0,
0.5*51252292925/65400821598,
0.5*-2691868925/45128329728,
0.5*187940372067/1594534317056,
0.5*-1776094331/19743644256,
0.5*11237099/235043384];
var c = [1/5,3/10,4/5,8/9,1,1];
var e = [-71/57600,0,71/16695,-71/1920,17253/339200,-22/525,1/40];
var i = 0,er,j;
var h = (x1-x0)/10;
var it = 0;
var add = numeric.add, mul = numeric.mul, y1,erinf;
var max = Math.max, min = Math.min, abs = Math.abs, norminf = numeric.norminf,pow = Math.pow;
var any = numeric.any, lt = numeric.lt, and = numeric.and, sub = numeric.sub;
var e0, e1, ev;
var ret = new numeric.Dopri(xs,ys,k1,ymid,-1,"");
if(typeof event === "function") e0 = event(x0,y0);
while(x0<x1 && it<maxit) {
++it;
if(x0+h>x1) h = x1-x0;
k2 = f(x0+c[0]*h, add(y0,mul( A2*h,k1[i])));
k3 = f(x0+c[1]*h, add(add(y0,mul(A3[0]*h,k1[i])),mul(A3[1]*h,k2)));
k4 = f(x0+c[2]*h, add(add(add(y0,mul(A4[0]*h,k1[i])),mul(A4[1]*h,k2)),mul(A4[2]*h,k3)));
k5 = f(x0+c[3]*h, add(add(add(add(y0,mul(A5[0]*h,k1[i])),mul(A5[1]*h,k2)),mul(A5[2]*h,k3)),mul(A5[3]*h,k4)));
k6 = f(x0+c[4]*h,add(add(add(add(add(y0,mul(A6[0]*h,k1[i])),mul(A6[1]*h,k2)),mul(A6[2]*h,k3)),mul(A6[3]*h,k4)),mul(A6[4]*h,k5)));
y1 = add(add(add(add(add(y0,mul(k1[i],h*b[0])),mul(k3,h*b[2])),mul(k4,h*b[3])),mul(k5,h*b[4])),mul(k6,h*b[5]));
k7 = f(x0+h,y1);
er = add(add(add(add(add(mul(k1[i],h*e[0]),mul(k3,h*e[2])),mul(k4,h*e[3])),mul(k5,h*e[4])),mul(k6,h*e[5])),mul(k7,h*e[6]));
if(typeof er === "number") erinf = abs(er);
else erinf = norminf(er);
if(erinf > tol) { // reject
h = 0.2*h*pow(tol/erinf,0.25);
if(x0+h === x0) {
ret.msg = "Step size became too small";
break;
}
continue;
}
ymid[i] = add(add(add(add(add(add(y0,
mul(k1[i],h*bm[0])),
mul(k3 ,h*bm[2])),
mul(k4 ,h*bm[3])),
mul(k5 ,h*bm[4])),
mul(k6 ,h*bm[5])),
mul(k7 ,h*bm[6]));
++i;
xs[i] = x0+h;
ys[i] = y1;
k1[i] = k7;
if(typeof event === "function") {
var yi,xl = x0,xr = x0+0.5*h,xi;
e1 = event(xr,ymid[i-1]);
ev = and(lt(e0,0),lt(0,e1));
if(!any(ev)) { xl = xr; xr = x0+h; e0 = e1; e1 = event(xr,y1); ev = and(lt(e0,0),lt(0,e1)); }
if(any(ev)) {
var xc, yc, en,ei;
var side=0, sl = 1.0, sr = 1.0;
while(1) {
if(typeof e0 === "number") xi = (sr*e1*xl-sl*e0*xr)/(sr*e1-sl*e0);
else {
xi = xr;
for(j=e0.length-1;j!==-1;--j) {
if(e0[j]<0 && e1[j]>0) xi = min(xi,(sr*e1[j]*xl-sl*e0[j]*xr)/(sr*e1[j]-sl*e0[j]));
}
}
if(xi <= xl || xi >= xr) break;
yi = ret._at(xi, i-1);
ei = event(xi,yi);
en = and(lt(e0,0),lt(0,ei));
if(any(en)) {
xr = xi;
e1 = ei;
ev = en;
sr = 1.0;
if(side === -1) sl *= 0.5;
else sl = 1.0;
side = -1;
} else {
xl = xi;
e0 = ei;
sl = 1.0;
if(side === 1) sr *= 0.5;
else sr = 1.0;
side = 1;
}
}
y1 = ret._at(0.5*(x0+xi),i-1);
ret.f[i] = f(xi,yi);
ret.x[i] = xi;
ret.y[i] = yi;
ret.ymid[i-1] = y1;
ret.events = ev;
ret.iterations = it;
return ret;
}
}
x0 += h;
y0 = y1;
e0 = e1;
h = min(0.8*h*pow(tol/erinf,0.25),4*h);
}
ret.iterations = it;
return ret;
}
// 11. Ax = b
numeric.LU = function(A, fast) {
fast = fast || false;
var abs = Math.abs;
var i, j, k, absAjk, Akk, Ak, Pk, Ai;
var max;
var n = A.length, n1 = n-1;
var P = new Array(n);
if(!fast) A = numeric.clone(A);
for (k = 0; k < n; ++k) {
Pk = k;
Ak = A[k];
max = abs(Ak[k]);
for (j = k + 1; j < n; ++j) {
absAjk = abs(A[j][k]);
if (max < absAjk) {
max = absAjk;
Pk = j;
}
}
P[k] = Pk;
if (Pk != k) {
A[k] = A[Pk];
A[Pk] = Ak;
Ak = A[k];
}
Akk = Ak[k];
for (i = k + 1; i < n; ++i) {
A[i][k] /= Akk;
}
for (i = k + 1; i < n; ++i) {
Ai = A[i];
for (j = k + 1; j < n1; ++j) {
Ai[j] -= Ai[k] * Ak[j];
++j;
Ai[j] -= Ai[k] * Ak[j];
}
if(j===n1) Ai[j] -= Ai[k] * Ak[j];
}
}
return {
LU: A,
P: P
};
}
numeric.LUsolve = function LUsolve(LUP, b) {
var i, j;
var LU = LUP.LU;
var n = LU.length;
var x = numeric.clone(b);
var P = LUP.P;
var Pi, LUi, LUii, tmp;
for (i=n-1;i!==-1;--i) x[i] = b[i];
for (i = 0; i < n; ++i) {
Pi = P[i];
if (P[i] !== i) {
tmp = x[i];
x[i] = x[Pi];
x[Pi] = tmp;
}
LUi = LU[i];
for (j = 0; j < i; ++j) {
x[i] -= x[j] * LUi[j];
}
}
for (i = n - 1; i >= 0; --i) {
LUi = LU[i];
for (j = i + 1; j < n; ++j) {
x[i] -= x[j] * LUi[j];
}
x[i] /= LUi[i];
}
return x;
}
numeric.solve = function solve(A,b,fast) { return numeric.LUsolve(numeric.LU(A,fast), b); }
// 12. Linear programming
numeric.echelonize = function echelonize(A) {
var s = numeric.dim(A), m = s[0], n = s[1];
var I = numeric.identity(m);
var P = Array(m);
var i,j,k,l,Ai,Ii,Z,a;
var abs = Math.abs;
var diveq = numeric.diveq;
A = numeric.clone(A);
for(i=0;i<m;++i) {
k = 0;
Ai = A[i];
Ii = I[i];
for(j=1;j<n;++j) if(abs(Ai[k])<abs(Ai[j])) k=j;
P[i] = k;
diveq(Ii,Ai[k]);
diveq(Ai,Ai[k]);
for(j=0;j<m;++j) if(j!==i) {
Z = A[j]; a = Z[k];
for(l=n-1;l!==-1;--l) Z[l] -= Ai[l]*a;
Z = I[j];
for(l=m-1;l!==-1;--l) Z[l] -= Ii[l]*a;
}
}
return {I:I, A:A, P:P};
}
numeric.__solveLP = function __solveLP(c,A,b,tol,maxit,x,flag) {
var sum = numeric.sum, log = numeric.log, mul = numeric.mul, sub = numeric.sub, dot = numeric.dot, div = numeric.div, add = numeric.add;
var m = c.length, n = b.length,y;
var unbounded = false, cb,i0=0;
var alpha = 1.0;
var f0,df0,AT = numeric.transpose(A), svd = numeric.svd,transpose = numeric.transpose,leq = numeric.leq, sqrt = Math.sqrt, abs = Math.abs;
var muleq = numeric.muleq;
var norm = numeric.norminf, any = numeric.any,min = Math.min;
var all = numeric.all, gt = numeric.gt;
var p = Array(m), A0 = Array(n),e=numeric.rep([n],1), H;
var solve = numeric.solve, z = sub(b,dot(A,x)),count;
var dotcc = dot(c,c);
var g;
for(count=i0;count<maxit;++count) {
var i,j,d;
for(i=n-1;i!==-1;--i) A0[i] = div(A[i],z[i]);
var A1 = transpose(A0);
for(i=m-1;i!==-1;--i) p[i] = (/*x[i]+*/sum(A1[i]));
alpha = 0.25*abs(dotcc/dot(c,p));
var a1 = 100*sqrt(dotcc/dot(p,p));
if(!isFinite(alpha) || alpha>a1) alpha = a1;
g = add(c,mul(alpha,p));
H = dot(A1,A0);
for(i=m-1;i!==-1;--i) H[i][i] += 1;
d = solve(H,div(g,alpha),true);
var t0 = div(z,dot(A,d));
var t = 1.0;
for(i=n-1;i!==-1;--i) if(t0[i]<0) t = min(t,-0.999*t0[i]);
y = sub(x,mul(d,t));
z = sub(b,dot(A,y));
if(!all(gt(z,0))) return { solution: x, message: "", iterations: count };
x = y;
if(alpha<tol) return { solution: y, message: "", iterations: count };
if(flag) {
var s = dot(c,g), Ag = dot(A,g);
unbounded = true;
for(i=n-1;i!==-1;--i) if(s*Ag[i]<0) { unbounded = false; break; }
} else {
if(x[m-1]>=0) unbounded = false;
else unbounded = true;
}
if(unbounded) return { solution: y, message: "Unbounded", iterations: count };
}
return { solution: x, message: "maximum iteration count exceeded", iterations:count };
}
numeric._solveLP = function _solveLP(c,A,b,tol,maxit) {
var m = c.length, n = b.length,y;
var sum = numeric.sum, log = numeric.log, mul = numeric.mul, sub = numeric.sub, dot = numeric.dot, div = numeric.div, add = numeric.add;
var c0 = numeric.rep([m],0).concat([1]);
var J = numeric.rep([n,1],-1);
var A0 = numeric.blockMatrix([[A , J ]]);
var b0 = b;
var y = numeric.rep([m],0).concat(Math.max(0,numeric.sup(numeric.neg(b)))+1);
var x0 = numeric.__solveLP(c0,A0,b0,tol,maxit,y,false);
var x = numeric.clone(x0.solution);
x.length = m;
var foo = numeric.inf(sub(b,dot(A,x)));
if(foo<0) { return { solution: NaN, message: "Infeasible", iterations: x0.iterations }; }
var ret = numeric.__solveLP(c, A, b, tol, maxit-x0.iterations, x, true);
ret.iterations += x0.iterations;
return ret;
};
numeric.solveLP = function solveLP(c,A,b,Aeq,beq,tol,maxit) {
if(typeof maxit === "undefined") maxit = 1000;
if(typeof tol === "undefined") tol = numeric.epsilon;
if(typeof Aeq === "undefined") return numeric._solveLP(c,A,b,tol,maxit);
var m = Aeq.length, n = Aeq[0].length, o = A.length;
var B = numeric.echelonize(Aeq);
var flags = numeric.rep([n],0);
var P = B.P;
var Q = [];
var i;
for(i=P.length-1;i!==-1;--i) flags[P[i]] = 1;
for(i=n-1;i!==-1;--i) if(flags[i]===0) Q.push(i);
var g = numeric.getRange;
var I = numeric.linspace(0,m-1), J = numeric.linspace(0,o-1);
var Aeq2 = g(Aeq,I,Q), A1 = g(A,J,P), A2 = g(A,J,Q), dot = numeric.dot, sub = numeric.sub;
var A3 = dot(A1,B.I);
var A4 = sub(A2,dot(A3,Aeq2)), b4 = sub(b,dot(A3,beq));
var c1 = Array(P.length), c2 = Array(Q.length);
for(i=P.length-1;i!==-1;--i) c1[i] = c[P[i]];
for(i=Q.length-1;i!==-1;--i) c2[i] = c[Q[i]];
var c4 = sub(c2,dot(c1,dot(B.I,Aeq2)));
var S = numeric._solveLP(c4,A4,b4,tol,maxit);
var x2 = S.solution;
if(x2!==x2) return S;
var x1 = dot(B.I,sub(beq,dot(Aeq2,x2)));
var x = Array(c.length);
for(i=P.length-1;i!==-1;--i) x[P[i]] = x1[i];
for(i=Q.length-1;i!==-1;--i) x[Q[i]] = x2[i];
return { solution: x, message:S.message, iterations: S.iterations };
}
numeric.MPStoLP = function MPStoLP(MPS) {
if(MPS instanceof String) { MPS.split('\n'); }
var state = 0;
var states = ['Initial state','NAME','ROWS','COLUMNS','RHS','BOUNDS','ENDATA'];
var n = MPS.length;
var i,j,z,N=0,rows = {}, sign = [], rl = 0, vars = {}, nv = 0;
var name;
var c = [], A = [], b = [];
function err(e) { throw new Error('MPStoLP: '+e+'\nLine '+i+': '+MPS[i]+'\nCurrent state: '+states[state]+'\n'); }
for(i=0;i<n;++i) {
z = MPS[i];
var w0 = z.match(/\S*/g);
var w = [];
for(j=0;j<w0.length;++j) if(w0[j]!=="") w.push(w0[j]);
if(w.length === 0) continue;
for(j=0;j<states.length;++j) if(z.substr(0,states[j].length) === states[j]) break;
if(j<states.length) {
state = j;
if(j===1) { name = w[1]; }
if(j===6) return { name:name, c:c, A:numeric.transpose(A), b:b, rows:rows, vars:vars };
continue;
}
switch(state) {
case 0: case 1: err('Unexpected line');
case 2:
switch(w[0]) {
case 'N': if(N===0) N = w[1]; else err('Two or more N rows'); break;
case 'L': rows[w[1]] = rl; sign[rl] = 1; b[rl] = 0; ++rl; break;
case 'G': rows[w[1]] = rl; sign[rl] = -1;b[rl] = 0; ++rl; break;
case 'E': rows[w[1]] = rl; sign[rl] = 0;b[rl] = 0; ++rl; break;
default: err('Parse error '+numeric.prettyPrint(w));
}
break;
case 3:
if(!vars.hasOwnProperty(w[0])) { vars[w[0]] = nv; c[nv] = 0; A[nv] = numeric.rep([rl],0); ++nv; }
var p = vars[w[0]];
for(j=1;j<w.length;j+=2) {
if(w[j] === N) { c[p] = parseFloat(w[j+1]); continue; }
var q = rows[w[j]];
A[p][q] = (sign[q]<0?-1:1)*parseFloat(w[j+1]);
}
break;
case 4:
for(j=1;j<w.length;j+=2) b[rows[w[j]]] = (sign[rows[w[j]]]<0?-1:1)*parseFloat(w[j+1]);
break;
case 5: /*FIXME*/ break;
case 6: err('Internal error');
}
}
err('Reached end of file without ENDATA');
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment