documentation.html

<!doctype html>
<html>
<head>
<link rel="SHORTCUT ICON" href="favicon.ico">
<link href='http://fonts.googleapis.com/css?family=Lato' rel='stylesheet' type='text/css'>
<link rel="stylesheet" type="text/css" href="resources/style.css">
<style>
table { font-size:14px; }
td { vertical-align:top; }
@media print
{
table { font-size:12px; }
td.navmain { font-size:26px; }
body { margin: 5mm 5mm 5mm 5mm; }
}
</style>
<title>Numeric Javascript: Documentation</title>
</head>
<body>
<!--#include file="resources/header.html" -->
<!--
This allows regression tests to run predictably:
<pre>
IN> numeric.seedrandom.seedrandom('1'); Math.random = numeric.seedrandom.random; Math.random();
OUT> 0.2694
</pre>
-->

<table cellspacing=5 style="border:5px solid black;">
<tr><td colspan=3 align="center" style="font-size:18px;">
    <b>Reference card for the <tt>numeric</tt> module</b>
<tr valign="top"><td valign="top" width="33%">
<table>
<tr><td><b>Function</b><td><b>Description</b>
<tr><td colspan=2><hr>
<tr><td><tt>abs</tt><td>Absolute value
<tr><td><tt>acos</tt><td>Arc-cosine
<tr><td><tt>add</tt><td>Pointwise sum x+y
<tr><td><tt>addeq</tt><td>Pointwise sum x+=y
<tr><td><tt>all</tt><td>All the components of x are true
<tr><td><tt>and</tt><td>Pointwise x &amp;&amp; y
<tr><td><tt>andeq</tt><td>Pointwise x &amp;= y
<tr><td><tt>any</tt><td>One or more of the components of x are true
<tr><td><tt>asin</tt><td>Arc-sine
<tr><td><tt>atan</tt><td>Arc-tangeant
<tr><td><tt>atan2</tt><td>Arc-tangeant (two parameters)
<tr><td><tt>band</tt><td>Pointwise x &amp; y
<tr><td><tt>bench</tt><td>Benchmarking routine
<tr><td><tt>bnot</tt><td>Binary negation ~x
<tr><td><tt>bor</tt><td>Binary or x|y
<tr><td><tt>bxor</tt><td>Binary xor x^y
<tr><td><tt>ccsDim</tt><td>Dimensions of sparse matrix
<tr><td><tt>ccsDot</tt><td>Sparse matrix-matrix product
<tr><td><tt>ccsFull</tt><td>Convert sparse to full
<tr><td><tt>ccsGather</tt><td>Gather entries of sparse matrix
<tr><td><tt>ccsGetBlock</tt><td>Get rows/columns of sparse matrix
<tr><td><tt>ccsLUP</tt><td>Compute LUP decomposition of sparse matrix
<tr><td><tt>ccsLUPSolve</tt><td>Solve Ax=b using LUP decomp
<tr><td><tt>ccsScatter</tt><td>Scatter entries of sparse matrix
<tr><td><tt>ccsSparse</tt><td>Convert from full to sparse
<tr><td><tt>ccsTSolve</tt><td>Solve upper/lower triangular system
<tr><td><tt>ccs&lt;op&gt;</td><td>Supported ops include: add/div/mul/geq/etc...
<tr><td><tt>cLU</tt><td>Coordinate matrix LU decomposition
<tr><td><tt>cLUsolve</tt><td>Coordinate matrix LU solve
<tr><td><tt>cdelsq</tt><td>Coordinate matrix Laplacian
<tr><td><tt>cdotMV</tt><td>Coordinate matrix-vector product
<tr><td><tt>ceil</tt><td>Pointwise Math.ceil(x)
<tr><td><tt>cgrid</tt><td>Coordinate grid for cdelsq
<tr><td><tt>clone</tt><td>Deep copy of Array
<tr><td><tt>cos</tt><td>Pointwise Math.cos(x)
<tr><td><tt>det</tt><td>Determinant
<tr><td><tt>diag</tt><td>Create diagonal matrix
<tr><td><tt>dim</tt><td>Get Array dimensions
<tr><td><tt>div</tt><td>Pointwise x/y
<tr><td><tt>diveq</tt><td>Pointwise x/=y
<tr><td><tt>dopri</tt><td>Numerical integration of ODE using Dormand-Prince RK method. Returns an object Dopri.
<tr><td><tt>Dopri.at</tt><td>Evaluate the ODE solution at a point
</table>
<td valign="top" width="33%">
<table>
<tr><td><b>Function</b><td><b>Description</b>
<tr><td colspan=2><hr>
<tr><td><tt>dot</tt><td>Matrix-Matrix, Matrix-Vector and Vector-Matrix product
<tr><td><tt>eig</tt><td>Eigenvalues and eigenvectors
<tr><td><tt>epsilon</tt><td>2.220446049250313e-16
<tr><td><tt>eq</tt><td>Pointwise comparison x === y
<tr><td><tt>exp</tt><td>Pointwise Math.exp(x)
<tr><td><tt>floor</tt><td>Poinwise Math.floor(x)
<tr><td><tt>geq</tt><td>Pointwise x&gt;=y
<tr><td><tt>getBlock</tt><td>Extract a block from a matrix
<tr><td><tt>getDiag</tt><td>Get the diagonal of a matrix
<tr><td><tt>gt</tt><td>Pointwise x&gt;y
<tr><td><tt>identity</tt><td>Identity matrix
<tr><td><tt>imageURL</tt><td>Encode a matrix as an image URL
<tr><td><tt>inv</tt><td>Matrix inverse
<tr><td><tt>isFinite</tt><td>Pointwise isFinite(x)
<tr><td><tt>isNaN</tt><td>Pointwise isNaN(x)
<tr><td><tt>largeArray</tt><td>Don't prettyPrint Arrays larger than this
<tr><td><tt>leq</tt><td>Pointwise x&lt;=y
<tr><td><tt>linspace</tt><td>Generate evenly spaced values
<tr><td><tt>log</tt><td>Pointwise Math.log(x)
<tr><td><tt>lshift</tt><td>Pointwise x&lt;&lt;y
<tr><td><tt>lshifteq</tt><td>Pointwise x&lt;&lt;=y
<tr><td><tt>lt</tt><td>Pointwise x&lt;y
<tr><td><tt>LU</tt><td>Dense LU decomposition
<tr><td><tt>LUsolve</tt><td>Dense LU solve
<tr><td><tt>mapreduce</tt><td>Make a pointwise map-reduce function
<tr><td><tt>mod</tt><td>Pointwise x%y
<tr><td><tt>modeq</tt><td>Pointwise x%=y
<tr><td><tt>mul</tt><td>Pointwise x*y
<tr><td><tt>neg</tt><td>Pointwise -x
<tr><td><tt>neq</tt><td>Pointwise x!==y
<tr><td><tt>norm2</tt><td>Square root of the sum of the square of the entries of x
<tr><td><tt>norm2Squared</tt><td>Sum of squares of entries of x
<tr><td><tt>norminf</tt><td>Largest modulus entry of x
<tr><td><tt>not</tt><td>Pointwise logical negation !x
<tr><td><tt>or</tt><td>Pointwise logical or x||y
<tr><td><tt>oreq</tt><td>Pointwise x|=y
<tr><td><tt>parseCSV</tt><td>Parse a CSV file into an Array
<tr><td><tt>parseDate</tt><td>Pointwise parseDate(x)
<tr><td><tt>parseFloat</tt><td>Pointwise parseFloat(x)
<tr><td><tt>pointwise</tt><td>Create a pointwise function
<tr><td><tt>pow</tt><td>Pointwise Math.pow(x)
<tr><td><tt>precision</tt><td>Number of digits to prettyPrint
<tr><td><tt>prettyPrint</tt><td>Pretty-prints x
<tr><td><tt>random</tt><td>Create an Array of random numbers
<tr><td><tt>rep</tt><td>Create an Array by duplicating values
</table>
<td valign="top" width="33%">
<table>
<tr><td><b>Function</b><td><b>Description</b>
<tr><td colspan=2><hr>
<tr><td><tt>round</tt><td>Pointwise Math.round(x)
<tr><td><tt>rrshift</tt><td>Pointwise x&gt;&gt;&gt;y
<tr><td><tt>rrshifteq</tt><td>Pointwise x&gt;&gt;&gt;=y
<tr><td><tt>rshift</tt><td>Pointwise x&gt;&gt;y
<tr><td><tt>rshifteq</tt><td>Pointwise x&gt;&gt;=y
<tr><td><tt>same</tt><td>x and y are entrywise identical
<tr><td><tt>seedrandom</tt><td>The seedrandom module
<tr><td><tt>setBlock</tt><td>Set a block of a matrix
<tr><td><tt>sin</tt><td>Pointwise Math.sin(x)
<tr><td><tt>solve</tt><td>Solve Ax=b
<tr><td><tt>solveLP</tt><td>Solve a linear programming problem
<tr><td><tt>solveQP</tt><td>Solve a quadratic programming problem
<tr><td><tt>spline</tt><td>Create a Spline object
<tr><td><tt>Spline.at</tt><td>Evaluate the Spline at a point
<tr><td><tt>Spline.diff</tt><td>Differentiate the Spline
<tr><td><tt>Spline.roots</tt><td>Find all the roots of the Spline
<tr><td><tt>sqrt</tt><td>Pointwise Math.sqrt(x)
<tr><td><tt>sub</tt><td>Pointwise x-y
<tr><td><tt>subeq</tt><td>Pointwise x-=y
<tr><td><tt>sum</tt><td>Sum all the entries of x
<tr><td><tt>svd</tt><td>Singular value decomposition
<tr><td><tt>t</tt><td>Create a tensor type T (may be complex-valued)
<tr><td><tt>T.&lt;numericfun&gt;</tt><td>Supported &lt;numericfun&gt; are: abs, add, cos, diag, div, dot, exp, getBlock, getDiag, inv, log, mul, neg, norm2, setBlock, sin, sub, transpose
<tr><td><tt>T.conj</tt><td>Pointwise complex conjugate
<tr><td><tt>T.fft</tt><td>Fast Fourier transform
<tr><td><tt>T.get</tt><td>Read an entry
<tr><td><tt>T.getRow</tt><td>Get a row
<tr><td><tt>T.getRows</tt><td>Get a range of rows
<tr><td><tt>T.ifft</tt><td>Inverse FFT
<tr><td><tt>T.reciprocal</tt><td>Pointwise 1/z
<tr><td><tt>T.set</tt><td>Set an entry
<tr><td><tt>T.setRow</tt><td>Set a row
<tr><td><tt>T.setRows</tt><td>Set a range of rows
<tr><td><tt>T.transjugate</tt><td>The conjugate-transpose of a matrix
<tr><td><tt>tan</tt><td>Pointwise Math.tan(x)
<tr><td><tt>tensor</tt><td>Tensor product ret[i][j] = x[i]*y[j]
<tr><td><tt>toCSV</tt><td>Make a CSV file
<tr><td><tt>transpose</tt><td>Matrix transpose
<tr><td><tt>uncmin</tt><td>Unconstrained optimization
<tr><td><tt>version</tt><td>Version string for the numeric library
<tr><td><tt>xor</tt><td>Pointwise x^y
<tr><td><tt>xoreq</tt><td>Pointwise x^=y
</table></table>

<br>
<h1>Numerical analysis in Javascript</h1>

<a href="http://www.numericjs.com/">Numeric Javascript</a> is
library that provides many useful functions for numerical
calculations, particularly for linear algebra (vectors and matrices).
You can create vectors and matrices and multiply them:
<pre>
IN> A = [[1,2,3],[4,5,6]];
OUT> [[1,2,3],
      [4,5,6]]
IN> x = [7,8,9]
OUT> [7,8,9]
IN> numeric.dot(A,x);
OUT> [50,122]
</pre>
The example shown above can be executed in the
<a href="http://www.numericjs.com/workshop.php">Javascript Workshop</a> or at any
Javascript prompt. The Workshop provides plotting capabilities:<br>
<img src="resources/workshop.png"><br>
The function <tt>workshop.plot()</tt> is essentially the <a href="http://code.google.com/p/flot/">flot</a>
plotting command.<br><br>

The <tt>numeric</tt> library provides functions that implement most of the usual Javascript
operators for vectors and matrices:
<pre>
IN> x = [7,8,9];
    y = [10,1,2];
    numeric['+'](x,y)
OUT> [17,9,11]
IN> numeric['>'](x,y)
OUT> [false,true,true]
</pre>
These operators can also be called with plain Javascript function names:
<pre>
IN> numeric.add([7,8,9],[10,1,2])
OUT> [17,9,11]
</pre>
You can also use these operators with three or more parameters:
<pre>
IN> numeric.add([1,2],[3,4],[5,6],[7,8])
OUT> [16,20]
</pre>

The function <tt>numeric.inv()</tt> can be used to compute the inverse of an invertible matrix:
<pre>
IN> A = [[1,2,3],[4,5,6],[7,1,9]]
OUT> [[1,2,3],
      [4,5,6],
      [7,1,9]]
IN> Ainv = numeric.inv(A);
OUT> [[-0.9286,0.3571,0.07143],
      [-0.1429,0.2857,-0.1429],
      [0.7381,-0.3095,0.07143]]
</pre>
The function <tt>numeric.prettyPrint()</tt> is used to print most of the examples in this documentation.
It formats objects, arrays and numbers so that they can be understood easily. All output is automatically
formatted using <tt>numeric.prettyPrint()</tt> when in the
<a href="http://www.numericjs.com/workshop.php">Workshop</a>. In order to present the information clearly and
succintly, the function <tt>numeric.prettyPrint()</tt> lays out matrices so that all the numbers align.
Furthermore, numbers are given approximately using the <tt>numeric.precision</tt> variable:
<pre>
IN> numeric.precision = 10; x = 3.141592653589793
OUT> 3.141592654
IN> numeric.precision = 4; x
OUT> 3.142
</pre>
The default precision is 4 digits. In addition to printing approximate numbers,
the function <tt>numeric.prettyPrint()</tt> will replace large arrays with the string <tt>...Large Array...</tt>:
<pre>
IN> numeric.identity(100)
OUT> ...Large Array...
</pre> 
By default, this happens with the Array's length is more than 50. This can be controlled by setting the
variable <tt>numeric.largeArray</tt> to an appropriate value:
<pre>
IN> numeric.largeArray = 2; A = numeric.identity(4)
OUT> ...Large Array...
IN> numeric.largeArray = 50; A
OUT> [[1,0,0,0],
      [0,1,0,0],
      [0,0,1,0],
      [0,0,0,1]]
</pre>
In particular, if you want to print all Arrays regardless of size, set <tt>numeric.largeArray = Infinity</tt>.
<br><br>


<h1>Math Object functions</h1>

The <tt>Math</tt> object functions have also been adapted to work on Arrays as follows:
<pre>
IN> numeric.exp([1,2]);
OUT> [2.718,7.389]
IN> numeric.exp([[1,2],[3,4]])
OUT> [[2.718, 7.389],
      [20.09, 54.6]]
IN> numeric.abs([-2,3])
OUT> [2,3]
IN> numeric.acos([0.1,0.2])
OUT> [1.471,1.369]
IN> numeric.asin([0.1,0.2])
OUT> [0.1002,0.2014]
IN> numeric.atan([1,2])
OUT> [0.7854,1.107]
IN> numeric.atan2([1,2],[3,4])
OUT> [0.3218,0.4636]
IN> numeric.ceil([-2.2,3.3])
OUT> [-2,4]
IN> numeric.floor([-2.2,3.3])
OUT> [-3,3]
IN> numeric.log([1,2])
OUT> [0,0.6931]
IN> numeric.pow([2,3],[0.25,7.1])
OUT> [1.189,2441]
IN> numeric.round([-2.2,3.3])
OUT> [-2,3]
IN> numeric.sin([1,2])
OUT> [0.8415,0.9093]
IN> numeric.sqrt([1,2])
OUT> [1,1.414]
IN> numeric.tan([1,2])
OUT> [1.557,-2.185]
</pre>


<h1>Utility functions</h1>

The function <tt>numeric.dim()</tt> allows you to compute the dimensions of an Array.

<pre>
IN> numeric.dim([1,2])
OUT> [2]
IN> numeric.dim([[1,2,3],[4,5,6]])
OUT> [2,3]
</pre>

You can perform a deep comparison of Arrays using <tt>numeric.same()</tt>:
<pre>
IN> numeric.same([1,2],[1,2])
OUT> true
IN> numeric.same([1,2],[1,2,3])
OUT> false
IN> numeric.same([1,2],[[1],[2]])
OUT> false
IN> numeric.same([[1,2],[3,4]],[[1,2],[3,4]])
OUT> true
IN> numeric.same([[1,2],[3,4]],[[1,2],[3,5]])
OUT> false
IN> numeric.same([[1,2],[2,4]],[[1,2],[3,4]])
OUT> false
</pre>

You can create a multidimensional Array from a given value using <tt>numeric.rep()</tt>
<pre>
IN> numeric.rep([3],5)
OUT> [5,5,5]
IN> numeric.rep([2,3],0)
OUT> [[0,0,0],
      [0,0,0]]
</pre>

You can loop over Arrays as you normally would. However, in order to quickly generate optimized
loops, the <tt>numeric</tt> library provides a few efficient loop-generation mechanisms. For example, the
<tt>numeric.mapreduce()</tt> function can be used to make a function that computes the sum of all the
entries of an Array.

<pre>
IN> sum = numeric.mapreduce('accum += xi','0'); sum([1,2,3])
OUT> 6
IN> sum([[1,2,3],[4,5,6]])
OUT> 21
</pre>

The functions <tt>numeric.any()</tt> and <tt>numeric.all()</tt> allow you to check whether any or all entries
of an Array are boolean true values.
<pre>
IN> numeric.any([false,true])
OUT> true
IN> numeric.any([[0,0,3.14],[0,false,0]])
OUT> true
IN> numeric.any([0,0,false])
OUT> false
IN> numeric.all([false,true])
OUT> false
IN> numeric.all([[1,4,3.14],["no",true,-1]])
OUT> true
IN> numeric.all([0,0,false])
OUT> false
</pre>

You can create a diagonal matrix using <tt>numeric.diag()</tt>
<pre>
IN> numeric.diag([1,2,3])
OUT> [[1,0,0],
      [0,2,0],
      [0,0,3]]
</pre>

The function <tt>numeric.identity()</tt> returns the identity matrix.
<pre>
IN> numeric.identity(3)
OUT> [[1,0,0],
      [0,1,0],
      [0,0,1]]
</pre>

Random Arrays can also be created:
<pre >
IN> numeric.random([2,3])
OUT> [[0.05303,0.1537,0.7280],
      [0.3839,0.08818,0.6316]]
</pre>

You can generate a vector of evenly spaced values:

<pre>
IN> numeric.linspace(1,5);
OUT> [1,2,3,4,5]
IN> numeric.linspace(1,3,5);
OUT> [1,1.5,2,2.5,3]
</pre>

<!--
<pre>
IN> numeric.blockMatrix([[[[1,2],[3,4]],[[5,6],[7,8]]],
                         [[[11,12],[13,14]],[[15,16],[17,18]]]])
OUT> [[ 1, 2, 5, 6],
      [ 3, 4, 7, 8],
      [11,12,15,16],
      [13,14,17,18]]
</pre>
-->


<h1>Arithmetic operations</h1>

The standard arithmetic operations have been vectorized:
<pre>
IN> numeric.addVV([1,2],[3,4])
OUT> [4,6]
IN> numeric.addVS([1,2],3)
OUT> [4,5]
</pre>

There are also polymorphic functions:
<pre>
IN> numeric.add(1,[2,3])
OUT> [3,4]
IN> numeric.add([1,2,3],[4,5,6])
OUT> [5,7,9]
</pre>

The other arithmetic operations are available:
<pre>
IN> numeric.sub([1,2],[3,4])
OUT> [-2,-2]
IN> numeric.mul([1,2],[3,4])
OUT> [3,8]
IN> numeric.div([1,2],[3,4])
OUT> [0.3333,0.5]
</pre>

The in-place operators (such as +=) are also available:
<pre>
IN> v = [1,2,3,4]; numeric.addeq(v,3); v
OUT> [4,5,6,7]
IN> numeric.subeq([1,2,3],[5,3,1])
OUT> [-4,-1,2]
</pre>

Unary operators:
<pre>
IN> numeric.neg([1,-2,3])
OUT> [-1,2,-3]
IN> numeric.isFinite([10,NaN,Infinity])
OUT> [true,false,false]
IN> numeric.isNaN([10,NaN,Infinity])
OUT> [false,true,false]
</pre>


<!--
<pre>
IN> n = 41; A = numeric.random([n,n]); numeric.norm2(numeric.sub(numeric.dotMMsmall(numeric.inv(A),A),numeric.identity(n)))<1e-12
OUT> true
IN> n = 42; A = numeric.random([n,n]); numeric.norm2(numeric.sub(numeric.dotMMsmall(numeric.inv(A),A),numeric.identity(n)))<1e-12
OUT> true
IN> n = 43; A = numeric.random([n,n]); numeric.norm2(numeric.sub(numeric.dotMMsmall(numeric.inv(A),A),numeric.identity(n)))<1e-12
OUT> true
IN> n = 44; A = numeric.random([n,n]); numeric.norm2(numeric.sub(numeric.dotMMbig(numeric.inv(A),A),numeric.identity(n)))<1e-12
OUT> true
IN> n = 45; A = numeric.random([n,n]); numeric.norm2(numeric.sub(numeric.dotMMbig(numeric.inv(A),A),numeric.identity(n)))<1e-12
OUT> true
IN> n = 46; A = numeric.random([n,n]); numeric.norm2(numeric.sub(numeric.dotMMbig(numeric.inv(A),A),numeric.identity(n)))<1e-12
OUT> true
</pre>
-->


<h1>Linear algebra</h1>

Matrix products are implemented in the functions
<tt>numeric.dotVV()</tt>
<tt>numeric.dotVM()</tt>
<tt>numeric.dotMV()</tt>
<tt>numeric.dotMM()</tt>:
<pre>
IN> numeric.dotVV([1,2],[3,4])
OUT> 11
IN> numeric.dotVM([1,2],[[3,4],[5,6]])
OUT> [13,16]
IN> numeric.dotMV([[1,2],[3,4]],[5,6])
OUT> [17,39]
IN> numeric.dotMMbig([[1,2],[3,4]],[[5,6],[7,8]])
OUT> [[19,22],
      [43,50]]
IN> numeric.dotMMsmall([[1,2],[3,4]],[[5,6],[7,8]])
OUT> [[19,22],
      [43,50]]
IN> numeric.dot([1,2,3,4,5,6,7,8,9],[1,2,3,4,5,6,7,8,9])
OUT> 285
</pre>

The function <tt>numeric.dot()</tt> is "polymorphic" and selects the appropriate Matrix product:

<pre>
IN> numeric.dot([1,2,3],[4,5,6])
OUT> 32
IN> numeric.dot([[1,2,3],[4,5,6]],[7,8,9])
OUT> [50,122]
</pre>

Solving the linear problem Ax=b (<a href="https://github.com/yuzeh">Dan Huang</a>):
<pre>
IN> numeric.solve([[1,2],[3,4]],[17,39])
OUT> [5,6]
</pre>
The algorithm is based on the LU decomposition:
<pre>
IN> LU = numeric.LU([[1,2],[3,4]])
OUT> {LU:[[3     ,4     ],
          [0.3333,0.6667]],
       P:[1,1]}
IN> numeric.LUsolve(LU,[17,39])
OUT> [5,6]
</pre>
<!--
Stress testing
<pre>
IN> ns = [5,6,10,16,25,40,41];
    for(j=0;j<ns.length;j++) {
        n = ns[j];
        for(k=0;k<10;++k) {
            A = numeric.random([n,n]);
            x = numeric.random([n]);
            b = numeric.dot(A,x);
            y = numeric.solve(A,b);
            if(!(numeric.norminf(numeric.sub(x,y))<1e-10)) throw new Error(
                "numeric.solve Stress Test:"+numeric.prettyPrint({j:j,k:k,A:A,x:x,b:b,y:y}));
        }
    }
    "numeric.solve Stress Test OK";
OUT> "numeric.solve Stress Test OK"
</pre>
-->

The determinant:
<pre>
IN> numeric.det([[1,2],[3,4]]);
OUT> -2
IN> numeric.det([[6,8,4,2,8,5],[3,5,2,4,9,2],[7,6,8,3,4,5],[5,5,2,8,1,6],[3,2,2,4,2,2],[8,3,2,2,4,1]]);
OUT> -1404
</pre>

The matrix inverse:
<pre>
IN> numeric.inv([[1,2],[3,4]])
OUT> [[   -2,    1],
      [  1.5, -0.5]]
</pre>

The transpose:
<pre>
IN> numeric.transpose([[1,2,3],[4,5,6]])
OUT> [[1,4],
      [2,5],
      [3,6]]
IN> numeric.transpose([[1,2,3,4,5,6,7,8,9,10,11,12]])
OUT> [[ 1],
      [ 2],
      [ 3],
      [ 4],
      [ 5],
      [ 6],
      [ 7],
      [ 8],
      [ 9],
      [10],
      [11],
      [12]]
</pre>

You can compute the 2-norm of an Array, which is the square root of the sum of the squares of the entries.
<pre>
IN> numeric.norm2([1,2])
OUT> 2.236
</pre>

Computing the tensor product of two vectors:
<pre>
IN> numeric.tensor([1,2],[3,4,5])
OUT> [[3,4,5],
      [6,8,10]]
</pre>


<h1>Data manipulation</h1>
There are also some data manipulation functions. You can parse dates:
<pre>
IN> numeric.parseDate(['1/13/2013','2001-5-9, 9:31']);
OUT> [1.358e12,9.894e11]
</pre>
Parse floating point quantities:
<pre>
IN> numeric.parseFloat(['12','0.1'])
OUT> [12,0.1]
</pre>
Parse CSV files:
<pre>
IN> numeric.parseCSV('a,b,c\n1,2.3,.3\n4e6,-5.3e-8,6.28e+4')
OUT> [[     "a",     "b",     "c"],
      [       1,     2.3,     0.3],
      [     4e6, -5.3e-8,   62800]]
IN> numeric.toCSV([[1.23456789123,2],[3,4]])
OUT> "1.23456789123,2
     3,4
     "
</pre>

You can also fetch a URL (a thin wrapper around XMLHttpRequest):
<pre>
IN> numeric.getURL('tools/helloworld.txt').responseText
OUT> "Hello, world!"
</pre>


<h1>Complex linear algebra</h1>
You can also manipulate complex numbers:
<pre>
IN> z = new numeric.T(3,4);
OUT> {x: 3, y: 4}
IN> z.add(5)
OUT> {x: 8, y: 4}
IN> w = new numeric.T(2,8);
OUT> {x: 2, y: 8}
IN> z.add(w)
OUT> {x: 5, y: 12}
IN> z.mul(w)
OUT> {x: -26, y: 32}
IN> z.div(w)
OUT> {x:0.5588,y:-0.2353}
IN> z.sub(w)
OUT> {x:1, y:-4}
</pre>

Complex vectors and matrices can also be handled:
<pre>
IN> z = new numeric.T([1,2],[3,4]);
OUT> {x: [1,2], y: [3,4]}
IN> z.abs()
OUT> {x:[3.162,4.472],y:}
IN> z.conj()
OUT> {x:[1,2],y:[-3,-4]}
IN> z.norm2()
OUT> 5.477
IN> z.exp()
OUT> {x:[-2.691,-4.83],y:[0.3836,-5.592]}
IN> z.cos()
OUT> {x:[-1.528,-2.459],y:[0.1658,-2.745]}
IN> z.sin()
OUT> {x:[0.2178,-2.847],y:[1.163,2.371]}
IN> z.log()
OUT> {x:[1.151,1.498],y:[1.249,1.107]}
</pre>

Complex matrices:
<pre>
IN> A = new numeric.T([[1,2],[3,4]],[[0,1],[2,-1]]);
OUT> {x:[[1, 2],
         [3, 4]],
      y:[[0, 1],
         [2,-1]]}
IN> A.inv();
OUT> {x:[[0.125,0.125],
         [  0.25,    0]],
      y:[[   0.5,-0.25],
         [-0.375,0.125]]}
IN> A.inv().dot(A)
OUT> {x:[[1,         0],
         [0,         1]],
      y:[[0,-2.776e-17],
         [0,         0]]}
IN> A.get([1,1])
OUT> {x: 4, y: -1}
IN> A.transpose()
OUT> { x: [[1, 3],
           [2, 4]],
       y: [[0, 2],
           [1,-1]] }
IN> A.transjugate()
OUT> { x: [[ 1, 3],
           [ 2, 4]],
       y: [[ 0,-2],
           [-1, 1]] }
IN> numeric.T.rep([2,2],new numeric.T(2,3));
OUT> { x: [[2,2],
           [2,2]],
       y: [[3,3],
           [3,3]] }
</pre>


<h1>Eigenvalues</h1>
Eigenvalues:
<pre>
IN> A = [[1,2,5],[3,5,-1],[7,-3,5]];
OUT> [[ 1,  2,  5],
      [ 3,  5, -1],
      [ 7, -3,  5]]
IN> B = numeric.eig(A);
OUT> {lambda:{x:[-4.284,9.027,6.257],y:},
      E:{x:[[ 0.7131,-0.5543,0.4019],
            [-0.2987,-0.2131,0.9139],
            [-0.6342,-0.8046,0.057 ]],
         y:}}
IN> C = B.E.dot(numeric.T.diag(B.lambda)).dot(B.E.inv());
OUT> {x: [[ 1,  2,  5],
          [ 3,  5, -1],
          [ 7, -3,  5]],
      y: }
</pre>
Note that eigenvalues and eigenvectors are returned as complex numbers (type <tt>numeric.T</tt>). This is because
eigenvalues are often complex even when the matrix is real.<br><br>


<h1>Singular value decomposition (Shanti Rao)</h1>

Shanti Rao kindly emailed me an implementation of the Golub and Reisch algorithm:

<pre>
IN> A=[[ 22, 10,  2,  3,  7],
       [ 14,  7, 10,  0,  8],
       [ -1, 13, -1,-11,  3],
       [ -3, -2, 13, -2,  4],
       [  9,  8,  1, -2,  4],
       [  9,  1, -7,  5, -1],
       [  2, -6,  6,  5,  1],
       [  4,  5,  0, -2,  2]];
    numeric.svd(A);
OUT> {U: 
[[    -0.7071,    -0.1581,     0.1768,     0.2494,     0.4625],
 [    -0.5303,    -0.1581,    -0.3536,     0.1556,    -0.4984],
 [    -0.1768,     0.7906,    -0.1768,    -0.1546,     0.3967],
 [ -1.506e-17,    -0.1581,    -0.7071,    -0.3277,        0.1],
 [    -0.3536,     0.1581,  1.954e-15,   -0.07265,    -0.2084],
 [    -0.1768,    -0.1581,     0.5303,    -0.5726,   -0.05555],
 [ -7.109e-18,    -0.4743,    -0.1768,    -0.3142,     0.4959],
 [    -0.1768,     0.1581,  1.915e-15,     -0.592,    -0.2791]],
S: 
[      35.33,         20,       19.6,          0,          0],
V: 
[[    -0.8006,    -0.3162,     0.2887,    -0.4191,          0],
 [    -0.4804,     0.6325,  7.768e-15,     0.4405,     0.4185],
 [    -0.1601,    -0.3162,     -0.866,     -0.052,     0.3488],
 [  4.684e-17,    -0.6325,     0.2887,     0.6761,     0.2442],
 [    -0.3203,  3.594e-15,    -0.2887,      0.413,    -0.8022]]}
</pre>

<!--
  Some further tests.
<pre>
IN> n = 31; A = numeric.random([n,n]); B = numeric.eig(A); !(B.E.dot(numeric.T.diag(B.lambda).dot(B.E.inv())).sub(A).norm2()>1e-12)
OUT> true
IN> n = 32; A = numeric.random([n,n]); B = numeric.eig(A); !(B.E.dot(numeric.T.diag(B.lambda).dot(B.E.inv())).sub(A).norm2()>1e-12)
OUT> true
IN> n = 33; A = numeric.random([n,n]); B = numeric.eig(A); !(B.E.dot(numeric.T.diag(B.lambda).dot(B.E.inv())).sub(A).norm2()>1e-12)
OUT> true
IN> n = 34; A = numeric.random([n,n]); B = numeric.eig(A); !(B.E.dot(numeric.T.diag(B.lambda).dot(B.E.inv())).sub(A).norm2()>1e-12)
OUT> true
IN> n = 35; A = numeric.random([n,n]); B = numeric.eig(A); !(B.E.dot(numeric.T.diag(B.lambda).dot(B.E.inv())).sub(A).norm2()>1e-12)
OUT> true
IN> m = 17; n = 12; A = numeric.random([m,n]); B = numeric.svd(A); U = new numeric.T(B.U); V = new numeric.T(B.V); !(U.dot(numeric.T.diag(B.S)).dot(V.transpose()).sub(A).norm2()>1e-12)
OUT> true
IN> m = 21; n = 19; A = numeric.random([m,n]); B = numeric.svd(A); U = new numeric.T(B.U); V = new numeric.T(B.V); !(U.dot(numeric.T.diag(B.S)).dot(V.transpose()).sub(A).norm2()>1e-12)
OUT> true
IN> m = 33; n = 33; A = numeric.random([m,n]); B = numeric.svd(A); U = new numeric.T(B.U); V = new numeric.T(B.V); !(U.dot(numeric.T.diag(B.S)).dot(V.transpose()).sub(A).norm2()>1e-12)
OUT> true
IN> m = 59; n = 42; A = numeric.random([m,n]); B = numeric.svd(A); U = new numeric.T(B.U); V = new numeric.T(B.V); !(U.dot(numeric.T.diag(B.S)).dot(V.transpose()).sub(A).norm2()>1e-12)
OUT> true
IN> numeric.eig([[1, 0, 0], [0, 0.7181, -0.6960], [0, 0.6960, 0.7181]]) // This was a bug found by bdmartin
OUT> {lambda: 
{x: 
[          1,     0.7181,     0.7181],
y: 
[          0,      0.696,     -0.696]},
E: 
{x: 
[[          1,          0,          0],
 [          0,          0,     0.7071],
 [          0,    -0.7071,          0]],
y: 
[[          0,          0,          0],
 [          0,    -0.7071,          0],
 [          0,          0,     0.7071]]}}
</pre>
-->


<h1>Sparse linear algebra</h1>

Sparse matrices are matrices that have a lot of zeroes. In numeric, sparse matrices are stored in the
"compressed column storage" ordering. Example:
<pre>
IN> A = [[1,2,0],
         [0,3,0],
         [2,0,5]];
    SA = numeric.ccsSparse(A);
OUT> [[0,2,4,5],
      [0,2,0,1,2],
      [1,2,2,3,5]]
</pre>
The relation between A and its sparse representation SA is:
<pre >
    A[i][SA[1][k]] = SA[2][k] with SA[0][i] &le; k &lt; SA[0][i+1]
</pre >
In other words, SA[2] stores the nonzero entries of A; SA[1] stores the row numbers and SA[0] stores the
offsets of the columns. See <i>I. DUFF, R. GRIMES, AND J. LEWIS, Sparse matrix test problems, ACM Trans. Math. Soft., 15 (1989), pp. 1-14.</i>
<pre>
IN> A = numeric.ccsSparse([[ 3, 5, 8,10, 8],
                          [ 7,10, 3, 5, 3],
                          [ 6, 3, 5, 1, 8],
                          [ 2, 6, 7, 1, 2],
                          [ 1, 2, 9, 3, 9]]);
OUT> [[0,5,10,15,20,25],
      [0,1,2,3,4,0,1,2,3,4,0,1,2,3,4,0,1,2,3,4,0,1,2,3,4],
      [3,7,6,2,1,5,10,3,6,2,8,3,5,7,9,10,5,1,1,3,8,3,8,2,9]]
IN> numeric.ccsFull(A);
OUT> [[ 3, 5, 8,10, 8],
      [ 7,10, 3, 5, 3],
      [ 6, 3, 5, 1, 8],
      [ 2, 6, 7, 1, 2],
      [ 1, 2, 9, 3, 9]]
IN> numeric.ccsDot(numeric.ccsSparse([[1,2,3],[4,5,6]]),numeric.ccsSparse([[7,8],[9,10],[11,12]]))
OUT> [[0,2,4],
      [0,1,0,1],
      [58,139,64,154]]
IN> M = [[0,1,3,6],[0,0,1,0,1,2],[3,-1,2,3,-2,4]];
    b = [9,3,2];
    x = numeric.ccsTSolve(M,b);
OUT> [3.167,2,0.5]
IN> numeric.ccsDot(M,[[0,3],[0,1,2],x])
OUT> [[0,3],[0,1,2],[9,3,2]]
</pre>
We provide an LU=PA decomposition:
<pre>
IN> A = [[0,5,10,15,20,25],
         [0,1,2,3,4,0,1,2,3,4,0,1,2,3,4,0,1,2,3,4,0,1,2,3,4],
         [3,7,6,2,1,5,10,3,6,2,8,3,5,7,9,10,5,1,1,3,8,3,8,2,9]];
    LUP = numeric.ccsLUP(A);
OUT> {L:[[0,5,9,12,14,15],
         [0,2,4,1,3,1,3,4,2,2,4,3,3,4,4],
         [1,0.1429,0.2857,0.8571,0.4286,1,-0.1282,-0.5641,-0.1026,1,0.8517,0.7965,1,-0.67,1]],
      U:[[0,1,3,6,10,15],
         [0,0,1,0,1,2,0,1,2,3,0,1,2,3,4],
         [7,10,-5.571,3,2.429,8.821,5,-3.286,1.949,5.884,3,5.429,9.128,0.1395,-3.476]],
      P:[1,2,4,0,3],
      Pinv:[3,0,1,4,2]}
IN> numeric.ccsFull(numeric.ccsDot(LUP.L,LUP.U))
OUT> [[ 7,10, 3, 5, 3],
      [ 6, 3, 5, 1, 8],
      [ 1, 2, 9, 3, 9],
      [ 3, 5, 8,10, 8],
      [ 2, 6, 7, 1, 2]]
IN> x = numeric.ccsLUPSolve(LUP,[96,63,82,51,89])
OUT> [3,1,4,1,5]
IN> X = numeric.trunc(numeric.ccsFull(numeric.ccsLUPSolve(LUP,A)),1e-15); // Solve LUX = PA
OUT> [[1,0,0,0,0],
      [0,1,0,0,0],
      [0,0,1,0,0],
      [0,0,0,1,0],
      [0,0,0,0,1]]
IN> numeric.ccsLUP(A,0.4).P;
OUT> [0,2,1,3,4]
</pre>
The LUP decomposition uses partial pivoting and has an optional thresholding argument.
With a threshold of 0.4, the pivots are [0,2,1,3,4] (only rows 1 and 2 have been exchanged) instead of the
"full partial pivoting" order above which was [1,2,4,0,3]. Threshold=0 gives no pivoting
and threshold=1 gives normal partial pivoting. Note that a small or zero threshold can result in numerical
instabilities and is normally used when the matrix A is already in some order that minimizes fill-in.
<!-- Stress test:
<pre>
IN> result = "Sparse LUP Stress Test OK";
    for(k=0;k<1000;++k) {
      A = numeric.ccsSparse(numeric.random([10,10]));
      LUP = numeric.ccsLUP(A);
      foo = numeric.ccsFull(numeric.ccsDot(LUP.L,LUP.U));
      PA = numeric.ccsFull(numeric.ccsGetBlock(A,LUP.P));
      res = numeric.norminf(numeric.sub(foo,PA));
      if(!isFinite(res) || res>1e-6) { 
        result = { 
          code: "Failed during 1000 sparse LUP",
          k:k,A:A,LUP:LUP,res:res
        };
        break;
      };
    };
    result;
OUT> "Sparse LUP Stress Test OK"
</pre>
-->

We also support arithmetic on CCS matrices:
<pre>
IN> A = numeric.ccsSparse([[1,2,0],[0,3,0],[0,0,5]]);
    B = numeric.ccsSparse([[2,9,0],[0,4,0],[-2,0,0]]);
    numeric.ccsadd(A,B);
OUT> [[0,2,4,5],
      [0,2,0,1,2],
      [3,-2,11,7,5]]
</pre>

We also have scatter/gather functions
<pre>
IN> X = [[0,0,1,1,2,2],[0,1,1,2,2,3],[1,2,3,4,5,6]];
    SX = numeric.ccsScatter(X);
OUT> [[0,1,3,5,6],
      [0,0,1,1,2,2],
      [1,2,3,4,5,6]]
IN> numeric.ccsGather(SX)
OUT> [[0,0,1,1,2,2],[0,1,1,2,2,3],[1,2,3,4,5,6]]
</pre>


<h1>Coordinate matrices</h1>

We also provide a banded matrix implementation using the coordinate encoding.<br><br>

LU decomposition:
<pre>
IN> lu = numeric.cLU([[0,0,1,1,1,2,2],[0,1,0,1,2,1,2],[2,-1,-1,2,-1,-1,2]])
OUT> {U:[[    0,    0,    1,      1,  2   ],
         [    0,    1,    1,      2,  2   ],
         [    2,   -1,  1.5,     -1, 1.333]],
      L:[[    0,    1,    1,      2,  2   ],
         [    0,    0,    1,      1,  2   ],
         [    1, -0.5,    1,-0.6667,  1   ]]}
IN> numeric.cLUsolve(lu,[5,-8,13])
OUT> [3,1,7]
</pre>
Note that <tt>numeric.cLU()</tt> does not have any pivoting.


<h1>Solving PDEs</h1>

The functions <tt>numeric.cgrid()</tt> and <tt>numeric.cdelsq()</tt> can be used to obtain a
numerical Laplacian for a domain.

<pre>
IN> g = numeric.cgrid(5)
OUT> 
[[-1,-1,-1,-1,-1],
 [-1, 0, 1, 2,-1],
 [-1, 3, 4, 5,-1],
 [-1, 6, 7, 8,-1],
 [-1,-1,-1,-1,-1]]
IN> coordL = numeric.cdelsq(g)
OUT> 
[[ 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8],
 [ 1, 3, 0, 0, 2, 4, 1, 1, 5, 2, 0, 4, 6, 3, 1, 3, 5, 7, 4, 2, 4, 8, 5, 3, 7, 6, 4, 6, 8, 7, 5, 7, 8],
 [-1,-1, 4,-1,-1,-1, 4,-1,-1, 4,-1,-1,-1, 4,-1,-1,-1,-1, 4,-1,-1,-1, 4,-1,-1, 4,-1,-1,-1, 4,-1,-1, 4]]
IN> L = numeric.sscatter(coordL); // Just to see what it looks like
OUT> 
[[          4,         -1,           ,         -1],
 [         -1,          4,         -1,           ,         -1],
 [           ,         -1,          4,           ,           ,         -1],
 [         -1,           ,           ,          4,         -1,           ,         -1],
 [           ,         -1,           ,         -1,          4,         -1,           ,         -1],
 [           ,           ,         -1,           ,         -1,          4,           ,           ,         -1],
 [           ,           ,           ,         -1,           ,           ,          4,         -1],
 [           ,           ,           ,           ,         -1,           ,         -1,          4,         -1],
 [           ,           ,           ,           ,           ,         -1,           ,         -1,          4]]
IN> lu = numeric.cLU(coordL); x = numeric.cLUsolve(lu,[1,1,1,1,1,1,1,1,1]);
OUT> [0.6875,0.875,0.6875,0.875,1.125,0.875,0.6875,0.875,0.6875]
IN> numeric.cdotMV(coordL,x)
OUT> [1,1,1,1,1,1,1,1,1]
IN> G = numeric.rep([5,5],0); for(i=0;i<5;i++) for(j=0;j<5;j++) if(g[i][j]>=0) G[i][j] = x[g[i][j]]; G
OUT> 
[[ 0     , 0     , 0     , 0     , 0     ],
 [ 0     , 0.6875, 0.875 , 0.6875, 0     ],
 [ 0     , 0.875 , 1.125 , 0.875 , 0     ],
 [ 0     , 0.6875, 0.875 , 0.6875, 0     ],
 [ 0     , 0     , 0     , 0     , 0     ]]
IN> workshop.html('&lt;img src="'+numeric.imageURL(numeric.mul([G,G,G],200))+'" width=100 /&gt;');
OUT> 
<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAUAAAAFCAIAAAACDbGyAAAAcElEQVQIHQARAO7/AAAAAAAAAAAAAAAAAAAAAAAAEADv/wAAAIqKiq+vr4mJiQAAAAAAEADv/wAAAK+vr+Hh4a+vrwAAAAAAEADv/wAAAIqKiq+vr4qKigAAAAABEADv/wAAAAAAAAAAAAAAAAAAAACRjRFNqL3leAAAAABJRU5ErkJggg==" width=100 />
</pre>

You can also work on an L-shaped or arbitrary-shape domain:
<pre>
IN> numeric.cgrid(6,'L')
OUT> 
[[-1,-1,-1,-1,-1,-1],
 [-1, 0, 1,-1,-1,-1],
 [-1, 2, 3,-1,-1,-1],
 [-1, 4, 5, 6, 7,-1],
 [-1, 8, 9,10,11,-1],
 [-1,-1,-1,-1,-1,-1]]
IN> numeric.cgrid(5,function(i,j) { return i!==2 || j!==2; })
OUT> 
[[-1,-1,-1,-1,-1],
 [-1, 0, 1, 2,-1],
 [-1, 3,-1, 4,-1],
 [-1, 5, 6, 7,-1],
 [-1,-1,-1,-1,-1]]
</pre>


<h1>Cubic splines</h1>

You can do some (natural) cubic spline interpolation:
<pre>
IN> numeric.spline([1,2,3,4,5],[1,2,1,3,2]).at(numeric.linspace(1,5,10))
OUT> [ 1, 1.731, 2.039, 1.604, 1.019, 1.294, 2.364, 3.085, 2.82, 2]
</pre>
For clamped splines:
<pre>
IN> numeric.spline([1,2,3,4,5],[1,2,1,3,2],0,0).at(numeric.linspace(1,5,10))
OUT> [ 1, 1.435, 1.98, 1.669, 1.034, 1.316, 2.442, 3.017, 2.482, 2]
</pre>
For periodic splines:
<pre>
IN> numeric.spline([1,2,3,4],[0.8415,0.04718,-0.8887,0.8415],'periodic').at(numeric.linspace(1,4,10))
OUT> [ 0.8415, 0.9024, 0.5788, 0.04718, -0.5106, -0.8919, -0.8887, -0.3918, 0.3131, 0.8415]
</pre>
Vector splines:
<pre>
IN> numeric.spline([1,2,3],[[0,1],[1,0],[0,1]]).at(1.78)
OUT> [0.9327,0.06728]
</pre>
Spline differentiation:
<pre>
IN> xs = [0,1,2,3]; numeric.spline(xs,numeric.sin(xs)).diff().at(1.5)
OUT> 0.07302
</pre>
Find all the roots:
<pre>
IN> xs = numeric.linspace(0,30,31); ys = numeric.sin(xs); s = numeric.spline(xs,ys).roots();
OUT> [0, 3.142, 6.284, 9.425, 12.57, 15.71, 18.85, 21.99, 25.13, 28.27]
</pre>


<h1>Fast Fourier Transforms</h1>
FFT and IFFT on numeric.T objects:
<pre>
IN> z = (new numeric.T([1,2,3,4,5],[6,7,8,9,10])).fft()
OUT> {x:[15,-5.941,-3.312,-1.688, 0.941],
      y:[40, 0.941,-1.688,-3.312,-5.941]}
IN> z.ifft()
OUT> {x:[1,2,3,4,5],
      y:[6,7,8,9,10]}
</pre>


<h1>Quadratic Programming (Alberto Santini)</h1>

The Quadratic Programming function <tt>numeric.solveQP()</tt> is based on <a href="https://github.com/albertosantini/node-quadprog">Alberto Santini's
quadprog</a>, which is itself a port of the corresponding
R routines.

<pre>
IN> numeric.solveQP([[1,0,0],[0,1,0],[0,0,1]],[0,5,0],[[-4,2,0],[-3,1,-2],[0,0,1]],[-8,2,0]);
OUT> { solution:              [0.4762,1.048,2.095],
       value:                 [-2.381],
       unconstrained_solution:[     0,    5,    0],
       iterations:            [     3,    0],
       iact:                  [     3,    2,    0],
       message:               "" }
</pre>


<h1>Unconstrained optimization</h1>

To minimize a function f(x) we provide the function <tt>numeric.uncmin(f,x0)</tt> where x0
is a starting point for the optimization.
Here are some demos from from Mor&eacute; et al., 1981:

<pre>
IN> sqr = function(x) { return x*x; }; 
    numeric.uncmin(function(x) { return sqr(10*(x[1]-x[0]*x[0])) + sqr(1-x[0]); },[-1.2,1]).solution
OUT> [1,1]
IN> f = function(x) { return sqr(-13+x[0]+((5-x[1])*x[1]-2)*x[1])+sqr(-29+x[0]+((x[1]+1)*x[1]-14)*x[1]); }; 
    x0 = numeric.uncmin(f,[0.5,-2]).solution
OUT> [11.41,-0.8968]
IN> f = function(x) { return sqr(1e4*x[0]*x[1]-1)+sqr(Math.exp(-x[0])+Math.exp(-x[1])-1.0001); }; 
    x0 = numeric.uncmin(f,[0,1]).solution
OUT> [1.098e-5,9.106]
IN> f = function(x) { return sqr(x[0]-1e6)+sqr(x[1]-2e-6)+sqr(x[0]*x[1]-2)}; 
    x0 = numeric.uncmin(f,[0,1]).solution
OUT> [1e6,2e-6]
IN> f = function(x) { 
       return sqr(1.5-x[0]*(1-x[1]))+sqr(2.25-x[0]*(1-x[1]*x[1]))+sqr(2.625-x[0]*(1-x[1]*x[1]*x[1]));
    }; 
    x0 = numeric.uncmin(f,[1,1]).solution
OUT> [3,0.5]
IN> f = function(x) { 
        var ret = 0,i; 
        for(i=1;i<=10;i++) ret+=sqr(2+2*i-Math.exp(i*x[0])-Math.exp(i*x[1]));
         return ret; 
    };
    x0 = numeric.uncmin(f,[0.3,0.4]).solution
OUT> [0.2578,0.2578]
IN> y = [0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,0.37,0.58,0.73,0.96,1.34,2.10,4.39]; 
    f = function(x) { 
        var ret = 0,i; 
        for(i=1;i<=15;i++) ret+=sqr(y[i-1]-(x[0]+i/((16-i)*x[1]+Math.min(i,16-i)*x[2]))); 
        return ret; 
    }; 
    x0 = numeric.uncmin(f,[1,1,1]).solution
OUT> [0.08241,1.133,2.344]
IN> y = [0.0009,0.0044,0.0175,0.0540,0.1295,0.2420,0.3521,0.3989,0.3521,0.2420,0.1295,0.0540,0.0175,0.0044,0.0009]; 
    f = function(x) { 
        var ret = 0,i; 
        for(i=1;i<=15;i++) 
        ret+=sqr(x[0]*Math.exp(-x[1]*sqr((8-i)/2-x[2])/2)-y[i-1]); 
        return ret; 
    }; 
    x0 = numeric.div(numeric.round(numeric.mul(numeric.uncmin(f,[1,1,1]).solution,1000)),1000)
OUT> [0.399,1,0]
IN> f = function(x) { return sqr(x[0]+10*x[1])+5*sqr(x[2]-x[3])+sqr(sqr(x[1]-2*x[2]))+10*sqr(x[0]-x[3]); }; 
    x0 = numeric.div(numeric.round(numeric.mul(numeric.uncmin(f,[3,-1,0,1]).solution,1e5)),1e5)
OUT> [0,0,0,0]
IN> f = function(x) { 
        return (sqr(10*(x[1]-x[0]*x[0]))+sqr(1-x[0])+
                90*sqr(x[3]-x[2]*x[2])+sqr(1-x[2])+
                10*sqr(x[1]+x[3]-2)+0.1*sqr(x[1]-x[3])); }; 
    x0 = numeric.uncmin(f,[-3,-1,-3,-1]).solution
OUT> [1,1,1,1]
IN> y = [0.1957,0.1947,0.1735,0.1600,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246]; 
    u = [4,2,1,0.5,0.25,0.167,0.125,0.1,0.0833,0.0714,0.0625]; 
    f = function(x) { 
        var ret=0, i; 
        for(i=0;i<11;++i) ret += sqr(y[i]-x[0]*(u[i]*u[i]+u[i]*x[1])/(u[i]*u[i]+u[i]*x[2]+x[3])); 
        return ret; 
    }; 
    x0 = numeric.uncmin(f,[0.25,0.39,0.415,0.39]).solution
OUT> [     0.1928,     0.1913,     0.1231,     0.1361]
IN> y = [0.844,0.908,0.932,0.936,0.925,0.908,0.881,0.850,0.818,0.784,0.751,0.718,
         0.685,0.658,0.628,0.603,0.580,0.558,0.538,0.522,0.506,0.490,0.478,0.467,
         0.457,0.448,0.438,0.431,0.424,0.420,0.414,0.411,0.406]; 
    f = function(x) { 
        var ret=0, i; 
        for(i=0;i<33;++i) ret += sqr(y[i]-(x[0]+x[1]*Math.exp(-10*i*x[3])+x[2]*Math.exp(-10*i*x[4]))); 
        return ret; 
    }; 
    x0 = numeric.uncmin(f,[0.5,1.5,-1,0.01,0.02]).solution
OUT> [     0.3754,      1.936,     -1.465,    0.01287,    0.02212]
IN> f = function(x) { 
        var ret=0, i,ti,yi,exp=Math.exp; 
        for(i=1;i<=13;++i) { 
            ti = 0.1*i; 
            yi = exp(-ti)-5*exp(-10*ti)+3*exp(-4*ti); 
            ret += sqr(x[2]*exp(-ti*x[0])-x[3]*exp(-ti*x[1])+x[5]*exp(-ti*x[4])-yi); 
        } 
        return ret; 
    }; 
    x0 = numeric.uncmin(f,[1,2,1,1,1,1],1e-14).solution;
    f(x0)<1e-20;
OUT> true
</pre>
There are optional parameters to <tt>numeric.uncmin(f,x0,tol,gradient,maxit,callback)</tt>. The iteration stops when
the gradient or step size is less than the optional parameter <tt>tol</tt>. The <tt>gradient()</tt> parameter is a function that computes the
gradient of <tt>f()</tt>. If it is not provided, a numerical gradient is used. The iteration stops when
<tt>maxit</tt> iterations have been performed. The optional <tt>callback()</tt> parameter, if provided, is called at each step:
<pre>
IN> z = []; 
    cb = function(i,x,f,g,H) { z.push({i:i, x:x, f:f, g:g, H:H }); }; 
    x0 = numeric.uncmin(function(x) { return Math.cos(2*x[0]); },
                        [1],1e-10,
                        function(x) { return [-2*Math.sin(2*x[0])]; },
                        100,cb);
OUT> {solution:    [1.571],
      f:           -1,
      gradient:    [2.242e-11],
      invHessian:  [[0.25]],
      iterations:  6,
      message:     "Newton step smaller than tol"}
IN> z
OUT> [{i:0, x:[1    ], f:-0.4161, g: [-1.819   ]  , H:[[1     ]] },
      {i:2, x:[1.909], f:-0.7795, g: [ 1.253   ]  , H:[[0.296 ]] },
      {i:3, x:[1.538], f:-0.9979, g: [-0.1296  ]  , H:[[0.2683]] },
      {i:4, x:[1.573], f:-1     , g: [ 9.392e-3]  , H:[[0.2502]] },
      {i:5, x:[1.571], f:-1     , g: [-6.105e-6]  , H:[[0.25  ]] },
      {i:6, x:[1.571], f:-1     , g: [ 2.242e-11] , H:[[0.25  ]] }]
</pre>


<h1>Linear programming</h1>

Linear programming is available:

<pre>
IN> x = numeric.solveLP([1,1],                   /* minimize [1,1]*x                */
                        [[-1,0],[0,-1],[-1,-2]], /* matrix of inequalities          */
                        [0,0,-3]                 /* right-hand-side of inequalities */
                        );       
    numeric.trunc(x.solution,1e-12);
OUT> [0,1.5]
</pre>

The function <tt>numeric.solveLP(c,A,b)</tt> minimizes dot(c,x) subject to dot(A,x) <= b.
The algorithm used is very basic. For alpha>0, define the ``barrier function''
f0 = dot(c,x) - alpha*sum(log(b-A*x)). The function numeric.solveLP calls numeric.uncmin
on f0 for smaller and smaller values of alpha. This is a basic ``interior point method''.

We also handle infeasible problems:
<pre>
IN> numeric.solveLP([1,1],[[1,0],[0,1],[-1,-1]],[-1,-1,-1])
OUT> { solution: NaN, message: "Infeasible", iterations: 5 }
</pre>

Unbounded problems:
<pre>
IN> numeric.solveLP([1,1],[[1,0],[0,1]],[0,0]).message;
OUT> "Unbounded"
</pre>

With an equality constraint:
<pre>
IN> numeric.solveLP([1,2,3],                      /* minimize [1,2,3]*x                */
                    [[-1,0,0],[0,-1,0],[0,0,-1]], /* matrix A of inequality constraint */
                    [0,0,0],                      /* RHS b of inequality constraint    */
                    [[1,1,1]],                    /* matrix Aeq of equality constraint */
                    [3]                           /* vector beq of equality constraint */
                    );
OUT> { solution:[3,1.685e-16,4.559e-19], message:"", iterations:12 }
</pre>

<!--
<pre>
IN> n = 7; m = 3;
    for(k=0;k<10;++k) {
        A = numeric.random([n,n]);
        x = numeric.rep([m],1).concat(numeric.rep([n-m],0));
        b = numeric.dot(A,x);
        J = numeric.diag(numeric.rep([n],-1));
        B = numeric.blockMatrix([[A                   , J     ],
                                 [numeric.neg(A)      , J     ],
                                 [numeric.rep([n,n],0), J     ]]);
        c = b.concat(numeric.neg(b)).concat(numeric.rep([n],0));
        d = numeric.rep([n],0).concat(numeric.rep([n],1));
        y = numeric.solveLP(d,B,c).solution;
        y.length = n;
        foo = numeric.norm2(numeric.sub(x,y));
        if(foo>1e-10) throw new Error("solveLP test fails: "+numeric.prettyPrint({A:A,x:x}));
    }
    "solveLP tests pass"
OUT> "solveLP tests pass"
</pre>
-->

<!-- 
Bug found by Michael J.
<pre>
IN> numeric.solveLP([1,2,3], /* minimize [1,2,3]*x */
    [[1, 0, 0], [0, 1, 0], [0, 0, 1],[-1,0,0],[0,-1,0],[0,0,-1]], /* matrix A of inequality constraint */
    [1,1,1,0,0,0], /* RHS b of inequality constraint */
    [[1,1,1]], /* matrix Aeq of equality constraint */
    [3] /* vector beq of equality constraint */
    );
OUT> { solution:NaN, message:"Infeasible", iterations:10 }
</pre>
-->

<h1>Solving ODEs</h1>

The function <tt>numeric.dopri()</tt> is an implementation of the Dormand-Prince-Runge-Kutta integrator with
adaptive time-stepping:
<pre>
IN> sol = numeric.dopri(0,1,1,function(t,y) { return y; })
OUT> { x:    [     0,   0.1, 0.1522, 0.361, 0.5792, 0.7843, 0.9813,     1],
       y:    [     1, 1.105,  1.164, 1.435,  1.785,  2.191,  2.668, 2.718],
       f:    [     1, 1.105,  1.164, 1.435,  1.785,  2.191,  2.668, 2.718],
       ymid: [ 1.051, 1.134,  1.293,   1.6,  1.977,  2.418,  2.693],
       iterations:8,
       events:,
       message:""}
IN> sol.at([0.3,0.7])
OUT> [1.35,2.014]
</pre>
The return value <tt>sol</tt> contains the x and y values of the solution.
If you need to know the value of the solution between the given x values, use the function
<tt>sol.at()</tt>, which uses the extra information contained in <tt>sol.ymid</tt> and <tt>sol.f</tt> to
produce approximations between these points.
The integrator is also able to handle vector equations. This is a harmonic oscillator:
<pre>
IN> sol = numeric.dopri(0,10,[3,0],function (x,y) { return [y[1],-y[0]]; });
    sol.at([0,0.5*Math.PI,Math.PI,1.5*Math.PI,2*Math.PI])
OUT> [[         3,         0],
      [ -9.534e-8,        -3],
      [        -3,  2.768e-7],
      [   3.63e-7,         3],
      [         3, -3.065e-7]]
</pre>
Van der Pol:
<pre>
IN> numeric.dopri(0,20,[2,0],function(t,y) { return [y[1], (1-y[0]*y[0])*y[1]-y[0]]; }).at([18,19,20])
OUT> [[ -1.208,   0.9916],
      [ 0.4258,    2.535],
      [  2.008, -0.04251]]
</pre>
You can also specify a tolerance, a maximum number of iterations and an event function. The integration stops if
the event function goes from negative to positive.
<pre>
IN> sol = numeric.dopri(0,2,1,function (x,y) { return y; },1e-8,100,function (x,y) { return y-1.3; });
OUT> { x:          [     0, 0.0181, 0.09051, 0.1822,  0.2624],
       y:          [     1,  1.018,   1.095,    1.2,     1.3],
       f:          [     1,  1.018,   1.095,    1.2,     1.3],
       ymid:       [ 1.009,  1.056,   1.146,  1.249],
       iterations: 5,
       events:     true,
       message:    "" }
</pre>
Note that at <tt>x=0.1822</tt>, the event function value was <tt>-0.1001</tt> while at <tt>x=0.2673</tt>, the event
value was <tt>6.433e-3</tt>. The integrator thus terminated at <tt>x=0.2673</tt> instead of continuing until the end
of the integration interval.<br><br>

Events can also be vector-valued:
<pre>
IN> sol = numeric.dopri(0,2,1,
                        function(x,y) { return y; },
                        undefined,50,
                        function(x,y) { return [y-1.5,Math.sin(y-1.5)]; });
OUT> { x:          [    0,   0.2, 0.4055],
       y:          [    1, 1.221,    1.5],
       f:          [    1, 1.221,    1.5],
       ymid:       [1.105, 1.354],
       iterations: 2,
       events:     [true,true],
       message:    ""}
</pre>

<!--
<pre>
IN> D = numeric.identity(8); d = numeric.rep([8],0); A = [[1, 1, -1, 0,  0, 0,  0, 0,  0, 0],[-1, 1,  0, 1,  0, 0,  0, 0,  0, 0],[1, 1,  0, 0, -1, 0,  0, 0,  0, 0],[-1, 1,  0, 0,  0, 1,  0, 0,  0, 0],[1, 1,  0, 0,  0, 0, -1, 0,  0, 0],[-1, 1,  0, 0,  0, 0,  0, 1,  0, 0],[1, 1,  0, 0,  0, 0,  0, 0, -1, 0],[-1, 1,  0, 0,  0, 0,  0, 0,  0, 1]]; b = [1,  1, -1,  0, -1,  0, -1,  0, -1,  0]; numeric.solveQP(D,d,A,b,undefined,2)
OUT> { solution:              [0.25,0,0.25,0,0.25,0,0.25,0],
       value:                 [0.125],
       unconstrained_solution:[0,0,0,0,0,0,0,0],
       iterations:            [3,0],
       iact:                  [1,2,0,0,0,0,0,0,0,0],
       message:               ""}
IN> numeric.imageURL(numeric.rep([3],[[1,2],[3,4]]));
OUT> "data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAIAAAACCAIAAAD91JpzAAAAH0lEQVQIHQAIAPf/AAEBAQICAgABBwD4/wMDAwQEBAAAogAfhs2H3QAAAABJRU5ErkJggg=="
IN> numeric.gradient(function(x) { return Math.exp(Math.sin(x[0])*(x[1]+2)); },[1,2]);
OUT> [62.59,24.37]
</pre>
-->


<h1>Seedrandom (David Bau)</h1>

The object <tt>numeric.seedrandom</tt> is based on 
<a href="http://davidbau.com/archives/2010/01/30/random_seeds_coded_hints_and_quintillions.html">David Bau's <tt>seedrandom.js</tt></a>.
This small library can be used to create better pseudorandom numbers than <tt>Math.random()</tt> which can
furthermore be "seeded".
<pre>
IN> numeric.seedrandom.seedrandom(3); numeric.seedrandom.random()
OUT> 0.7569
IN> numeric.seedrandom.random()
OUT> 0.6139
IN> numeric.seedrandom.seedrandom(3); numeric.seedrandom.random()
OUT> 0.7569
</pre>
For performance reasons, <tt>numeric.random()</tt> uses the default <tt>Math.random()</tt>. If you want to use
the seedrandom version, just do <tt>Math.random = numeric.seedrandom.random</tt>. Note that this may slightly
decrease the performance of all <tt>Math</tt> operations.


<br><br><br>