Skip to content

Instantly share code, notes, and snippets.

@dribnet

dribnet/.block

Forked from uwcc/.block
Last active October 18, 2020 21:50
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save dribnet/547506869d1051ffb76b380c141932a1 to your computer and use it in GitHub Desktop.
Save dribnet/547506869d1051ffb76b380c141932a1 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

Animating eyebrows and mouth based on relationship between face points.

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 = 8;
// other variables can be in here too
// these control the colors used
const bg_color = [225, 206, 187];
const fg_color = [151, 102, 52];
const stroke_color = [95, 52, 8];
function Face() {
// these are state variables for a face
// (your variables may be different)
this.bread = 3; // 1-5
this.jam = 5; // 0-10
/*
* Draw a 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) {
// let bread = int(map(s1, 0, 100, 1, 5));
// let mouth = map(s2, 0, 100, 0, 10);
// let blush = map(s3, 0, 100, 0, 10);
// let jam = map(s4, 0, 100, 0, 10);
// let jam_col = map(s5, 0, 100, 0, 10);
// let jam_size = int(map(s6, 0, 100, 0, 10));
// let eyebrows_rotate = map(s7, 0, 100, 0, 20);
// let glasses = int(map(s8, 0, 100, 0, 10));
// let pupil = int(map(s9, 0, 100, 1, 10));
let bread = this.bread;
let mouth = 5;
let blush = 5;
let jam = this.jam;
let jam_col = 5;
let jam_size = 5;
let eyebrows_rotate=10;
let glasses = 5;
let pupil = 5;
// Colour Sets
let bread1 = '#fffaf5';
let bread2 = '#f5f2ea';
let bread3 = '#dbb58c';
let bread4 = '#c48d52';
let bread5 = '#36210d';
let breadColour;
let left_eye = average_point(positions.left_eye);
let right_eye = average_point(positions.right_eye);
let left_eyebrow = average_point(positions.left_eyebrow);
let right_eyebrow = average_point(positions.right_eyebrow);
let left_d = dist(left_eye[0], left_eye[1], left_eyebrow[0], left_eyebrow[1]);
let right_d = dist(right_eye[0], right_eye[1], right_eyebrow[0], right_eyebrow[1]);
let left_eb_move = map(left_d, 0.4, 0.7, 0, 2, true);
let right_eb_move = map(right_d, 0.4, 0.7, 0, 2, true);
// print(left_d);
left_eye[0] *= 3;
left_eye[1] *= 3;
right_eye[0] *= 3;
right_eye[1] *= 3;
push();
scale(0.33);
// Bread (Brown)
fill(221,168,128);
strokeWeight(0.3);
stroke(0);
beginShape();
vertex(-6.45, -5);
quadraticVertex(-4, -8, -0.5, -6);
quadraticVertex(3.5, -9, 6.5, -5.6);
quadraticVertex(7.8, -3.8, 6.5, -1.5);
quadraticVertex(6.8, 3, 6.3, 5.5);
quadraticVertex(6, 6, 5.1, 6.5);
endShape();
// Bread Colour
if(bread == 1){
breadColour = bread1; //***ALMOST WHITE***//
}else if (bread == 2){
breadColour = bread2; //***IVORY***//
}else if (bread == 3){
breadColour = bread3; //***WELL DONE***//
}else if (bread == 4){
breadColour = bread4; //***OVERCOOKED***//
}else{
breadColour = bread5; //***BURNT***//
}
fill(breadColour);
strokeWeight(0.3);
stroke(0);
rect(-6, -4.5, 11, 11);
push();
angleMode(RADIANS);
arc(-3.9, -3, 6.8, 6, -4, -0.4);
arc(2.4, -3.2, 6.8, 6, -2.8, 0.65);
pop();
// Jam (Colour, Spread)
let amount = map(jam_col, 0, 10, 0, 1);
let colourPink = color('#ed3267');
let colourPurple = color('#edd532');
let jamColour = lerpColor(colourPink, colourPurple, amount);
let jamSpread = map(jam, 0, 100, 0, 100);
let jamS = map(jam_size, 0, 100, 1.3, 15);
strokeWeight(jamS);
stroke(jamColour);
if (jamSpread == 0){
}
if (jamSpread >= 2){
line(-3.5, -2, 2.5, -3);
}
if (jamSpread >=4){
line(2.5, -3, -3, 1);
}
if (jamSpread >=6){
line(-3, 1, 2, 0);
}
if (jamSpread >=8){
line(2, 0, -3.5, 4);
}
if (jamSpread >=10){
line(-3.5, 4, 2.5, 3);
}
// Mouth
strokeWeight(0.4);
stroke(0);
line(-1.55, -1, -0.05, -1);
let top_lip_point = positions.top_lip[9];
let bottom_lip_point = positions.bottom_lip[9];
// fill(255, 0, 0);
let d = dist(top_lip_point[0], top_lip_point[1], bottom_lip_point[0], bottom_lip_point[1])
// print(d);
// Mouth Open
if(d < 0.1) {
d = 0;
}
mouth = map(d, 0, 0.5, 0, 10);
let mouth_size = map(mouth, 0, 10, 0, 3.5);
strokeWeight(0.3);
fill('#fffcf7');
rect(-1.6, -1, 1.6, mouth_size, 0.08);
/*
push();
scale(5);
noStroke();
textSize(0.2);
for (let i=0; i<positions.bottom_lip.length; i++) {
let p = positions.bottom_lip[i]
fill(255, 0, 0);
ellipse(p[0], p[1], 0.2);
fill(0);
text(i, p[0], p[1]);
}
pop();
push();
scale(5);
noStroke();
textSize(0.2);
for (let i=0; i<positions.top_lip.length; i++) {
let p = positions.top_lip[i]
fill(0, 0, 255);
ellipse(p[0], p[1], 0.2);
fill(0);
text(i, p[0], p[1]);
}
pop();
*/
// Blush
noStroke();
if (blush > 3 && blush < 6){ //***PINK BLUSH***//
fill('#ff87b1');
ellipse(-2.6, -1.1, 1.5, 0.8);
ellipse(1.1, -1.1, 1.5, 0.8);
}
else if (blush > 5 && blush <= 8){ //***RED BLUSH***//
stroke('#bf304d');
strokeWeight(0.2);
line(-3.5, -0.9, -3.2, -1.2); //***LEFT***//
line(-2.9, -0.9, -2.6, -1.2);
line(-2.3, -0.9, -2, -1.2);
line(0.5, -0.9, 0.8, -1.2); //***RIGHT***//
line(1.1, -0.9, 1.4, -1.2);
line(1.7, -0.9, 2, -1.2);
}
else if (blush >7 && blush <= 10){ //***TEARS***//
noStroke();
fill('#49a3e6');
triangle(-3.3, -1.19, -2.9, -0.85, -2.75, -1.55);
ellipse(-3.2, -0.9, 0.6, 0.6);
triangle(1.81, -1.24, 1.36, -0.87, 1.2, -1.7);
ellipse(1.65, -1, 0.6, 0.6);
}
// Eyes
noStroke();
fill(0);
// ellipse(-2.2, -2, 1.1, 1.1);
// ellipse(0.5, -2, 1.1, 1.1);
ellipse(left_eye[0], left_eye[1], 1.1, 1.1);
ellipse(right_eye[0], right_eye[1], 1.1, 1.1);
// Pupil
if (pupil > 3 && pupil < 6) {
fill('black');
ellipse(left_eye[0]-0.1, left_eye[1]-0.1, 0.3, 0.3); //***LEFT TOP***//
ellipse(left_eye[0]+0.15, left_eye[1]+0.1, 0.2, 0.2); //***LEFT BOTTOM***//
ellipse(left_eye[0]-0.1, left_eye[1]-0.1, 0.3, 0.3); //***RIGHT TOP***//
ellipse(left_eye[0]+0.15, left_eye[1]+0.1, 0.2, 0.2); //***RIGHT BOTTOM***//
// ellipse(-2.3, -2.1, 0.3, 0.3); //***LEFT TOP***//
// ellipse(-2.05, -1.9, 0.2, 0.2); //***LEFT BOTTOM***//
// ellipse(0.4, -2.1, 0.3, 0.3); //***RIGHT TOP***//
// ellipse(0.65, -1.9, 0.2, 0.2); //***RIGHT BOTTOM***//
}
else if (pupil > 5 && pupil <= 10) {
// fill('white');
// ellipse(-2.3, -2.1, 0.3, 0.3); //***LEFT TOP***//
// ellipse(-2.05, -1.9, 0.2, 0.2); //***LEFT BOTTOM***//
// ellipse(0.4, -2.1, 0.3, 0.3); //***RIGHT TOP***//
// ellipse(0.65, -1.9, 0.2, 0.2); //***RIGHT BOTTOM***//
}
// Eyebrows
strokeWeight(0.4);
stroke(0);
push();
angleMode(DEGREES);
translate(left_eye[0]-0.1, left_eye[1]-1);
translate(0, -left_eb_move);
rotate(eyebrows_rotate);
line(-0.3, 0, 0.5, 0);
pop();
push();
angleMode(DEGREES);
translate(right_eye[0]-0.2, right_eye[1]-1);
translate(0, -right_eb_move);
rotate(-eyebrows_rotate);
line(-0.3, 0, 0.5, 0);
pop();
// Glasses
if (glasses > 5 && glasses < 8){
stroke(0);
strokeWeight(0.35);
noFill();
ellipse(-2.3, -2.4, 2.6, 2.6);
ellipse(0.6, -2.4, 2.6, 2.6);
}
// Sunglasses
else if (glasses > 7 && glasses <= 10){
stroke(0);
strokeWeight(0.3);
fill('white');
ellipse(-2.6, -2.4, 3, 2.6);
ellipse(1, -2.4, 3, 2.6);
noStroke();
fill('black');
ellipse(-2.6, -2.37, 2.4, 2);
ellipse(1, -2.37, 2.4, 2);
stroke(0);
strokeWeight(0.3);
line(-1.1, -2.5, -0.65, -2.5);
}
pop();
}
/* set internal properties based on list numbers 0-100 */
this.setProperties = function(settings) {
this.bread = int(map(settings[0], 0, 100, 1, 5));
this.jam = int(map(settings[1], 0, 100, 0, 10));
}
/* get internal properties as list of numbers 0-100 */
this.getProperties = function() {
let settings = new Array(2);
settings[0] = map(this.bread, 1, 5, 0, 100);
settings[1] = map(this.jam, 0, 10, 0, 100);
return settings;
}
}
// given an array of [x,y] points, return the average
function average_point(list) {
var sum_x = 0;
var sum_y = 0;
var num_points = 0;
for(var i=0; i<list.length; i++) {
sum_x += list[i][0];
sum_y += list[i][1];
num_points += 1;
}
return [sum_x / num_points, sum_y / num_points];
}
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
[
"five_faces.jpg",
"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
training_values.json
{
"000001": [
0,
100
],
"000002": [
0,
100
],
"000005": [
0,
100
],
"000006": [
50,
100
],
"000007": [
25,
0
],
"000009": [
25,
100
],
"000010": [
25,
100
],
"000013": [
25,
0
],
"000014": [
75,
100
],
"000015": [
25,
0
],
"000016": [
25,
0
],
"000018": [
25,
100
],
"000020": [
25,
0
],
"000023": [
25,
0
],
"000025": [
25,
0
],
"000028": [
25,
100
],
"000029": [
25,
100
],
"000030": [
25,
0
],
"000031": [
0,
100
],
"000032": [
0,
0
],
"000035": [
0,
100
],
"000037": [
50,
0
],
"000038": [
0,
0
],
"000040": [
0,
90
],
"000041": [
0,
0
],
"000042": [
0,
100
],
"000043": [
0,
100
],
"000044": [
50,
100
],
"000045": [
0,
100
],
"000047": [
50,
100
],
"000048": [
0,
0
],
"000050": [
0,
0
],
"000051": [
0,
0
],
"000052": [
0,
0
],
"000054": [
0,
100
],
"000055": [
0,
0
],
"000056": [
0,
100
],
"000058": [
0,
100
],
"000060": [
100,
0
],
"000064": [
0,
0
],
"000065": [
50,
0
],
"000068": [
0,
0
],
"000069": [
0,
0
],
"000071": [
0,
100
],
"000073": [
0,
100
],
"000076": [
0,
0
],
"000077": [
0,
100
],
"000078": [
50,
100
],
"000079": [
0,
0
],
"000080": [
0,
0
],
"000081": [
0,
0
],
"000083": [
0,
100
],
"000085": [
0,
100
],
"000086": [
25,
100
],
"000088": [
0,
100
],
"000091": [
0,
0
],
"000092": [
0,
100
],
"000096": [
0,
100
],
"000097": [
0,
100
],
"000099": [
25,
100
],
"000100": [
0,
100
],
"000103": [
25,
100
],
"000104": [
25,
0
],
"000106": [
25,
0
],
"000108": [
25,
100
],
"000109": [
0,
0
],
"000110": [
25,
100
],
"000111": [
50,
100
],
"000114": [
0,
0
],
"000115": [
25,
0
],
"000116": [
0,
0
],
"000117": [
100,
100
],
"000118": [
25,
100
],
"000121": [
25,
100
],
"000122": [
0,
100
],
"000125": [
0,
0
],
"000126": [
0,
100
],
"000129": [
0,
0
],
"000131": [
25,
100
],
"000132": [
50,
100
],
"000133": [
0,
100
],
"000134": [
100,
0
],
"000135": [
75,
0
],
"000137": [
0,
0
],
"000140": [
0,
100
],
"000142": [
0,
100
],
"000143": [
0,
0
],
"000145": [
0,
100
],
"000146": [
0,
100
],
"000147": [
0,
100
],
"000148": [
0,
100
],
"000150": [
0,
0
],
"000151": [
50,
70
],
"000152": [
0,
0
],
"000153": [
0,
0
],
"000155": [
25,
100
],
"000156": [
0,
100
],
"000157": [
0,
100
],
"000160": [
0,
0
],
"000161": [
0,
90
]
}
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 ,
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