Skip to content

Instantly share code, notes, and snippets.

@uwcc

uwcc/.block

Last active October 19, 2020 23:35
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 41 You must be signed in to fork a gist
  • Save uwcc/9eb15a4c2898f83906cf9964eead8bea to your computer and use it in GitHub Desktop.
Save uwcc/9eb15a4c2898f83906cf9964eead8bea to your computer and use it in GitHub Desktop.
Download ZIP
2020 MDDN342 Assignment 3: Face Mappings
Raw
.block
license: mit
Raw
README.md

2020 MDDN342 Assignment 3: Face Mappings

REPLACE ALL TEXT IN THIS FILE

This README should eventually explain your different paramaterised faces.

Raw
face.js
/*
* FaceMap class - holds all informaiton about one mapped
* face and is able to draw itself.
*/
// remove this or set to false to enable full program (load will be slower)
var DEBUG_MODE = true;
// this can be used to set the number of sliders to show
var NUM_SLIDERS = 3;
// other variables can be in here too
// here's some examples for colors used
const bg_color = [225, 206, 187];
const fg_color = [151, 102, 52];
const stroke_color = [95, 52, 8];
// example of a global function
// given a segment, this returns the average point [x, y]
function segment_average(segment) {
let sum_x = 0;
let sum_y = 0;
let s_len = segment.length;
for (let i=0; i<s_len; i++) {
sum_x = sum_x + segment[i][0];
sum_y = sum_y + segment[i][1];
}
return [sum_x / s_len , sum_y / s_len ];
}
// This where you define your own face object
function Face() {
// these are state variables for a face
// (your variables should be different!)
this.num_eyes = 2; // can be either 1 (cyclops) or 2 (two eyes)
this.eye_shift = -1; // range is -10 to 10
this.mouth_value = 1; // range is 0.5 to 8
// example of a function *inside* the face object.
// this draws a segment, and do_loop will connect the ends if true
this.draw_segment = function(segment, do_loop) {
for(let i=0; i<segment.length; i++) {
let px = segment[i][0];
let py = segment[i][1];
ellipse(px, py, 0.1);
if(i < segment.length - 1) {
let nx = segment[i+1][0];
let ny = segment[i+1][1];
line(px, py, nx, ny);
}
else if(do_loop) {
let nx = segment[0][0];
let ny = segment[0][1];
line(px, py, nx, ny);
}
}
};
/*
* Draw the face with position lists that include:
* chin, right_eye, left_eye, right_eyebrow, left_eyebrow
* bottom_lip, top_lip, nose_tip, nose_bridge,
*/
this.draw = function(positions) {
// head
stroke(stroke_color);
fill(fg_color);
ellipse(0, 0, 3, 4);
noStroke();
// mouth
fill(bg_color);
ellipse(0, 0.64, 1.36, 0.25 * this.mouth_value);
// eyebrows
fill(0);
stroke(0);
strokeWeight(0.08);
this.draw_segment(positions.left_eyebrow);
this.draw_segment(positions.right_eyebrow);
// draw segments of face using points
fill(128);
stroke(128);
this.draw_segment(positions.chin);
fill(100, 0, 100);
stroke(100, 0, 100);
this.draw_segment(positions.nose_bridge);
this.draw_segment(positions.nose_tip);
strokeWeight(0.03);
fill(200, 0, 0);
stroke(200, 0, 0);
this.draw_segment(positions.top_lip);
this.draw_segment(positions.bottom_lip);
fill(255);
stroke(255);
let left_eye_pos = segment_average(positions.left_eye);
let right_eye_pos = segment_average(positions.right_eye);
// eyes
noStroke();
let curEyeShift = 0.04 * this.eye_shift;
if(this.num_eyes == 2) {
fill(bg_color);
ellipse(left_eye_pos[0], left_eye_pos[1], 0.45, 0.27);
ellipse(right_eye_pos[0], right_eye_pos[1], 0.45, 0.27);
fill(fg_color);
ellipse(left_eye_pos[0] + curEyeShift, left_eye_pos[1], 0.18);
ellipse(right_eye_pos[0] + curEyeShift, right_eye_pos[1], 0.18);
}
else {
let eyePosX = (left_eye_pos[0] + right_eye_pos[0]) / 2;
let eyePosY = (left_eye_pos[1] + right_eye_pos[1]) / 2;
fill(bg_color);
ellipse(eyePosX, eyePosY, 0.45, 0.27);
fill(fg_color);
ellipse(eyePosX - 0.1 + curEyeShift, eyePosY, 0.18);
}
}
/* set internal properties based on list numbers 0-100 */
this.setProperties = function(settings) {
this.num_eyes = int(map(settings[0], 0, 100, 1, 2));
this.eye_shift = map(settings[1], 0, 100, -2, 2);
this.mouth_value = map(settings[2], 0, 100, 0.5, 8);
}
/* get internal properties as list of numbers 0-100 */
this.getProperties = function() {
let settings = new Array(3);
settings[0] = map(this.num_eyes, 1, 2, 0, 100);
settings[1] = map(this.eye_shift, -2, 2, 0, 100);
settings[2] = map(this.mouth_value, 0.5, 8, 0, 100);
return settings;
}
}
Raw
index.html
<head>
<script src="https://cdnjs.cloudflare.com/ajax/libs/p5.js/1.1.9/p5.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/seedrandom/2.4.3/seedrandom.min.js"></script>
<script src="https://d3js.org/d3-random.v1.min.js"></script>
<script src="z_clmtrackr.js"></script>
<script src="z_model_pca_20_svm.js"></script>
<script language="javascript" type="text/javascript" src="z_purview_helper.js"></script>
<script language="javascript" type="text/javascript" src="z_focused_random.js"></script>
<script language="javascript" type="text/javascript" src="z_face_system.js"></script>
<script language="javascript" type="text/javascript" src="z_kdTree.js"></script>
<script language="javascript" type="text/javascript" src="z_face-api.js"></script>
<script language="javascript" type="text/javascript" src="z_training_images.js"></script>
<script language="javascript" type="text/javascript" src="z_testing_images.js"></script>
<script language="javascript" type="text/javascript" src="face.js"></script>
<script language="javascript" type="text/javascript" src="system_runner.js"></script>
<style>
body { padding: 0; margin: 0; }
.inner { position: absolute; }
#controls {
font: 300 12px "Helvetica Neue";
padding: 5;
margin: 5;
background: #f0f0f0;
opacity: 0.0;
-webkit-transition: opacity 0.2s ease;
-moz-transition: opacity 0.2s ease;
-o-transition: opacity 0.2s ease;
-ms-transition: opacity 0.2s ease;
}
#controls:hover { opacity: 0.9; }
</style>
</head>
<body style="background-color:white">
<div class="outer">
<div class="inner" id="controls" height="500">
<table>
<tr>
<td>setting 1</td>
<td id="slider1Container"></td>
</tr>
<tr>
<td>setting 2</td>
<td id="slider2Container"></td>
</tr>
<tr>
<td>setting 3</td>
<td id="slider3Container"></td>
</tr>
<tr>
<td>setting 4</td>
<td id="slider4Container"></td>
</tr>
<tr>
<td>setting 5</td>
<td id="slider5Container"></td>
</tr>
<tr>
<td>setting 6</td>
<td id="slider6Container"></td>
</tr>
<tr>
<td>setting 7</td>
<td id="slider7Container"></td>
</tr>
<tr>
<td>setting 8</td>
<td id="slider8Container"></td>
</tr>
<!-- YOU CAN ADD MORE SLIDERS HERE, LIKE THIS
<tr>
<td>setting 9</td>
<td id="slider9Container"></td>
</tr>
<tr>
<td>setting 10</td>
<td id="slider10Container"></td>
</tr>
<tr>
<td>setting 11</td>
<td id="slider11Container"></td>
</tr>
<tr>
<td>setting 12</td>
<td id="slider12Container"></td>
</tr>
... AND KEEP GOING
-->
<tr>
</tr>
<tr>
<td>show target</td>
<td id="sliderTintContainer"></td>
</tr>
<tr>
<td>Draw function</td>
<td id="selector1Container"></td>
</tr>
<tr>
<td>Face Draw</td>
<td id="checkbox1Container"></td>
</tr>
<tr>
<td>Face Targets</td>
<td id="checkbox2Container"></td>
</tr>
<tr>
<td>Face Points</td>
<td id="checkbox3Container"></td>
</tr>
<tr>
<td></td>
<td id="button1Container"></td>
</tr>
<tr>
<td></td>
<td id="button2Container"></td>
</tr>
<tr>
<td></td>
<td id="button3Container"></td>
</tr>
<tr>
<td></td>
<td id="button4Container"></td>
</tr>
</table>
</div>
<div>
<div id="canvasContainer"></div>
<a href="face.js">face code</a><br>
<a href="sketch.html">sketches</a>
</div>
</div>
<pre>
<p id="output">
</p>
</pre>
</body>
Raw
sample_images.json
[
"williams.jpg",
"oscar_selfie.jpg"
]
Raw
sketch.html
<head>
<style> body {padding: 0; margin: 0;} </style>
</head>
<body style="background-color:white">
<img src="same_pose.jpg" width="960" height="500"/><br>
Same Pose
<hr>
<img src="same_subject.jpg" width="960" height="500"/><br>
Same Subject
<p>
<a href="index.html">program</a>
</body>
Raw
system_runner.js
// recent work from my @VicUniWgtn 3rd year media design students that combines @nulhom's dlib face library with @p5xjs
var canvasWidth = 960;
var canvasHeight = 500;
var button;
var curRandomSeed;
var mainFace;
var faceImages = [];
var curFaceIndex = 0;
var curTrainIndex = 0;
var curValidIndex = 0;
var main_canvas;
var video_buffer;
var faceSelector;
var drawFaceCheckbox;
var faceGuideCheckbox;
var facePointsCheckbox;
var sliders = [];
var sliderTint;
var trainDataKeys = []
var trainValues = {}
var validDataKeys = []
var faceMapping = null;
let model_loaded = false;
let faces_processing = false;
let faces_processed = false;
var sample_images;
var selfieData = []
var video_capture = null;
if (typeof DEBUG_MODE === 'undefined' || DEBUG_MODE === null) {
var DEBUG_MODE = false;
}
if (typeof NUM_SLIDERS === 'undefined' || NUM_SLIDERS === null) {
var NUM_SLIDERS = 12;
}
async function preload () {
sample_images = loadJSON('sample_images.json')
trainValues = loadJSON('training_values.json');
if (!DEBUG_MODE) {
await faceapi.loadSsdMobilenetv1Model("./");
await faceapi.loadFaceLandmarkModel("./");
await faceapi.loadFaceRecognitionModel("./");
}
model_loaded = true;
}
var allEmbeddingsTree;
var allEmbeddings = [];
var embeddingDimensions;
var curNeighbors = [];
function squaredDistance(a, b) {
var sum = 0;
var length = 128;
for(var i=0; i<128; i++) {
var diff = a[i] - b[i];
sum += diff * diff;
}
// print(a.length,diff);
// print(sum, a==b);
return sum;
}
var haveStarted = false;
function setup () {
let keys = Object.keys(sample_images);
for (let i=0; i<keys.length; i++) {
let obj = {};
obj.url = sample_images[keys[i]];
selfieData.push(obj);
}
// create the drawing canvas, save the canvas element
main_canvas = createCanvas(canvasWidth, canvasHeight);
main_canvas.parent('canvasContainer');
curRandomSeed = int(focusedRandom(0, 100));
mainFace = new Face();
littleFace = new Face();
for(var i=0; i<selfieData.length; i++) {
var data = selfieData[i];
data.image = loadImage(data.url);
}
trainDataKeys = Object.keys(trainData);
for(var i=0; i<trainDataKeys.length; i++) {
var curKey = trainDataKeys[i];
var data = trainData[curKey];
var curEmbedding = data.embedding[0];
if(curEmbedding.length == 128) {
curEmbedding.push(curKey);
allEmbeddings.push(curEmbedding);
}
else {
print("rejected embedding ", curEmbedding.length, curEmbedding);
}
data.image = loadImage(data.url);
}
validDataKeys = Object.keys(validData);
for(var i=0; i<validDataKeys.length; i++) {
var curKey = validDataKeys[i];
var data = validData[curKey];
data.image = loadImage(data.url);
}
// print("Length: ", allEmbeddings.length);
// setup k-d tree
var N = allEmbeddings[0].length - 1;
embeddingDimensions = Array.apply(null, {length: N}).map(Number.call, Number);
// print(embeddingDimensions)
allEmbeddingsTree = new kdTree(allEmbeddings, squaredDistance, embeddingDimensions);
// print(allEmbeddingsTree)
faceSelector = createSelect();
if(DEBUG_MODE) {
faceSelector.option('Train');
faceSelector.value('Train');
}
else {
faceSelector.option('Faces');
faceSelector.option('Train');
faceSelector.option('NearestNeighbors');
faceSelector.option('TrainingQuiz');
faceSelector.option('InterpolationQuiz');
faceSelector.option('Video');
faceSelector.value('Faces');
}
faceSelector.parent('selector1Container');
/* create the sliders */
for(i=1; i<=NUM_SLIDERS; i++) {
var slider = createSlider(0, 100, 50);
var parentStr = 'slider' + i + 'Container';
slider.parent(parentStr);
sliders.push(slider);
}
drawFaceCheckbox = createCheckbox('', true);
drawFaceCheckbox.parent('checkbox1Container');
faceGuideCheckbox = createCheckbox('', false);
faceGuideCheckbox.parent('checkbox2Container');
facePointsCheckbox = createCheckbox('', false);
facePointsCheckbox.parent('checkbox3Container');
if(!DEBUG_MODE) {
sliderTint = createSlider(0, 100, 10);
sliderTint.parent("sliderTintContainer");
var interpButton = createButton('interpolate');
interpButton.mousePressed(interpolateCurrent);
interpButton.parent('button1Container');
}
/* and the buttons */
var loadButton = createButton('load');
loadButton.mousePressed(loadCurrentSettings);
loadButton.parent('button1Container');
var saveButton = createButton('save');
saveButton.mousePressed(saveCurrentSettings);
saveButton.parent('button2Container');
var getValuesButton = createButton('get values');
getValuesButton.mousePressed(getSingleJson);
getValuesButton.parent('button3Container');
var getAllValuesButton = createButton('get all values');
getAllValuesButton.mousePressed(getAllJson);
getAllValuesButton.parent('button4Container');
updateSlidersForTraining();
// rotation in degrees
angleMode(DEGREES);
background(255);
fill(0);
textSize(50);
textAlign(CENTER);
text("(waiting for models to load...)", width/2, height/2);
haveStarted = true;
}
function saveCurrentSettings() {
var curKey = trainDataKeys[curTrainIndex];
obj = mainFace.getProperties();
trainValues[curKey] = obj;
// for(var i=0; i<obj.length; i++) {
// trainData[curKey][i] = obj[i];
// }
var text = select('#output');
text.html("Storing values for " + curKey);
// getAllJson();
}
function getSingleJson() {
obj = mainFace.getProperties();
var text = select('#output');
var json = JSON.stringify(obj, null, 2);
text.html(json)
}
function getAllJson() {
obj = trainValues;
var text = select('#output');
var json = JSON.stringify(obj, null, 2);
// alert(json);
text.html(json)
}
// global variables for colors
var bg_color1 = [50, 50, 50];
var lastSwapTime = 0;
var millisPerSwap = 5000;
function changeRandomSeed() {
curRandomSeed = curRandomSeed + 1;
lastSwapTime = millis();
}
function mouseClicked() {
// changeRandomSeed();
}
var quiz_done = true;
var guessed_answer = 0;
function num_dist(e1, e2) {
let dist = 0;
for(let i=0; i<e1.length; i++) {
print(i, e1[i], e2[i]);
dist = dist + Math.abs(e1[i] - e2[i]);
}
return dist;
}
var processing_vid_face = false;
var lastProcessedVidFace = null;
async function draw () {
if (!model_loaded) {
return;
}
if (!DEBUG_MODE) {
if (!faces_processing) {
faces_processing = true;
background(255);
fill(0);
textSize(50);
textAlign(CENTER);
text("(processing faces...)", width/2, height/2);
for(var i=0; i<selfieData.length; i++) {
var data = selfieData[i];
let fullFaceDescriptions = await faceapi.detectAllFaces(data.image.canvas).withFaceLandmarks().withFaceDescriptors();
data.landmarks = get_landmarks(fullFaceDescriptions);
data.embedding = get_latents(fullFaceDescriptions);
}
// print("Some distances")
// var data = selfieData[0];
// for(var i=0; i<data.landmarks.length; i++) {
// for(var j=0; j<data.landmarks.length; j++) {
// print("dist ", i, "old to ", j, " new -> ", num_dist(data.embedding[i], data.latents[j]));
// }
// }
faces_processed = true;
return;
}
if(!faces_processed) {
return;
}
}
var mode = faceSelector.value();
if(millis() > lastSwapTime + millisPerSwap) {
lastSwapTime = millis();
// changeRandomSeed();
}
resetFocusedRandom(curRandomSeed);
noStroke();
var textDisplay = "unknown";
var params = [];
for(var i=0; i<NUM_SLIDERS; i++) {
params.push(sliders[i].value());
}
if (mode == 'Faces' || mode == 'FacePair' || mode == 'Train' || mode == 'Video') {
var do_train = (mode == 'Train');
var is_face = (mode == 'Faces');
var is_video = (mode == 'Video');
var show_face_guide = faceGuideCheckbox.checked();
var show_face_points = facePointsCheckbox.checked();
var do_draw_face = drawFaceCheckbox.checked();
if(do_train) {
// var keys = Object.keys(trainData);
var curKey = trainDataKeys[curTrainIndex];
var data = trainData[curKey];
}
else {
var data = selfieData[curFaceIndex];
}
const v_w = 640;
const v_h = 480;
const v_w_p = v_w / this._pixelDensity;
const v_h_p = v_h / this._pixelDensity;
// setup "video image" if in video mode
if (is_video) {
if (video_capture == null) {
video_capture = createCapture(VIDEO);
video_capture.size(canvasWidth, canvasHeight);
video_buffer = createGraphics(v_w, v_h);
return;
}
// print(processing_vid_face);
if(processing_vid_face == false) {
video_buffer.image(video_capture, 0, 0, v_w_p, v_h_p);
// print("grabbing video");
processing_vid_face = true;
lastProcessedVidFace = await faceapi.detectAllFaces(video_buffer.canvas).withFaceLandmarks().withFaceDescriptors();
processing_vid_face = false;
// print("grabbed video");
return;
}
if (lastProcessedVidFace == null) {
background(255);
fill(0);
textSize(50);
textAlign(CENTER);
text("(setting up camera...)", width/2, height/2);
return
}
data = {};
data["image"] = video_buffer;
data["landmarks"] = get_landmarks(lastProcessedVidFace);
data["embedding"] = get_latents(lastProcessedVidFace);
// print("Found faces: ", data.landmarks.length)
}
// we are not bailing, draw background
background(bg_color1);
// Displays the image at its actual size at point (0,0)
var img = data.image
if (is_face) {
x2 = 0;
y1 = 0;
x1 = 0;
var im_w = canvasWidth;
var im_h = canvasHeight;
var rect_w = canvasWidth;
var rect_h = canvasHeight;
image(img, x2, y1, im_w, im_h);
}
else if(is_video) {
// let vw = video_capture.elt.videoWidth;
// let vh = video_capture.elt.videoHeight;
// x2 = 0;
// y1 = 0;
// x1 = 0;
x1 = (canvasWidth - v_w) / 2.0;
x2 = x1;
y1 = (canvasHeight - v_h) / 2.0;
var im_w = v_w;
var im_h = v_h;
var rect_w = v_w;
var rect_h = v_h;
// print(im_w, im_h);
image(img, x2, y1, im_w, im_h, 0, 0, v_w_p, v_h_p);
}
else {
var x1 = (width/4-400/2);
var x2 = (3*width/4-400/2);
var y1 = (height/2-400/2);
var rect_w = 400;
var rect_h = 400;
var im_w = 400;
var im_h = 400;
image(img, x1, y1, 400, 400);
}
if(do_train) {
if (curKey in trainValues) {
fill(0, 200, 0);
}
else {
fill(200, 0, 0);
}
ellipse(x1+400/2, y1+400+15, 10, 10);
}
if(!DEBUG_MODE) {
noStroke();
var curSliderTintValue = sliderTint.value();
var overlayAlpha = map(curSliderTintValue, 0, 100, 255, 0);
fill(bg_color1[0], bg_color1[1], bg_color1[2], overlayAlpha);
rect(x2, y1, rect_w, rect_h);
}
stroke(0);
fill(255);
for(var i=0; i<data.landmarks.length; i++) {
// get array of face marker positions [x, y] format
var positions = data.landmarks[i];
var shifted_positions = JSON.parse(JSON.stringify(positions))
var data_mean = [0.0, 0.0];
var data_scale = 1.0;
var data_angle = 0.0;
if ('transform' in positions) {
data_mean = positions.transform.center;
data_scale = positions.transform.scale;
data_angle = positions.transform.angle;
delete shifted_positions.transform
}
var scale_x = im_w / img.width;
var scale_y = im_h / img.height;
push();
translate(x1, y1)
translate(scale_x*data_mean[0], scale_y*data_mean[1]);
scale(scale_x*data_scale, scale_y*data_scale);
rotate(degrees(data_angle));
stroke(0);
fill(255);
strokeWeight(1/data_scale);
Object.keys(positions).forEach(function(key) {
if (key=='transform') {
return;
}
var curSection = positions[key];
var shiftedSection = shifted_positions[key];
for (var i=0; i<curSection.length; i++) {
var cur_x = curSection[i][0];
var cur_y = curSection[i][1];
if (show_face_points) {
ellipse(cur_x, cur_y, 5/data_scale, 5/data_scale);
}
// get ready for drawing the face
shiftedSection[i][0] = cur_x;
shiftedSection[i][1] = cur_y;
}
});
noFill();
if(show_face_guide) {
stroke(0, 0, 255);
ellipse(0, 0, 4, 4);
line(0, -2, 0, 2);
line(-2, 0, 2, 0);
}
// ellipse(x1+data_mean[0], y1+data_mean[1], 4*data_scale, 4*data_scale);
// line(x1+data_mean[0], y1+data_mean[1]-2*data_scale, x1+data_mean[0], y1+data_mean[1]+2*data_scale);
pop();
var settings = params;
if (Object.keys(trainValues).length > 0) {
// NOT NOW
if ((typeof data.embedding !== 'undefined') &&
(data.embedding != null) &&
(data.embedding.length > i) &&
(data.embedding[i] != null) &&
(typeof data.embedding[i].length !== 'undefined') &&
(data.embedding[i].length == 128)) {
// print("Using embedding ", i)
var curEmbedding = data.embedding[i];
results = getAverageSettingsFrom(curEmbedding);
settings = results[0];
}
}
push();
translate(x2, y1)
translate(scale_x*data_mean[0], scale_y*data_mean[1]);
scale(scale_x*data_scale, scale_y*data_scale);
rotate(degrees(data_angle));
strokeWeight(1/data_scale);
mainFace.setProperties(settings);
if (do_draw_face) {
mainFace.draw(shifted_positions);
}
pop();
}
if(do_train) {
textDisplay = "Train: " + curKey;
}
else {
textDisplay = "";
}
}
else if (mode == 'NearestNeighbors') {
background(bg_color1);
// var keys = Object.keys(trainData);
var curKey = trainDataKeys[curTrainIndex];
var data = trainData[curKey];
// Displays the image at its actual size at point (0,0)
var img = data.image
var x1 = (width/4-250/2);
var y1 = (height/3-250/2);
image(img, x1, y1, 250, 250);
if (curKey in trainValues) {
fill(0, 200, 0);
}
else {
fill(200, 0, 0);
}
ellipse(x1+250/2, y1+250+15, 10, 10);
var y2 = (3*height/4-80/2);
for(var i=0; i<4; i++) {
// var keys = Object.keys(trainData);
var curKey = curNeighbors[i];
var nearData = trainData[curKey];
// Displays the image at its actual size at point (0,0)
var img = nearData.image
var x2 = (width/4 - 200 + i*100);
image(img, x2, y2, 80, 80);
}
for(var i=0; i<1; i++) {
// get array of face marker positions [x, y] format
var positions = data.landmarks[i];
var shifted_positions = JSON.parse(JSON.stringify(positions))
var data_mean = [0.0, 0.0];
var data_scale = 1.0;
var data_angle = 0.0;
if ('transform' in positions) {
data_mean = positions.transform.center;
data_scale = positions.transform.scale;
data_angle = positions.transform.angle;
delete shifted_positions.transform
}
var scale_x = 400.0 / img.width;
var scale_y = 400.0 / img.height;
Object.keys(positions).forEach(function(key) {
if (key=='transform') {
return;
}
var curSection = positions[key];
var shiftedSection = shifted_positions[key];
for (var i=0; i<curSection.length; i++) {
var cur_x = curSection[i][0];
var cur_y = curSection[i][1];
// get ready for drawing the face
shiftedSection[i][0] = cur_x;
shiftedSection[i][1] = cur_y;
}
});
var scale_x = 250.0 / img.width;
var scale_y = 250.0 / img.height;
var x2 = (3*width/4-250/2);
push();
translate(x2, y1);
translate(scale_x*data_mean[0], scale_y*data_mean[1]);
scale(scale_x*data_scale, scale_y*data_scale);
rotate(degrees(data_angle));
strokeWeight(1/data_scale);
mainFace.setProperties(params);
mainFace.draw(shifted_positions);
pop();
var scale_x = 80.0 / img.width;
var scale_y = 80.0 / img.height;
for(var j=0; j<4; j++) {
// var keys = Object.keys(trainData);
var curKey = curNeighbors[j];
var x2 = (3*width/4 - 200 + j*100);
push();
translate(x2, y2);
if (curKey in trainValues) {
var settings = trainValues[curKey]
translate(scale_x*data_mean[0], scale_y*data_mean[1]);
scale(scale_x*data_scale, scale_y*data_scale);
rotate(degrees(data_angle));
strokeWeight(1/data_scale);
littleFace.setProperties(settings);
littleFace.draw(shifted_positions);
}
else {
noFill();
stroke(100);
rect(10, 10, 70, 70);
}
pop();
}
}
textDisplay = "Neighbors: " + trainDataKeys[curTrainIndex];
}
else if (mode == 'TrainingQuiz' || mode == 'InterpolationQuiz') {
background(bg_color1);
var curKey = trainDataKeys[curTrainIndex];
var data = trainData[curKey];
var valid_mode = false;
if (mode == 'InterpolationQuiz') {
valid_mode = true;
curKey = validDataKeys[curValidIndex];
data = validData[curKey];
}
// Displays the image at its actual size at point (0,0)
var img = data.image
var x1 = (width/2-200/2);
var y1 = (height/3-300/2);
image(img, x1, y1, 200, 200);
if(valid_mode) {
fill(0, 0, 200);
}
else if (curKey in trainValues) {
fill(0, 200, 0);
}
else {
fill(200, 0, 0);
}
ellipse(x1+200/2, y1+200+15, 10, 10);
var y2 = (3*height/5-80/2);
var y3 = (4*height/5-80/2);
/*
for(var i=0; i<4; i++) {
// var keys = Object.keys(trainData);
var curKey = curNeighbors[i];
var nearData = trainData[curKey];
// Displays the image at its actual size at point (0,0)
var img = nearData.image
var x2 = (width/4 - 200 + i*100);
image(img, x2, y2, 80, 80);
}
*/
for(var i=0; i<1; i++) {
// get array of face marker positions [x, y] format
var positions = data.landmarks[i];
var shifted_positions = JSON.parse(JSON.stringify(positions))
var data_mean = [0.0, 0.0];
var data_scale = 1.0;
var data_angle = 0.0;
if ('transform' in positions) {
data_mean = positions.transform.center;
data_scale = positions.transform.scale;
data_angle = positions.transform.angle;
delete shifted_positions.transform
}
var scale_x = 400.0 / img.width;
var scale_y = 400.0 / img.height;
Object.keys(positions).forEach(function(key) {
if (key=='transform') {
return;
}
var curSection = positions[key];
var shiftedSection = shifted_positions[key];
for (var i=0; i<curSection.length; i++) {
var cur_x = curSection[i][0];
var cur_y = curSection[i][1];
// get ready for drawing the face
shiftedSection[i][0] = cur_x;
shiftedSection[i][1] = cur_y;
}
});
/*
var scale_x = 250.0 / img.width;
var scale_y = 250.0 / img.height;
var x2 = (3*width/4-250/2);
push();
translate(x2, y1);
translate(scale_x*data_mean[0], scale_y*data_mean[1]);
scale(scale_x*data_scale, scale_y*data_scale);
rotate(degrees(data_angle));
strokeWeight(1/data_scale);
mainFace.setProperties(params);
mainFace.draw(shifted_positions);
pop();
*/
var scale_x = 80.0 / img.width;
var scale_y = 80.0 / img.height;
var otherKeys = Object.keys(trainValues);
var index = otherKeys.indexOf(trainDataKeys[curTrainIndex]);
if(index >= 0) {
otherKeys.splice(index, 1);
}
var answerSlot = int(focusedRandom(0, 4));
var answerKeys = Array(4);
for(var j=0; j<4; j++) {
if(j == answerSlot) {
curKey = trainDataKeys[curTrainIndex];
}
else {
var guess = int(focusedRandom(0, otherKeys.length));
// if(otherKeys.length > j+2) {
// while(answerKeys.indexOf(guess) == -1) {
// guess = int(focusedRandom(0, otherKeys.length));
// }
// }
curKey = otherKeys[guess];
}
answerKeys[j] = curKey;
// print("Answer", j, " is ", curKey);
var x2 = (width/2 - 200 + j*100);
var settings = params;
if (valid_mode && j == answerSlot) {
var curEmbedding = data.embedding[0];
results = getAverageSettingsFrom(curEmbedding);
settings = results[0];
var validTrainKeys = results[1];
}
else if (curKey in trainValues) {
settings = trainValues[curKey];
}
push();
translate(x2, y2);
translate(scale_x*data_mean[0], scale_y*data_mean[1]);
scale(scale_x*data_scale, scale_y*data_scale);
rotate(degrees(data_angle));
strokeWeight(1/data_scale);
if (typeof settings !== 'undefined') {
littleFace.setProperties(settings);
}
littleFace.draw(shifted_positions);
pop();
if(quiz_done && guessed_answer == (j+1)) {
push();
translate(x2, y2);
noFill();
strokeWeight(4);
if(guessed_answer == (answerSlot+1)) {
stroke(0, 100, 0);
}
else {
stroke(100, 0, 0);
}
rect(-10, -10, 100, 100);
pop();
}
}
if(quiz_done) {
for(var j=0; j<4; j++) {
if (valid_mode && (answerSlot+1) == (j+1)) {
for(var k=0; k<4; k++) {
var curKey = validTrainKeys[k];
var nearData = trainData[curKey];
// Displays the image at its actual size at point (0,0)
var img = nearData.image
var x2 = (width/2 - 200 + j*100 + (k%2)*40);
var y4 = y3 + (int(k/2))*40;
image(img, x2, y4, 40, 40);
}
}
else {
var curKey = answerKeys[j];
var nearData = trainData[curKey];
// Displays the image at its actual size at point (0,0)
if (typeof nearData !== 'undefined') {
var img = nearData.image
var x2 = (width/2 - 200 + j*100);
image(img, x2, y3, 80, 80);
}
}
}
}
}
if(valid_mode) {
if(quiz_done) {
textDisplay = "InterpolationQuiz: hit a number to continue";
}
else {
textDisplay = "InterpolationQuiz: hit 1, 2, 3, or 4 to guess";
}
}
else {
if(quiz_done) {
textDisplay = "TrainingQuiz: hit a number to continue";
}
else {
textDisplay = "TrainingQuiz: hit 1, 2, 3, or 4 to guess";
}
}
}
fill(255);
textSize(32);
textAlign(CENTER);
text(textDisplay, width/2, height-12);
}
async function keyTyped() {
if(!haveStarted) {
return;
}
var mode = faceSelector.value();
if (key == 'q' && mode != 'Faces') {
faceSelector.value('Faces');
}
else if (key == 'w' && mode != 'Train') {
faceSelector.value('Train');
}
else if (key == 'e' && mode != 'NearestNeighbors') {
faceSelector.value('NearestNeighbors');
}
else if (key == 'r' && mode != 'TrainingQuiz') {
faceSelector.value('TrainingQuiz');
}
else if (key == 't' && mode != 'InterpolationQuiz') {
faceSelector.value('InterpolationQuiz');
}
else if (key == 'y' && mode != 'Video') {
faceSelector.value('Video');
}
if (key >= '0' && key <= '9' &&
(mode == 'TrainingQuiz' || mode == 'InterpolationQuiz') && quiz_done) {
quiz_done = false;
if(mode == 'TrainingQuiz') {
curTrainIndex = (curTrainIndex + 1) % trainDataKeys.length;
}
else {
curValidIndex = (curValidIndex + 1) % validDataKeys.length;
}
changeRandomSeed();
}
else if ((mode == 'TrainingQuiz' || mode == 'InterpolationQuiz') && quiz_done == false) {
if(key >= '1' && key <= '4') {
guessed_answer = key - '0';
quiz_done = true;
}
}
if (key == 's') {
saveCurrentSettings();
}
else if (key == 'i') {
interpolateCurrent();
}
else if (key == 'l') {
loadCurrentSettings();
}
if (key == '!') {
saveBlocksImages();
}
else if (key == '@') {
saveBlocksImages(true);
}
/* THIS CODE IS USED TO UPDATE TRAINING_IMAGES.JS and TESTING_IMAGES.JS
if (key == '/') {
print("loading new train[0]");
for(let i=0; i<trainDataKeys.length; i++) {
var curKey = trainDataKeys[i];
var data = trainData[curKey];
let fullFaceDescriptions = await faceapi.detectAllFaces(data.image.canvas).withFaceLandmarks().withFaceDescriptors();
data.landmarks = get_landmarks(fullFaceDescriptions);
data.embedding = get_latents(fullFaceDescriptions);
}
print("loaded new trains");
return;
}
if (key == '?') {
print("running diagnostic on training");
let obj = {};
for(let i=0; i<trainDataKeys.length; i++) {
var curKey = trainDataKeys[i];
var data = trainData[curKey];
let fullFaceDescriptions = await faceapi.detectAllFaces(data.image.canvas).withFaceLandmarks().withFaceDescriptors();
let found_embed = get_latents(fullFaceDescriptions);
let found_landmarks = get_landmarks(fullFaceDescriptions);
// print(data.embedding[0], found_embed[0]);
// print("dist old to new -> ", num_dist(found_embed[0], data.embedding[0]));
let lst = [];
for(let i=0; i<128; i++) {
lst[i] = found_embed[0][i];
}
obj[curKey] = {"url": data.url, "embedding": [lst], "landmarks": found_landmarks};
}
// print(obj);
// obj = found_embed;
var text = select('#output');
var json = JSON.stringify(obj, null);
// alert(json);
text.html(json)
return;
}
if (key == ':') {
print("running diagnostic on validation");
let obj = {};
for(let i=0; i<validDataKeys.length; i++) {
var curKey = validDataKeys[i];
var data = validData[curKey];
let fullFaceDescriptions = await faceapi.detectAllFaces(data.image.canvas).withFaceLandmarks().withFaceDescriptors();
let found_embed = get_latents(fullFaceDescriptions);
let found_landmarks = get_landmarks(fullFaceDescriptions);
// print(data.embedding[0], found_embed[0]);
// print("dist old to new -> ", num_dist(found_embed[0], data.embedding[0]));
let lst = [];
for(let i=0; i<128; i++) {
lst[i] = found_embed[0][i];
}
obj[curKey] = {"url": data.url, "embedding": [lst], "landmarks": found_landmarks};
}
// print(obj);
// obj = found_embed;
var text = select('#output');
var json = JSON.stringify(obj, null);
// alert(json);
text.html(json)
return;
}
*/
}
function interpolateCurrent() {
var curNeighborSettings = [];
for(var i=0; i<4; i++) {
neighborKey = curNeighbors[i]
if(neighborKey in trainValues) {
curNeighborSettings.push(trainValues[neighborKey]);
}
}
for(var i=0; i<NUM_SLIDERS; i++) {
sliders[i].value(50);
}
if(curNeighborSettings.length > 0) {
settings = curNeighborSettings[0];
for(var i=0; i<settings.length; i++) {
var sum = 0;
for(j=0; j<curNeighborSettings.length; j++) {
sum += curNeighborSettings[j][i];
}
var avg = int(sum / curNeighborSettings.length)
sliders[i].value(avg);
}
}
}
function loadCurrentSettings() {
var curKey = trainDataKeys[curTrainIndex];
var mainSettings = mainFace.getProperties();
for(var i=0; i<NUM_SLIDERS; i++) {
if(i < mainSettings.length) {
sliders[i].value(mainSettings[i])
}
else {
sliders[i].value(50);
}
}
if (curKey in trainValues) {
var settings = trainValues[curKey]
for(var i=0; i<settings.length; i++) {
sliders[i].value(settings[i]);
}
}
}
function updateSlidersForTraining() {
var mode = faceSelector.value();
var curKey = trainDataKeys[curTrainIndex];
// first find the closest neighbors
var nearest = allEmbeddingsTree.nearest(trainData[curKey].embedding[0], 5);
curNeighbors = [];
curNeighborSettings = [];
for(var i=0; i<5; i++) {
if(nearest[i][0][128] != curKey) {
var neighborKey = nearest[i][0][128];
curNeighbors.push(neighborKey);
if(neighborKey in trainValues) {
curNeighborSettings.push(trainValues[neighborKey]);
}
}
}
loadCurrentSettings();
// if(mode == 'NearestNeighbors') {
// interpolateCurrent();
// }
// else {
// loadCurrentSettings();
// }
}
function getAverageSettingsFrom(e) {
// first find the closest neighbors
var nearest = allEmbeddingsTree.nearest(e, 4);
curNeighbors = [];
curNeighborSettings = [];
for(var i=0; i<4; i++) {
var neighborKey = nearest[i][0][128];
curNeighbors.push(neighborKey);
if(neighborKey in trainValues) {
curNeighborSettings.push(trainValues[neighborKey]);
}
}
for(var i=0; i<4; i++) {
neighborKey = curNeighbors[i]
if(neighborKey in trainValues) {
curNeighborSettings.push(trainValues[neighborKey]);
}
}
var trainValueKeys = Object.keys(trainValues);
var props = trainValues[trainValueKeys[0]];
if(curNeighborSettings.length > 0) {
settings = curNeighborSettings[0];
for(var i=0; i<settings.length; i++) {
var sum = 0;
for(j=0; j<curNeighborSettings.length; j++) {
sum += curNeighborSettings[j][i];
}
var avg = int(sum / curNeighborSettings.length)
props[i] = avg;
}
}
return [props, curNeighbors];
}
function keyPressed() {
if(!haveStarted) {
return;
}
var mode = faceSelector.value();
if (mode == 'Faces') {
if (keyCode == LEFT_ARROW || keyCode == UP_ARROW) {
curFaceIndex = (curFaceIndex + selfieData.length - 1) % selfieData.length;
} else if (keyCode === RIGHT_ARROW || keyCode == DOWN_ARROW) {
curFaceIndex = (curFaceIndex + 1) % selfieData.length;
}
}
else if (mode == 'Train' || mode == 'NearestNeighbors') {
if (keyCode == LEFT_ARROW || keyCode == UP_ARROW) {
curTrainIndex = (curTrainIndex + trainDataKeys.length - 1) % trainDataKeys.length;
updateSlidersForTraining();
} else if (keyCode == RIGHT_ARROW || keyCode == DOWN_ARROW) {
curTrainIndex = (curTrainIndex + 1) % trainDataKeys.length;
updateSlidersForTraining();
}
}
}
Raw
z_clmtrackr.js
This file has been truncated, but you can view the full file.
(function (global, factory) {
typeof exports === 'object' && typeof module !== 'undefined' ? module.exports = factory() :
typeof define === 'function' && define.amd ? define(factory) :
(global.clm = factory());
}(this, (function () { 'use strict';
var commonjsGlobal = typeof window !== 'undefined' ? window : typeof global !== 'undefined' ? global : typeof self !== 'undefined' ? self : {};
function createCommonjsModule(fn, module) {
return module = { exports: {} }, fn(module, module.exports), module.exports;
}
var numeric1_2_6 = createCommonjsModule(function (module, exports) {
"use strict";
var numeric = exports;
if(typeof commonjsGlobal !== "undefined") { commonjsGlobal.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; }
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');
};
// seedrandom.js version 2.0.
// Author: David Bau 4/2/2011
//
// Defines a method Math.seedrandom() that, when called, substitutes
// an explicitly seeded RC4-based algorithm for Math.random(). Also
// supports automatic seeding from local or network sources of entropy.
//
// Usage:
//
// <script src=http://davidbau.com/encode/seedrandom-min.js></script>
//
// Math.seedrandom('yipee'); Sets Math.random to a function that is
// initialized using the given explicit seed.
//
// Math.seedrandom(); Sets Math.random to a function that is
// seeded using the current time, dom state,
// and other accumulated local entropy.
// The generated seed string is returned.
//
// Math.seedrandom('yowza', true);
// Seeds using the given explicit seed mixed
// together with accumulated entropy.
//
// <script src="http://bit.ly/srandom-512"></script>
// Seeds using physical random bits downloaded
// from random.org.
//
// <script src="https://jsonlib.appspot.com/urandom?callback=Math.seedrandom">
// </script> Seeds using urandom bits from call.jsonlib.com,
// which is faster than random.org.
//
// Examples:
//
// Math.seedrandom("hello"); // Use "hello" as the seed.
// document.write(Math.random()); // Always 0.5463663768140734
// document.write(Math.random()); // Always 0.43973793770592234
// var rng1 = Math.random; // Remember the current prng.
//
// var autoseed = Math.seedrandom(); // New prng with an automatic seed.
// document.write(Math.random()); // Pretty much unpredictable.
//
// Math.random = rng1; // Continue "hello" prng sequence.
// document.write(Math.random()); // Always 0.554769432473455
//
// Math.seedrandom(autoseed); // Restart at the previous seed.
// document.write(Math.random()); // Repeat the 'unpredictable' value.
//
// Notes:
//
// Each time seedrandom('arg') is called, entropy from the passed seed
// is accumulated in a pool to help generate future seeds for the
// zero-argument form of Math.seedrandom, so entropy can be injected over
// time by calling seedrandom with explicit data repeatedly.
//
// On speed - This javascript implementation of Math.random() is about
// 3-10x slower than the built-in Math.random() because it is not native
// code, but this is typically fast enough anyway. Seeding is more expensive,
// especially if you use auto-seeding. Some details (timings on Chrome 4):
//
// Our Math.random() - avg less than 0.002 milliseconds per call
// seedrandom('explicit') - avg less than 0.5 milliseconds per call
// seedrandom('explicit', true) - avg less than 2 milliseconds per call
// seedrandom() - avg about 38 milliseconds per call
//
// LICENSE (BSD):
//
// Copyright 2010 David Bau, all rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of this module nor the names of its contributors may
// be used to endorse or promote products derived from this software
// without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
/**
* All code is in an anonymous closure to keep the global namespace clean.
*
* @param {number=} overflow
* @param {number=} startdenom
*/
// Patched by Seb so that seedrandom.js does not pollute the Math object.
// My tests suggest that doing Math.trouble = 1 makes Math lookups about 5%
// slower.
numeric.seedrandom = { pow:Math.pow, random:Math.random };
(function (pool, math, width, chunks, significance, overflow, startdenom) {
//
// seedrandom()
// This is the seedrandom function described above.
//
math['seedrandom'] = function seedrandom(seed, use_entropy) {
var key = [];
var arc4;
// Flatten the seed string or build one from local entropy if needed.
seed = mixkey(flatten(
use_entropy ? [seed, pool] :
arguments.length ? seed :
[new Date().getTime(), pool, window], 3), key);
// Use the seed to initialize an ARC4 generator.
arc4 = new ARC4(key);
// Mix the randomness into accumulated entropy.
mixkey(arc4.S, pool);
// Override Math.random
// This function returns a random double in [0, 1) that contains
// randomness in every bit of the mantissa of the IEEE 754 value.
math['random'] = function random() { // Closure to return a random double:
var n = arc4.g(chunks); // Start with a numerator n < 2 ^ 48
var d = startdenom; // and denominator d = 2 ^ 48.
var x = 0; // and no 'extra last byte'.
while (n < significance) { // Fill up all significant digits by
n = (n + x) * width; // shifting numerator and
d *= width; // denominator and generating a
x = arc4.g(1); // new least-significant-byte.
}
while (n >= overflow) { // To avoid rounding up, before adding
n /= 2; // last byte, shift everything
d /= 2; // right using integer math until
x >>>= 1; // we have exactly the desired bits.
}
return (n + x) / d; // Form the number within [0, 1).
};
// Return the seed that was used
return seed;
};
//
// ARC4
//
// An ARC4 implementation. The constructor takes a key in the form of
// an array of at most (width) integers that should be 0 <= x < (width).
//
// The g(count) method returns a pseudorandom integer that concatenates
// the next (count) outputs from ARC4. Its return value is a number x
// that is in the range 0 <= x < (width ^ count).
//
/** @constructor */
function ARC4(key) {
var t, u, me = this, keylen = key.length;
var i = 0, j = me.i = me.j = me.m = 0;
me.S = [];
me.c = [];
// The empty key [] is treated as [0].
if (!keylen) { key = [keylen++]; }
// Set up S using the standard key scheduling algorithm.
while (i < width) { me.S[i] = i++; }
for (i = 0; i < width; i++) {
t = me.S[i];
j = lowbits(j + t + key[i % keylen]);
u = me.S[j];
me.S[i] = u;
me.S[j] = t;
}
// The "g" method returns the next (count) outputs as one number.
me.g = function getnext(count) {
var s = me.S;
var i = lowbits(me.i + 1); var t = s[i];
var j = lowbits(me.j + t); var u = s[j];
s[i] = u;
s[j] = t;
var r = s[lowbits(t + u)];
while (--count) {
i = lowbits(i + 1); t = s[i];
j = lowbits(j + t); u = s[j];
s[i] = u;
s[j] = t;
r = r * width + s[lowbits(t + u)];
}
me.i = i;
me.j = j;
return r;
};
// For robust unpredictability discard an initial batch of values.
// See http://www.rsa.com/rsalabs/node.asp?id=2009
me.g(width);
}
//
// flatten()
// Converts an object tree to nested arrays of strings.
//
/** @param {Object=} result
* @param {string=} prop
* @param {string=} typ */
function flatten(obj, depth, result, prop, typ) {
result = [];
typ = typeof(obj);
if (depth && typ == 'object') {
for (prop in obj) {
if (prop.indexOf('S') < 5) { // Avoid FF3 bug (local/sessionStorage)
try { result.push(flatten(obj[prop], depth - 1)); } catch (e) {}
}
}
}
return (result.length ? result : obj + (typ != 'string' ? '\0' : ''));
}
//
// mixkey()
// Mixes a string seed into a key that is an array of integers, and
// returns a shortened string seed that is equivalent to the result key.
//
/** @param {number=} smear
* @param {number=} j */
function mixkey(seed, key, smear, j) {
seed += ''; // Ensure the seed is a string
smear = 0;
for (j = 0; j < seed.length; j++) {
key[lowbits(j)] =
lowbits((smear ^= key[lowbits(j)] * 19) + seed.charCodeAt(j));
}
seed = '';
for (j in key) { seed += String.fromCharCode(key[j]); }
return seed;
}
//
// lowbits()
// A quick "n mod width" for width a power of 2.
//
function lowbits(n) { return n & (width - 1); }
//
// The following constants are related to IEEE 754 limits.
//
startdenom = math.pow(width, chunks);
significance = math.pow(2, significance);
overflow = significance * 2;
//
// When seedrandom.js is loaded, we immediately mix a few bits
// from the built-in RNG into the entropy pool. Because we do
// not want to intefere with determinstic PRNG state later,
// seedrandom will not call math.random on its own again after
// initialization.
//
mixkey(math.random(), pool);
// End anonymous scope, and pass initial values.
}(
[], // pool: entropy pool starts empty
numeric.seedrandom, // math: package containing random, pow, and seedrandom
256, // width: each RC4 output is 0 <= x < 256
6, // chunks: at least six RC4 outputs for each double
52 // significance: there are 52 significant digits in a double
));
/* This file is a slightly modified version of quadprog.js from Alberto Santini.
* It has been slightly modified by Sébastien Loisel to make sure that it handles
* 0-based Arrays instead of 1-based Arrays.
* License is in resources/LICENSE.quadprog */
(function(exports) {
function base0to1(A) {
if(typeof A !== "object") { return A; }
var ret = [], i,n=A.length;
for(i=0;i<n;i++) ret[i+1] = base0to1(A[i]);
return ret;
}
function base1to0(A) {
if(typeof A !== "object") { return A; }
var ret = [], i,n=A.length;
for(i=1;i<n;i++) ret[i-1] = base1to0(A[i]);
return ret;
}
function dpori(a, lda, n) {
var i, j, k, kp1, t;
for (k = 1; k <= n; k = k + 1) {
a[k][k] = 1 / a[k][k];
t = -a[k][k];
//~ dscal(k - 1, t, a[1][k], 1);
for (i = 1; i < k; i = i + 1) {
a[i][k] = t * a[i][k];
}
kp1 = k + 1;
if (n < kp1) {
break;
}
for (j = kp1; j <= n; j = j + 1) {
t = a[k][j];
a[k][j] = 0;
//~ daxpy(k, t, a[1][k], 1, a[1][j], 1);
for (i = 1; i <= k; i = i + 1) {
a[i][j] = a[i][j] + (t * a[i][k]);
}
}
}
}
function dposl(a, lda, n, b) {
var i, k, kb, t;
for (k = 1; k <= n; k = k + 1) {
//~ t = ddot(k - 1, a[1][k], 1, b[1], 1);
t = 0;
for (i = 1; i < k; i = i + 1) {
t = t + (a[i][k] * b[i]);
}
b[k] = (b[k] - t) / a[k][k];
}
for (kb = 1; kb <= n; kb = kb + 1) {
k = n + 1 - kb;
b[k] = b[k] / a[k][k];
t = -b[k];
//~ daxpy(k - 1, t, a[1][k], 1, b[1], 1);
for (i = 1; i < k; i = i + 1) {
b[i] = b[i] + (t * a[i][k]);
}
}
}
function dpofa(a, lda, n, info) {
var i, j, jm1, k, t, s;
for (j = 1; j <= n; j = j + 1) {
info[1] = j;
s = 0;
jm1 = j - 1;
if (jm1 < 1) {
s = a[j][j] - s;
if (s <= 0) {
break;
}
a[j][j] = Math.sqrt(s);
} else {
for (k = 1; k <= jm1; k = k + 1) {
//~ t = a[k][j] - ddot(k - 1, a[1][k], 1, a[1][j], 1);
t = a[k][j];
for (i = 1; i < k; i = i + 1) {
t = t - (a[i][j] * a[i][k]);
}
t = t / a[k][k];
a[k][j] = t;
s = s + t * t;
}
s = a[j][j] - s;
if (s <= 0) {
break;
}
a[j][j] = Math.sqrt(s);
}
info[1] = 0;
}
}
function qpgen2(dmat, dvec, fddmat, n, sol, crval, amat,
bvec, fdamat, q, meq, iact, nact, iter, work, ierr) {
var i, j, l, l1, info, it1, iwzv, iwrv, iwrm, iwsv, iwuv, nvl, r, iwnbv,
temp, sum, t1, tt, gc, gs, nu,
t1inf, t2min,
vsmall, tmpa, tmpb,
go;
r = Math.min(n, q);
l = 2 * n + (r * (r + 5)) / 2 + 2 * q + 1;
vsmall = 1.0e-60;
do {
vsmall = vsmall + vsmall;
tmpa = 1 + 0.1 * vsmall;
tmpb = 1 + 0.2 * vsmall;
} while (tmpa <= 1 || tmpb <= 1);
for (i = 1; i <= n; i = i + 1) {
work[i] = dvec[i];
}
for (i = n + 1; i <= l; i = i + 1) {
work[i] = 0;
}
for (i = 1; i <= q; i = i + 1) {
iact[i] = 0;
}
info = [];
if (ierr[1] === 0) {
dpofa(dmat, fddmat, n, info);
if (info[1] !== 0) {
ierr[1] = 2;
return;
}
dposl(dmat, fddmat, n, dvec);
dpori(dmat, fddmat, n);
} else {
for (j = 1; j <= n; j = j + 1) {
sol[j] = 0;
for (i = 1; i <= j; i = i + 1) {
sol[j] = sol[j] + dmat[i][j] * dvec[i];
}
}
for (j = 1; j <= n; j = j + 1) {
dvec[j] = 0;
for (i = j; i <= n; i = i + 1) {
dvec[j] = dvec[j] + dmat[j][i] * sol[i];
}
}
}
crval[1] = 0;
for (j = 1; j <= n; j = j + 1) {
sol[j] = dvec[j];
crval[1] = crval[1] + work[j] * sol[j];
work[j] = 0;
for (i = j + 1; i <= n; i = i + 1) {
dmat[i][j] = 0;
}
}
crval[1] = -crval[1] / 2;
ierr[1] = 0;
iwzv = n;
iwrv = iwzv + n;
iwuv = iwrv + r;
iwrm = iwuv + r + 1;
iwsv = iwrm + (r * (r + 1)) / 2;
iwnbv = iwsv + q;
for (i = 1; i <= q; i = i + 1) {
sum = 0;
for (j = 1; j <= n; j = j + 1) {
sum = sum + amat[j][i] * amat[j][i];
}
work[iwnbv + i] = Math.sqrt(sum);
}
nact = 0;
iter[1] = 0;
iter[2] = 0;
function fn_goto_50() {
iter[1] = iter[1] + 1;
l = iwsv;
for (i = 1; i <= q; i = i + 1) {
l = l + 1;
sum = -bvec[i];
for (j = 1; j <= n; j = j + 1) {
sum = sum + amat[j][i] * sol[j];
}
if (Math.abs(sum) < vsmall) {
sum = 0;
}
if (i > meq) {
work[l] = sum;
} else {
work[l] = -Math.abs(sum);
if (sum > 0) {
for (j = 1; j <= n; j = j + 1) {
amat[j][i] = -amat[j][i];
}
bvec[i] = -bvec[i];
}
}
}
for (i = 1; i <= nact; i = i + 1) {
work[iwsv + iact[i]] = 0;
}
nvl = 0;
temp = 0;
for (i = 1; i <= q; i = i + 1) {
if (work[iwsv + i] < temp * work[iwnbv + i]) {
nvl = i;
temp = work[iwsv + i] / work[iwnbv + i];
}
}
if (nvl === 0) {
return 999;
}
return 0;
}
function fn_goto_55() {
for (i = 1; i <= n; i = i + 1) {
sum = 0;
for (j = 1; j <= n; j = j + 1) {
sum = sum + dmat[j][i] * amat[j][nvl];
}
work[i] = sum;
}
l1 = iwzv;
for (i = 1; i <= n; i = i + 1) {
work[l1 + i] = 0;
}
for (j = nact + 1; j <= n; j = j + 1) {
for (i = 1; i <= n; i = i + 1) {
work[l1 + i] = work[l1 + i] + dmat[i][j] * work[j];
}
}
t1inf = true;
for (i = nact; i >= 1; i = i - 1) {
sum = work[i];
l = iwrm + (i * (i + 3)) / 2;
l1 = l - i;
for (j = i + 1; j <= nact; j = j + 1) {
sum = sum - work[l] * work[iwrv + j];
l = l + j;
}
sum = sum / work[l1];
work[iwrv + i] = sum;
if (iact[i] < meq) {
// continue;
break;
}
if (sum < 0) {
// continue;
break;
}
t1inf = false;
it1 = i;
}
if (!t1inf) {
t1 = work[iwuv + it1] / work[iwrv + it1];
for (i = 1; i <= nact; i = i + 1) {
if (iact[i] < meq) {
// continue;
break;
}
if (work[iwrv + i] < 0) {
// continue;
break;
}
temp = work[iwuv + i] / work[iwrv + i];
if (temp < t1) {
t1 = temp;
it1 = i;
}
}
}
sum = 0;
for (i = iwzv + 1; i <= iwzv + n; i = i + 1) {
sum = sum + work[i] * work[i];
}
if (Math.abs(sum) <= vsmall) {
if (t1inf) {
ierr[1] = 1;
// GOTO 999
return 999;
} else {
for (i = 1; i <= nact; i = i + 1) {
work[iwuv + i] = work[iwuv + i] - t1 * work[iwrv + i];
}
work[iwuv + nact + 1] = work[iwuv + nact + 1] + t1;
// GOTO 700
return 700;
}
} else {
sum = 0;
for (i = 1; i <= n; i = i + 1) {
sum = sum + work[iwzv + i] * amat[i][nvl];
}
tt = -work[iwsv + nvl] / sum;
t2min = true;
if (!t1inf) {
if (t1 < tt) {
tt = t1;
t2min = false;
}
}
for (i = 1; i <= n; i = i + 1) {
sol[i] = sol[i] + tt * work[iwzv + i];
if (Math.abs(sol[i]) < vsmall) {
sol[i] = 0;
}
}
crval[1] = crval[1] + tt * sum * (tt / 2 + work[iwuv + nact + 1]);
for (i = 1; i <= nact; i = i + 1) {
work[iwuv + i] = work[iwuv + i] - tt * work[iwrv + i];
}
work[iwuv + nact + 1] = work[iwuv + nact + 1] + tt;
if (t2min) {
nact = nact + 1;
iact[nact] = nvl;
l = iwrm + ((nact - 1) * nact) / 2 + 1;
for (i = 1; i <= nact - 1; i = i + 1) {
work[l] = work[i];
l = l + 1;
}
if (nact === n) {
work[l] = work[n];
} else {
for (i = n; i >= nact + 1; i = i - 1) {
if (work[i] === 0) {
// continue;
break;
}
gc = Math.max(Math.abs(work[i - 1]), Math.abs(work[i]));
gs = Math.min(Math.abs(work[i - 1]), Math.abs(work[i]));
if (work[i - 1] >= 0) {
temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
} else {
temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
}
gc = work[i - 1] / temp;
gs = work[i] / temp;
if (gc === 1) {
// continue;
break;
}
if (gc === 0) {
work[i - 1] = gs * temp;
for (j = 1; j <= n; j = j + 1) {
temp = dmat[j][i - 1];
dmat[j][i - 1] = dmat[j][i];
dmat[j][i] = temp;
}
} else {
work[i - 1] = temp;
nu = gs / (1 + gc);
for (j = 1; j <= n; j = j + 1) {
temp = gc * dmat[j][i - 1] + gs * dmat[j][i];
dmat[j][i] = nu * (dmat[j][i - 1] + temp) - dmat[j][i];
dmat[j][i - 1] = temp;
}
}
}
work[l] = work[nact];
}
} else {
sum = -bvec[nvl];
for (j = 1; j <= n; j = j + 1) {
sum = sum + sol[j] * amat[j][nvl];
}
if (nvl > meq) {
work[iwsv + nvl] = sum;
} else {
work[iwsv + nvl] = -Math.abs(sum);
if (sum > 0) {
for (j = 1; j <= n; j = j + 1) {
amat[j][nvl] = -amat[j][nvl];
}
bvec[nvl] = -bvec[nvl];
}
}
// GOTO 700
return 700;
}
}
return 0;
}
function fn_goto_797() {
l = iwrm + (it1 * (it1 + 1)) / 2 + 1;
l1 = l + it1;
if (work[l1] === 0) {
// GOTO 798
return 798;
}
gc = Math.max(Math.abs(work[l1 - 1]), Math.abs(work[l1]));
gs = Math.min(Math.abs(work[l1 - 1]), Math.abs(work[l1]));
if (work[l1 - 1] >= 0) {
temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
} else {
temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
}
gc = work[l1 - 1] / temp;
gs = work[l1] / temp;
if (gc === 1) {
// GOTO 798
return 798;
}
if (gc === 0) {
for (i = it1 + 1; i <= nact; i = i + 1) {
temp = work[l1 - 1];
work[l1 - 1] = work[l1];
work[l1] = temp;
l1 = l1 + i;
}
for (i = 1; i <= n; i = i + 1) {
temp = dmat[i][it1];
dmat[i][it1] = dmat[i][it1 + 1];
dmat[i][it1 + 1] = temp;
}
} else {
nu = gs / (1 + gc);
for (i = it1 + 1; i <= nact; i = i + 1) {
temp = gc * work[l1 - 1] + gs * work[l1];
work[l1] = nu * (work[l1 - 1] + temp) - work[l1];
work[l1 - 1] = temp;
l1 = l1 + i;
}
for (i = 1; i <= n; i = i + 1) {
temp = gc * dmat[i][it1] + gs * dmat[i][it1 + 1];
dmat[i][it1 + 1] = nu * (dmat[i][it1] + temp) - dmat[i][it1 + 1];
dmat[i][it1] = temp;
}
}
return 0;
}
function fn_goto_798() {
l1 = l - it1;
for (i = 1; i <= it1; i = i + 1) {
work[l1] = work[l];
l = l + 1;
l1 = l1 + 1;
}
work[iwuv + it1] = work[iwuv + it1 + 1];
iact[it1] = iact[it1 + 1];
it1 = it1 + 1;
if (it1 < nact) {
// GOTO 797
return 797;
}
return 0;
}
function fn_goto_799() {
work[iwuv + nact] = work[iwuv + nact + 1];
work[iwuv + nact + 1] = 0;
iact[nact] = 0;
nact = nact - 1;
iter[2] = iter[2] + 1;
return 0;
}
go = 0;
while (true) {
go = fn_goto_50();
if (go === 999) {
return;
}
while (true) {
go = fn_goto_55();
if (go === 0) {
break;
}
if (go === 999) {
return;
}
if (go === 700) {
if (it1 === nact) {
fn_goto_799();
} else {
while (true) {
fn_goto_797();
go = fn_goto_798();
if (go !== 797) {
break;
}
}
fn_goto_799();
}
}
}
}
}
function solveQP(Dmat, dvec, Amat, bvec, meq, factorized) {
Dmat = base0to1(Dmat);
dvec = base0to1(dvec);
Amat = base0to1(Amat);
var i, n, q,
nact, r,
crval = [], iact = [], sol = [], work = [], iter = [],
message;
meq = meq || 0;
factorized = factorized ? base0to1(factorized) : [undefined, 0];
bvec = bvec ? base0to1(bvec) : [];
// In Fortran the array index starts from 1
n = Dmat.length - 1;
q = Amat[1].length - 1;
if (!bvec) {
for (i = 1; i <= q; i = i + 1) {
bvec[i] = 0;
}
}
for (i = 1; i <= q; i = i + 1) {
iact[i] = 0;
}
nact = 0;
r = Math.min(n, q);
for (i = 1; i <= n; i = i + 1) {
sol[i] = 0;
}
crval[1] = 0;
for (i = 1; i <= (2 * n + (r * (r + 5)) / 2 + 2 * q + 1); i = i + 1) {
work[i] = 0;
}
for (i = 1; i <= 2; i = i + 1) {
iter[i] = 0;
}
qpgen2(Dmat, dvec, n, n, sol, crval, Amat,
bvec, n, q, meq, iact, nact, iter, work, factorized);
message = "";
if (factorized[1] === 1) {
message = "constraints are inconsistent, no solution!";
}
if (factorized[1] === 2) {
message = "matrix D in quadratic function is not positive definite!";
}
return {
solution: base1to0(sol),
value: base1to0(crval),
unconstrained_solution: base1to0(dvec),
iterations: base1to0(iter),
iact: base1to0(iact),
message: message
};
}
exports.solveQP = solveQP;
}(numeric));
/*
Shanti Rao sent me this routine by private email. I had to modify it
slightly to work on Arrays instead of using a Matrix object.
It is apparently translated from http://stitchpanorama.sourceforge.net/Python/svd.py
*/
numeric.svd= function svd(A) {
var temp;
//Compute the thin SVD from G. H. Golub and C. Reinsch, Numer. Math. 14, 403-420 (1970)
var prec= numeric.epsilon; //Math.pow(2,-52) // assumes double prec
var tolerance= 1.e-64/prec;
var itmax= 50;
var c=0;
var i=0;
var j=0;
var k=0;
var l=0;
var u= numeric.clone(A);
var m= u.length;
var n= u[0].length;
if (m < n) throw "Need more rows than columns"
var e = new Array(n);
var q = new Array(n);
for (i=0; i<n; i++) e[i] = q[i] = 0.0;
var v = numeric.rep([n,n],0);
// v.zero();
function pythag(a,b)
{
a = Math.abs(a);
b = Math.abs(b);
if (a > b)
return a*Math.sqrt(1.0+(b*b/a/a))
else if (b == 0.0)
return a
return b*Math.sqrt(1.0+(a*a/b/b))
}
//Householder's reduction to bidiagonal form
var f= 0.0;
var g= 0.0;
var h= 0.0;
var x= 0.0;
var y= 0.0;
var z= 0.0;
var s= 0.0;
for (i=0; i < n; i++)
{
e[i]= g;
s= 0.0;
l= i+1;
for (j=i; j < m; j++)
s += (u[j][i]*u[j][i]);
if (s <= tolerance)
g= 0.0;
else
{
f= u[i][i];
g= Math.sqrt(s);
if (f >= 0.0) g= -g;
h= f*g-s;
u[i][i]=f-g;
for (j=l; j < n; j++)
{
s= 0.0;
for (k=i; k < m; k++)
s += u[k][i]*u[k][j];
f= s/h;
for (k=i; k < m; k++)
u[k][j]+=f*u[k][i];
}
}
q[i]= g;
s= 0.0;
for (j=l; j < n; j++)
s= s + u[i][j]*u[i][j];
if (s <= tolerance)
g= 0.0;
else
{
f= u[i][i+1];
g= Math.sqrt(s);
if (f >= 0.0) g= -g;
h= f*g - s;
u[i][i+1] = f-g;
for (j=l; j < n; j++) e[j]= u[i][j]/h;
for (j=l; j < m; j++)
{
s=0.0;
for (k=l; k < n; k++)
s += (u[j][k]*u[i][k]);
for (k=l; k < n; k++)
u[j][k]+=s*e[k];
}
}
y= Math.abs(q[i])+Math.abs(e[i]);
if (y>x)
x=y;
}
// accumulation of right hand gtransformations
for (i=n-1; i != -1; i+= -1)
{
if (g != 0.0)
{
h= g*u[i][i+1];
for (j=l; j < n; j++)
v[j][i]=u[i][j]/h;
for (j=l; j < n; j++)
{
s=0.0;
for (k=l; k < n; k++)
s += u[i][k]*v[k][j];
for (k=l; k < n; k++)
v[k][j]+=(s*v[k][i]);
}
}
for (j=l; j < n; j++)
{
v[i][j] = 0;
v[j][i] = 0;
}
v[i][i] = 1;
g= e[i];
l= i;
}
// accumulation of left hand transformations
for (i=n-1; i != -1; i+= -1)
{
l= i+1;
g= q[i];
for (j=l; j < n; j++)
u[i][j] = 0;
if (g != 0.0)
{
h= u[i][i]*g;
for (j=l; j < n; j++)
{
s=0.0;
for (k=l; k < m; k++) s += u[k][i]*u[k][j];
f= s/h;
for (k=i; k < m; k++) u[k][j]+=f*u[k][i];
}
for (j=i; j < m; j++) u[j][i] = u[j][i]/g;
}
else
for (j=i; j < m; j++) u[j][i] = 0;
u[i][i] += 1;
}
// diagonalization of the bidiagonal form
prec= prec*x;
for (k=n-1; k != -1; k+= -1)
{
for (var iteration=0; iteration < itmax; iteration++)
{ // test f splitting
var test_convergence = false;
for (l=k; l != -1; l+= -1)
{
if (Math.abs(e[l]) <= prec)
{ test_convergence= true;
break
}
if (Math.abs(q[l-1]) <= prec)
break
}
if (!test_convergence)
{ // cancellation of e[l] if l>0
c= 0.0;
s= 1.0;
var l1= l-1;
for (i =l; i<k+1; i++)
{
f= s*e[i];
e[i]= c*e[i];
if (Math.abs(f) <= prec)
break
g= q[i];
h= pythag(f,g);
q[i]= h;
c= g/h;
s= -f/h;
for (j=0; j < m; j++)
{
y= u[j][l1];
z= u[j][i];
u[j][l1] = y*c+(z*s);
u[j][i] = -y*s+(z*c);
}
}
}
// test f convergence
z= q[k];
if (l== k)
{ //convergence
if (z<0.0)
{ //q[k] is made non-negative
q[k]= -z;
for (j=0; j < n; j++)
v[j][k] = -v[j][k];
}
break //break out of iteration loop and move on to next k value
}
if (iteration >= itmax-1)
throw 'Error: no convergence.'
// shift from bottom 2x2 minor
x= q[l];
y= q[k-1];
g= e[k-1];
h= e[k];
f= ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g= pythag(f,1.0);
if (f < 0.0)
f= ((x-z)*(x+z)+h*(y/(f-g)-h))/x;
else
f= ((x-z)*(x+z)+h*(y/(f+g)-h))/x;
// next QR transformation
c= 1.0;
s= 1.0;
for (i=l+1; i< k+1; i++)
{
g= e[i];
y= q[i];
h= s*g;
g= c*g;
z= pythag(f,h);
e[i-1]= z;
c= f/z;
s= h/z;
f= x*c+g*s;
g= -x*s+g*c;
h= y*s;
y= y*c;
for (j=0; j < n; j++)
{
x= v[j][i-1];
z= v[j][i];
v[j][i-1] = x*c+z*s;
v[j][i] = -x*s+z*c;
}
z= pythag(f,h);
q[i-1]= z;
c= f/z;
s= h/z;
f= c*g+s*y;
x= -s*g+c*y;
for (j=0; j < m; j++)
{
y= u[j][i-1];
z= u[j][i];
u[j][i-1] = y*c+z*s;
u[j][i] = -y*s+z*c;
}
}
e[l]= 0.0;
e[k]= f;
q[k]= x;
}
}
//vt= transpose(v)
//return (u,q,vt)
for (i=0;i<q.length; i++)
if (q[i] < prec) q[i] = 0;
//sort eigenvalues
for (i=0; i< n; i++)
{
//writeln(q)
for (j=i-1; j >= 0; j--)
{
if (q[j] < q[i])
{
// writeln(i,'-',j)
c = q[j];
q[j] = q[i];
q[i] = c;
for(k=0;k<u.length;k++) { temp = u[k][i]; u[k][i] = u[k][j]; u[k][j] = temp; }
for(k=0;k<v.length;k++) { temp = v[k][i]; v[k][i] = v[k][j]; v[k][j] = temp; }
// u.swapCols(i,j)
// v.swapCols(i,j)
i = j;
}
}
}
return {U:u,S:q,V:v}
};
});
var performanceNow = createCommonjsModule(function (module) {
// Generated by CoffeeScript 1.12.2
(function() {
var getNanoSeconds, hrtime, loadTime, moduleLoadTime, nodeLoadTime, upTime;
if ((typeof performance !== "undefined" && performance !== null) && performance.now) {
module.exports = function() {
return performance.now();
};
} else if ((typeof process !== "undefined" && process !== null) && process.hrtime) {
module.exports = function() {
return (getNanoSeconds() - nodeLoadTime) / 1e6;
};
hrtime = process.hrtime;
getNanoSeconds = function() {
var hr;
hr = hrtime();
return hr[0] * 1e9 + hr[1];
};
moduleLoadTime = getNanoSeconds();
upTime = process.uptime() * 1e9;
nodeLoadTime = moduleLoadTime - upTime;
} else if (Date.now) {
module.exports = function() {
return Date.now() - loadTime;
};
loadTime = Date.now();
} else {
module.exports = function() {
return new Date().getTime() - loadTime;
};
loadTime = new Date().getTime();
}
}).call(commonjsGlobal);
});
var root = typeof window === 'undefined' ? commonjsGlobal : window;
var vendors = ['moz', 'webkit'];
var suffix = 'AnimationFrame';
var raf = root['request' + suffix];
var caf = root['cancel' + suffix] || root['cancelRequest' + suffix];
for(var i = 0; !raf && i < vendors.length; i++) {
raf = root[vendors[i] + 'Request' + suffix];
caf = root[vendors[i] + 'Cancel' + suffix]
|| root[vendors[i] + 'CancelRequest' + suffix];
}
// Some versions of FF have rAF but not cAF
if(!raf || !caf) {
var last = 0
, id = 0
, queue = []
, frameDuration = 1000 / 60;
raf = function(callback) {
if(queue.length === 0) {
var _now = performanceNow()
, next = Math.max(0, frameDuration - (_now - last));
last = next + _now;
setTimeout(function() {
var cp = queue.slice(0);
// Clear queue here to prevent
// callbacks from appending listeners
// to the current frame's queue
queue.length = 0;
for(var i = 0; i < cp.length; i++) {
if(!cp[i].cancelled) {
try{
cp[i].callback(last);
} catch(e) {
setTimeout(function() { throw e }, 0);
}
}
}
}, Math.round(next));
}
queue.push({
handle: ++id,
callback: callback,
cancelled: false
});
return id
};
caf = function(handle) {
for(var i = 0; i < queue.length; i++) {
if(queue[i].handle === handle) {
queue[i].cancelled = true;
}
}
};
}
var index = function(fn) {
// Wrap in a new function to prevent
// `cancel` potentially being assigned
// to the native rAF function
return raf.call(root, fn)
};
var cancel = function() {
caf.apply(root, arguments);
};
var polyfill = function() {
root.requestAnimationFrame = raf;
root.cancelAnimationFrame = caf;
};
index.cancel = cancel;
index.polyfill = polyfill;
var promise = createCommonjsModule(function (module) {
(function (root) {
// Store setTimeout reference so promise-polyfill will be unaffected by
// other code modifying setTimeout (like sinon.useFakeTimers())
var setTimeoutFunc = setTimeout;
function noop() {}
// Polyfill for Function.prototype.bind
function bind(fn, thisArg) {
return function () {
fn.apply(thisArg, arguments);
};
}
function Promise(fn) {
if (typeof this !== 'object') throw new TypeError('Promises must be constructed via new');
if (typeof fn !== 'function') throw new TypeError('not a function');
this._state = 0;
this._handled = false;
this._value = undefined;
this._deferreds = [];
doResolve(fn, this);
}
function handle(self, deferred) {
while (self._state === 3) {
self = self._value;
}
if (self._state === 0) {
self._deferreds.push(deferred);
return;
}
self._handled = true;
Promise._immediateFn(function () {
var cb = self._state === 1 ? deferred.onFulfilled : deferred.onRejected;
if (cb === null) {
(self._state === 1 ? resolve : reject)(deferred.promise, self._value);
return;
}
var ret;
try {
ret = cb(self._value);
} catch (e) {
reject(deferred.promise, e);
return;
}
resolve(deferred.promise, ret);
});
}
function resolve(self, newValue) {
try {
// Promise Resolution Procedure: https://github.com/promises-aplus/promises-spec#the-promise-resolution-procedure
if (newValue === self) throw new TypeError('A promise cannot be resolved with itself.');
if (newValue && (typeof newValue === 'object' || typeof newValue === 'function')) {
var then = newValue.then;
if (newValue instanceof Promise) {
self._state = 3;
self._value = newValue;
finale(self);
return;
} else if (typeof then === 'function') {
doResolve(bind(then, newValue), self);
return;
}
}
self._state = 1;
self._value = newValue;
finale(self);
} catch (e) {
reject(self, e);
}
}
function reject(self, newValue) {
self._state = 2;
self._value = newValue;
finale(self);
}
function finale(self) {
if (self._state === 2 && self._deferreds.length === 0) {
Promise._immediateFn(function() {
if (!self._handled) {
Promise._unhandledRejectionFn(self._value);
}
});
}
for (var i = 0, len = self._deferreds.length; i < len; i++) {
handle(self, self._deferreds[i]);
}
self._deferreds = null;
}
function Handler(onFulfilled, onRejected, promise) {
this.onFulfilled = typeof onFulfilled === 'function' ? onFulfilled : null;
this.onRejected = typeof onRejected === 'function' ? onRejected : null;
this.promise = promise;
}
/**
* Take a potentially misbehaving resolver function and make sure
* onFulfilled and onRejected are only called once.
*
* Makes no guarantees about asynchrony.
*/
function doResolve(fn, self) {
var done = false;
try {
fn(function (value) {
if (done) return;
done = true;
resolve(self, value);
}, function (reason) {
if (done) return;
done = true;
reject(self, reason);
});
} catch (ex) {
if (done) return;
done = true;
reject(self, ex);
}
}
Promise.prototype['catch'] = function (onRejected) {
return this.then(null, onRejected);
};
Promise.prototype.then = function (onFulfilled, onRejected) {
var prom = new (this.constructor)(noop);
handle(this, new Handler(onFulfilled, onRejected, prom));
return prom;
};
Promise.all = function (arr) {
var args = Array.prototype.slice.call(arr);
return new Promise(function (resolve, reject) {
if (args.length === 0) return resolve([]);
var remaining = args.length;
function res(i, val) {
try {
if (val && (typeof val === 'object' || typeof val === 'function')) {
var then = val.then;
if (typeof then === 'function') {
then.call(val, function (val) {
res(i, val);
}, reject);
return;
}
}
args[i] = val;
if (--remaining === 0) {
resolve(args);
}
} catch (ex) {
reject(ex);
}
}
for (var i = 0; i < args.length; i++) {
res(i, args[i]);
}
});
};
Promise.resolve = function (value) {
if (value && typeof value === 'object' && value.constructor === Promise) {
return value;
}
return new Promise(function (resolve) {
resolve(value);
});
};
Promise.reject = function (value) {
return new Promise(function (resolve, reject) {
reject(value);
});
};
Promise.race = function (values) {
return new Promise(function (resolve, reject) {
for (var i = 0, len = values.length; i < len; i++) {
values[i].then(resolve, reject);
}
});
};
// Use polyfill for setImmediate for performance gains
Promise._immediateFn = (typeof setImmediate === 'function' && function (fn) { setImmediate(fn); }) ||
function (fn) {
setTimeoutFunc(fn, 0);
};
Promise._unhandledRejectionFn = function _unhandledRejectionFn(err) {
if (typeof console !== 'undefined' && console) {
console.warn('Possible Unhandled Promise Rejection:', err); // eslint-disable-line no-console
}
};
/**
* Set the immediate function to execute callbacks
* @param fn {function} Function to execute
* @deprecated
*/
Promise._setImmediateFn = function _setImmediateFn(fn) {
Promise._immediateFn = fn;
};
/**
* Change the function to execute on unhandled rejection
* @param {function} fn Function to execute on unhandled rejection
* @deprecated
*/
Promise._setUnhandledRejectionFn = function _setUnhandledRejectionFn(fn) {
Promise._unhandledRejectionFn = fn;
};
if ('object' !== 'undefined' && module.exports) {
module.exports = Promise;
} else if (!root.Promise) {
root.Promise = Promise;
}
})(commonjsGlobal);
});
// IE polyfill from MDN
(function () {
if (typeof window.CustomEvent === 'function') return false;
function CustomEvent (event, params) {
params = params || {bubbles : false, cancelable : false, detail: undefined};
var evt = document.createEvent('CustomEvent');
evt.initCustomEvent(event, params.bubbles, params.cancelable, params.detail);
return evt;
}
CustomEvent.prototype = window.Event.prototype;
window.CustomEvent = CustomEvent;
})();
function emitEvent(eventName) {
var evt = new CustomEvent(eventName, {'bubbles': true, 'cancelable': true});
document.dispatchEvent(evt);
}
/**
* Fast Fourier Transform
* 1D-FFT/IFFT, 2D-FFT/IFFT (radix-2)
*
* @author ryo / github.com/wellflat
* Based on https://github.com/wellflat/jslib with some tiny optimizations
*/
function FFT() {
var _n = 0, // order
_bitrev = null, // bit reversal table
_cstb = null; // sin/cos table
var _tre, _tim;
this.init = function (n) {
if(n !== 0 && (n & (n - 1)) === 0) {
_n = n;
_setVariables();
_makeBitReversal();
_makeCosSinTable();
} else {
throw new Error("init: radix-2 required");
}
};
// 1D-FFT
this.fft1d = function (re, im) {
fft(re, im, 1);
};
// 1D-IFFT
this.ifft1d = function (re, im) {
var n = 1/_n;
fft(re, im, -1);
for(var i=0; i<_n; i++) {
re[i] *= n;
im[i] *= n;
}
};
// 2D-FFT
this.fft2d = function (re, im) {
var i = 0;
// x-axis
for(var y=0; y<_n; y++) {
i = y*_n;
for(var x1=0; x1<_n; x1++) {
_tre[x1] = re[x1 + i];
_tim[x1] = im[x1 + i];
}
this.fft1d(_tre, _tim);
for(var x2=0; x2<_n; x2++) {
re[x2 + i] = _tre[x2];
im[x2 + i] = _tim[x2];
}
}
// y-axis
for(var x=0; x<_n; x++) {
for(var y1=0; y1<_n; y1++) {
i = x + y1*_n;
_tre[y1] = re[i];
_tim[y1] = im[i];
}
this.fft1d(_tre, _tim);
for(var y2=0; y2<_n; y2++) {
i = x + y2*_n;
re[i] = _tre[y2];
im[i] = _tim[y2];
}
}
};
// 2D-IFFT
this.ifft2d = function (re, im) {
var i = 0;
// x-axis
for(var y=0; y<_n; y++) {
i = y*_n;
for(var x1=0; x1<_n; x1++) {
_tre[x1] = re[x1 + i];
_tim[x1] = im[x1 + i];
}
this.ifft1d(_tre, _tim);
for(var x2=0; x2<_n; x2++) {
re[x2 + i] = _tre[x2];
im[x2 + i] = _tim[x2];
}
}
// y-axis
for(var x=0; x<_n; x++) {
for(var y1=0; y1<_n; y1++) {
i = x + y1*_n;
_tre[y1] = re[i];
_tim[y1] = im[i];
}
this.ifft1d(_tre, _tim);
for(var y2=0; y2<_n; y2++) {
i = x + y2*_n;
re[i] = _tre[y2];
im[i] = _tim[y2];
}
}
};
// core operation of FFT
function fft(re, im, inv) {
var d, h, ik, m, tmp, wr, wi, xr, xi,
n4 = _n >> 2;
// bit reversal
for(var l=0; l<_n; l++) {
m = _bitrev[l];
if(l < m) {
tmp = re[l];
re[l] = re[m];
re[m] = tmp;
tmp = im[l];
im[l] = im[m];
im[m] = tmp;
}
}
// butterfly operation
for(var k=1; k<_n; k<<=1) {
h = 0;
d = _n/(k << 1);
for(var j=0; j<k; j++) {
wr = _cstb[h + n4];
wi = inv*_cstb[h];
for(var i=j; i<_n; i+=(k<<1)) {
ik = i + k;
xr = wr*re[ik] + wi*im[ik];
xi = wr*im[ik] - wi*re[ik];
re[ik] = re[i] - xr;
re[i] += xr;
im[ik] = im[i] - xi;
im[i] += xi;
}
h += d;
}
}
}
// set variables
function _setVariables() {
if(typeof Uint8Array !== 'undefined') {
_bitrev = new Uint8Array(_n);
} else {
_bitrev = new Array(_n);
}
if(typeof Float64Array !== 'undefined') {
_cstb = new Float64Array(_n*1.25);
_tre = new Float64Array(_n*_n);
_tim = new Float64Array(_n*_n);
} else {
_cstb = new Array(_n*1.25);
_tre = new Array(_n*_n);
_tim = new Array(_n*_n);
}
}
// make bit reversal table
function _makeBitReversal() {
var i = 0,
j = 0,
k = 0;
_bitrev[0] = 0;
while(++i < _n) {
k = _n >> 1;
while(k <= j) {
j -= k;
k >>= 1;
}
j += k;
_bitrev[i] = j;
}
}
// make trigonometric function table
function _makeCosSinTable() {
var n2 = _n >> 1,
n4 = _n >> 2,
n8 = _n >> 3,
n2p4 = n2 + n4,
t = Math.sin(Math.PI/_n),
dc = 2*t*t,
ds = Math.sqrt(dc*(2 - dc)),
c = _cstb[n4] = 1,
s = _cstb[0] = 0;
t = 2*dc;
for(var i=1; i<n8; i++) {
c -= dc;
dc += t*c;
s += ds;
ds -= t*s;
_cstb[i] = s;
_cstb[n4 - i] = c;
}
if(n8 !== 0) {
_cstb[n8] = Math.sqrt(0.5);
}
for(var j=0; j<n4; j++) {
_cstb[n2 - j] = _cstb[j];
}
for(var k=0; k<n2p4; k++) {
_cstb[k + n2] = -_cstb[k];
}
}
}
var left_eye_filter = {"real": [1.5419219943717721, 0.40010880110578706, -0.79043641265342957, -1.2685464969238938, 0.39878117336167285, -1.0673489992245377, -0.079880838229404019, -0.45374680224191505, -0.043474097938900787, -0.31125662385352687, 0.17092430376098702, -0.29613086164846153, 0.5616469648110296, -1.559786848789493, 0.64513037997492662, -1.2899747976234162, 1.1761667998175334, -1.2899747976233551, 0.64513037997490474, -1.5597868487894897, 0.56164696481102505, -0.29613086164845964, 0.17092430376099094, -0.31125662385352959, -0.043474097938900787, -0.45374680224191177, -0.079880838229404658, -1.0673489992245357, 0.39878117336167307, -1.2685464969238942, -0.79043641265343012, 0.40010880110578717, -1.3820969331049027, 0.069560471269205768, -1.9786339579213206, -1.9807415717551982, -0.78667274410450883, -1.2217002325587256, -0.19150029104902774, -0.35131617290773243, -0.17646388464205803, -0.16672095020503441, -0.092298612924566523, -0.028899376452253527, -0.1314555696102146, -0.32892265898101813, -0.40987148655061206, 0.11741827111366547, -0.67254330182605138, -0.46007833291519956, -0.67215259521101001, -0.44871907432473013, -0.034749316729184583, 0.0055639281302433969, -0.17675902360981591, -0.26196208085032191, -0.36301254306387037, -0.33546767337818123, -0.6458889740799838, -1.1981932989987978, 0.12372650763830917, -1.4996172161865935, -2.4084298023013888, -2.0505291279591722, -1.7249706159518585, -2.277646289702639, -3.1259631743419591, -2.9656385065342015, -2.8480835086962011, -1.4260964500310189, -0.61792590829173544, -0.2611655301498782, -0.38519889843539723, -0.17511899827006483, -0.32808050503227176, 0.0076800871037463036, -0.18710828510427668, 0.1976534820339281, -0.55444453100465052, 0.14583567590328381, -0.69844971117515287, -0.90188577233526623, -0.53500016384583371, -0.044420751861669799, 0.014727914354086128, -0.28084584584371913, -0.29890408748685848, -0.39431380149336548, -0.39569215798819307, -0.74351999988258299, -0.82502198370631752, -1.851491897104155, -0.74302378668934244, 0.21156442062863762, -3.3061472495599986, -1.7990472945779568, -2.2193764251732282, -2.3438802466919251, -3.3615971067123311, -3.5383249085863708, -2.2639673745086588, -2.0271757806780748, -0.75242583405872232, -0.30143411016839378, -0.3625272253546275, -0.25489431004647689, -0.18928491561467081, -0.1179891518538482, 0.027920290231533224, -0.035472107498143821, -0.29008721857562259, -0.3604588674139817, -0.39156143807433802, -0.82222257402876564, -0.44979914971695928, -0.098136330355476253, 0.065628582466229365, -0.33607304327303128, -0.32161201323497779, -0.41856090178723965, -0.64028425429629054, -0.7766428172010218, -1.3946448661671447, -2.2603422126144683, -0.38769722219534525, -0.95341593939478653, -1.412952994959813, -2.3602336858020432, -1.2756392437278019, -2.0983496132652038, -2.5682454610054268, -2.8791053946930378, -2.1809972632688095, -0.84281293847776861, -0.75998936793718697, -0.18584599820380068, -0.30105748355308259, -0.16098142942852958, -0.13792125740417191, -0.089790022871128708, -0.12321821342876504, -0.1128661923016878, -0.3924098378001975, -0.5780902167586397, -0.48685989567066695, -0.53565359443296234, -0.051036689850526382, -0.0068547033925117689, -0.18963405157839419, -0.22514761090777807, -0.35555823460888908, -0.46670603976585517, -0.56179541485257889, -0.7495095888115163, -1.4772075422260349, -1.5836466114968029, -2.3846549454186694, -1.4884613952536236, -1.8237453905245253, -1.6712324532934877, -1.5169157844507295, -1.6930052820597281, -2.1023566589276004, -2.2062031109308458, -1.7945281756942255, -0.26457398838912649, 0.22038139379151148, -0.43479836723775234, -0.19830827357221226, -0.18018565146479498, -0.097060879184795737, -0.10088329756370379, -0.063069709957272527, -0.17970932516041177, -0.1943040732581543, -0.37970560392277619, -0.47302301606251812, -0.30366967948052181, -0.064732391018915397, -0.08902516330269715, -0.082000200083027344, -0.22965854401457736, -0.32035624605031326, -0.31836783196552437, -0.40132058236311119, -0.65601747033470859, -0.59040483751417483, -1.8503084663080034, -1.8694842425148914, -1.9326778896298584, -1.6301578422923519, -1.4332006785118301, -1.305707665299106, -1.364200787821644, -1.5357935460809622, -1.6161992336951241, -0.74003518668370516, -0.29423824173210689, 0.025934598230976654, -0.043349004411304674, -0.25408021803022468, -0.066965686484977499, -0.075717498698635255, 0.007057189465364498, -0.042171356658338113, -0.036938315661768008, -0.34221561581756049, -0.20400167508805764, -0.37417116097079772, -0.25039909487805356, -0.070874531394524931, -0.0569972852039487, -0.067238206950403182, -0.17397285212300442, -0.20428337307808273, -0.23651154356493315, -0.33356498933276568, -0.07339749754226077, -0.70367959806681601, -0.82403680021595049, -1.6058616381755235, -1.6192427030685497, -1.5705638815427956, -1.4659201063980019, -0.95504179549951018, -0.97237526162739873, -1.0460191987834688, -0.91465668941265721, -0.60548232361398524, 0.01898438364933451, -0.19419044456729498, -0.039627851124307223, 0.0012357796666701798, -0.078110822445325079, 0.0048626364920250518, -0.040449089662379589, -0.0035054269587873454, -0.13387544724730729, -0.10031131456276647, -0.25968674675684189, -0.20555329767005767, -0.26509289948725284, -0.038788452621647145, -0.076999891872251258, -0.071661433038976499, -0.14182240789719938, -0.1654673053291095, -0.19859450279267193, -0.053382326365810369, -0.2156585383674445, -0.045097357284793499, -0.62449818579949512, -0.92624906744917224, -1.0411254782363617, -1.122035196738675, -1.0607692164246043, -0.57723811773534028, -0.63187735896388075, -0.54813311204421922, -0.55320252101738743, -0.30197299587482401, -0.047213249757838388, 0.082808930467383288, -0.067715134483222431, -0.01022881748368659, 0.042038311258956552, -0.063371767399980669, 0.029161890169972702, -0.091396316586836127, -0.0034600735070754811, -0.12424052925006424, -0.24432996418012101, -0.26521664175359499, -0.22745980283820413, -0.14361316535317664, -0.00075904203100577935, -0.020936168457862139, -0.14205665196423617, -0.19024248288823023, -0.079686122362245204, -0.15016133237735926, 0.049598910651295514, -0.11760486834511712, -0.1837522251545049, -0.38594205494114608, -0.53542516436999843, -0.57340991730807989, -0.52753621424018138, -0.23151163972118355, -0.22295096919949259, -0.33704349161770436, -0.26165852514054583, -0.13898866968588663, 0.034596483191139484, -0.012631210076789067, 0.047371310076345617, -0.038651839330751551, -0.0019970761454430018, 0.063048845258375494, -0.1124891762554399, 0.08556992539656616, -0.21043659051868199, -0.19223333969456, -0.39082994830035861, -0.19294368007162721, -0.41025595439938572, -0.17178084419175166, -0.010933041190555012, -0.089512936152074493, -0.21569610281495066, -0.09144756671688016, -0.19525258909505316, -0.029753598134641936, -0.021307245660079924, 0.029087127940551009, 0.037511290653097842, -0.20600990120705839, -0.26967580750352926, -0.21000923681194664, -0.28209018858285628, -0.11925518789339556, -0.24869348141289982, -0.21025892926356746, -0.15567029136726124, -0.040546729108395907, -0.0050266153100547101, 0.030710887069787196, -0.0061104340245858278, 0.0369376092260571, -0.054862661367900321, 0.013297880203253048, 0.19659447375886394, -0.2499491329142558, -0.062959699002865757, -0.53055029095956008, -0.38784811281629444, -0.53891285075962392, -0.41886712861154285, -0.099230097260325875, -0.16474199810952628, -0.28693665642627014, -0.0095667980850221105, -0.32619954993450928, -0.08627491478166284, -0.073253161755714766, 0.015634174038690329, 0.082440536547531792, 0.025411878261881942, -0.11318909242737961, -0.1270560226842935, -0.21657212936164139, -0.13993873549611191, -0.37510275237622831, -0.26472923111076219, -0.24460131567533192, -0.14127652303494026, -0.050428686591045178, 0.041347840374190772, -0.0061780445153000636, 0.0073990345210250153, -0.014062739037014381, 0.14348925152561878, -0.015321787554403667, 0.0017746672356015968, 0.25165135427361052, -0.626463828190993, -0.48167134330805639, -1.045863293770664, -0.69512591788493194, -0.44532127384388254, -0.28479724025368391, -0.39470955087317983, 0.20227228344720469, -0.53909912073488953, -0.12025629051789474, -0.1899243750597305, -0.048474806721595133, 0.060764771353227762, 0.090648151782516159, 0.091608208912697275, 0.0036582478916540977, -0.22492530005263131, -0.27295314658024766, -0.35559738025257359, -0.62902925014412947, -0.57166411974881004, -0.37258895173129181, -0.22157638610464933, 0.022494427132080854, 0.014769425415166171, 0.003526808789406817, -0.011346909674078769, 0.050921170848348289, 0.090308541799219627, 0.37260817254533324, -0.25909871392159911, -0.42379280974334355, -0.095380647808568128, -1.1906083748893519, -0.78599914414892469, -0.95277914352730275, -0.63659778359422337, -0.98026015008952749, 0.48173198285916102, -0.60092009018055192, -0.10265418316164113, -0.39913639006279306, -0.17310908908773887, -0.0194191171632387, 0.054047965289179878, 0.1388529643463832, 0.15661099050145999, -0.10898263774416243, -0.33291231456737602, -0.59569027865888713, -0.69353081584948972, -1.0999707493347484, -0.74392084753736687, -0.49074781214158159, -0.065190556733852961, 0.012289768389229717, 0.024577513704595676, 0.0040302804696096322, 0.036047756292976456, 0.058236765637246286, 0.13893846256790621, 0.036944676036934632, 0.41686279554239464, -0.85232286388185818, -1.2988315127624981, -0.47352779677305168, -0.81763632541546793, -0.77384457803621831, -1.4256240004519281, 0.52588993532360684, -0.89821724022902683, 0.1591911967653899, -0.55046596772346867, -0.30980016041271019, -0.16709614007114884, -0.046029700131955266, 0.044793268150423983, 0.1689242242845459, 0.14412365934528507, -0.0088250071313367359, -0.36778545124666312, -0.79393844517732104, -1.1610479066529615, -0.76523210008850662, -0.63009858032048405, -0.13947023057344932, -0.017173105577524262, 0.039030007688455846, 0.014491273083805401, 0.039792542943837252, 0.054072846696920814, 0.11729310469925348, 0.053609281522667675, 0.0081549498718087084, -0.30910813452845548, 0.25944224899607843, -1.3584842180322938, -1.5885570490138659, -0.65759582794618221, -1.139869490652734, 0.70928264080594694, -1.9674198903133462, 0.37712664425406606, -0.84336038390578949, -0.47788074719428036, -0.18342000086663721, -0.18811394573901796, -0.055050027645985648, 0.045043056834335606, 0.11486303559854361, 0.22023958716404868, 0.14735402009444676, -0.27894427087197998, -0.73080536953129638, -0.76794305693297227, -0.37355919765840223, 0.12353986794322802, 0.090505348376311842, 0.14069908672094206, 0.087373214380278855, 0.023353946735568523, 0.031400559920396587, 0.079550230446202241, 0.084927161382185437, 0.040777158255349423, -0.16274954314482293, -0.41184413435479567, -0.71871288822574875, 0.55302907456342854, -1.5309493464500674, -2.9026104205694736, 0.42043303599508353, -1.7138106264793671, 0.29513888249127102, -1.2517216433630918, -0.66769942176516839, -0.28576739334390183, -0.24127777006787937, -0.10778095858902549, -0.036092425009198861, 0.021519213385077923, 0.13414694961717147, 0.16917378957839613, 0.17307922682581758, 0.076246758829015673, -0.047904835134272621, -0.27544262702406924, 0.61826249566563185, 0.26987423123693399, 0.2085883517320696, 0.26073426210721973, 0.12070625812911842, 0.062945582093309679, 0.083649573916505432, 0.049688095345785867, 0.019564357607843069, -0.046035817476596949, -0.13409074070830324, -0.49027201814294552, -0.47756457321420159, -0.74403675135427549, -0.3080068432033089, -0.043712438842705037, -4.735594317158907, -0.043712438842706695, -0.30800684320330962, -0.74403675135427572, -0.47756457321420304, -0.49027201814294813, -0.13409074070830412, -0.046035817476598156, 0.019564357607843069, 0.049688095345786006, 0.083649573916506056, 0.062945582093310845, 0.12070625812911921, 0.26073426210722073, 0.20858835173207019, 0.26987423123693399, -0.37355919765836759, -0.27544262702403433, -0.047904835134273127, 0.076246758829012523, 0.17307922682581853, 0.16917378957839499, 0.13414694961716844, 0.02151921338507657, -0.036092425009199861, -0.1077809585890261, -0.24127777006787943, -0.2857673933439015, -0.66769942176516905, -1.2517216433630949, 0.29513888249127429, -1.7138106264793713, 0.42043303599507681, -2.902610420569474, -1.5309493464500692, 0.55302907456342232, -0.71871288822575019, -0.41184413435479833, -0.16274954314482265, 0.04077715825534866, 0.08492716138218645, 0.079550230446203143, 0.031400559920398419, 0.023353946735571576, 0.08737321438028138, 0.14069908672095732, 0.090505348376334033, 0.1235398679432393, -0.76523210008847808, -0.76794305693296139, -0.73080536953128505, -0.27894427087197604, 0.1473540200944477, 0.22023958716404682, 0.11486303559854165, 0.045043056834333829, -0.055050027645986453, -0.18811394573901843, -0.18342000086663854, -0.47788074719428042, -0.84336038390579149, 0.37712664425406617, -1.9674198903133469, 0.70928264080593695, -1.1398694906527307, -0.65759582794619398, -1.588557049013867, -1.3584842180322987, 0.25944224899607732, -0.30910813452845781, 0.0081549498718086911, 0.053609281522667279, 0.11729310469925426, 0.054072846696921202, 0.039792542943838709, 0.014491273083807311, 0.039030007688458185, -0.017173105577517028, -0.13947023057343994, -0.63009858032045107, -1.0999707493347308, -1.1610479066529467, -0.79393844517731305, -0.3677854512466584, -0.0088250071313340107, 0.14412365934528559, 0.16892422428454401, 0.044793268150420118, -0.046029700131956147, -0.16709614007115095, -0.30980016041271097, -0.55046596772347045, 0.15919119676539073, -0.8982172402290286, 0.52588993532360329, -1.4256240004519327, -0.77384457803621687, -0.8176363254154656, -0.47352779677305679, -1.2988315127625027, -0.85232286388185829, 0.41686279554239525, 0.036944676036935756, 0.13893846256790574, 0.058236765637246675, 0.036047756292977066, 0.0040302804696111128, 0.02457751370459911, 0.012289768389232913, -0.065190556733844662, -0.49074781214156804, -0.74392084753735632, -0.62902925014412903, -0.69353081584948562, -0.59569027865888302, -0.33291231456737491, -0.10898263774416028, 0.15661099050145985, 0.13885296434638142, 0.054047965289177706, -0.019419117163239467, -0.17310908908773912, -0.399136390062
Raw
z_face-api.js
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_face-api.js.map
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_face_landmark_68_model-shard1
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_face_landmark_68_model-weights_manifest.json
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_face_recognition_model-shard1
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_face_recognition_model-shard2
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_face_recognition_model-weights_manifest.json
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_face_system.js
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_focused_random.js
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_kdTree.js
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_model_pca_20_svm.js
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_purview_helper.js
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_ssd_mobilenetv1_model-shard1
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_ssd_mobilenetv1_model-shard2
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_ssd_mobilenetv1_model-weights_manifest.json
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_testing_images.js
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
z_training_images.js
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment