Remove JAMA example files
This commit is contained in:
Mark Baker
parent
1f601e0ecf
commit
f2803e8690
|
|
@ -1,116 +0,0 @@
|
||||||
<?php
|
|
||||||
/**
|
|
||||||
* quadratic (p-o)'S'S(p-o)
|
|
||||||
* solve for o, S
|
|
||||||
* S is a single scale factor
|
|
||||||
*/
|
|
||||||
class LMQuadTest {
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @param array[] $x
|
|
||||||
* @param array[] $a
|
|
||||||
*/
|
|
||||||
function val($x, $a) {
|
|
||||||
if (count($a) != 3) die ("Wrong number of elements in array a");
|
|
||||||
if (count($x) != 2) die ("Wrong number of elements in array x");
|
|
||||||
|
|
||||||
$ox = $a[0];
|
|
||||||
$oy = $a[1];
|
|
||||||
$s = $a[2];
|
|
||||||
|
|
||||||
$sdx = $s * ($x[0] - $ox);
|
|
||||||
$sdy = $s * ($x[1] - $oy);
|
|
||||||
|
|
||||||
return ($sdx * $sdx) + ($sdy * $sdy);
|
|
||||||
} // function val()
|
|
||||||
|
|
||||||
|
|
||||||
/**
|
|
||||||
* z = (p-o)'S'S(p-o)
|
|
||||||
* dz/dp = 2S'S(p-o)
|
|
||||||
*
|
|
||||||
* z = (s*(px-ox))^2 + (s*(py-oy))^2
|
|
||||||
* dz/dox = -2(s*(px-ox))*s
|
|
||||||
* dz/ds = 2*s*[(px-ox)^2 + (py-oy)^2]
|
|
||||||
*
|
|
||||||
* z = (s*dx)^2 + (s*dy)^2
|
|
||||||
* dz/ds = 2(s*dx)*dx + 2(s*dy)*dy
|
|
||||||
*
|
|
||||||
* @param array[] $x
|
|
||||||
* @param array[] $a
|
|
||||||
* @param int $a_k
|
|
||||||
* @param array[] $a
|
|
||||||
*/
|
|
||||||
function grad($x, $a, $a_k) {
|
|
||||||
if (count($a) != 3) die ("Wrong number of elements in array a");
|
|
||||||
if (count($x) != 2) die ("Wrong number of elements in array x");
|
|
||||||
if ($a_k < 3) die ("a_k=".$a_k);
|
|
||||||
|
|
||||||
$ox = $a[0];
|
|
||||||
$oy = $a[1];
|
|
||||||
$s = $a[2];
|
|
||||||
|
|
||||||
$dx = ($x[0] - $ox);
|
|
||||||
$dy = ($x[1] - $oy);
|
|
||||||
|
|
||||||
if ($a_k == 0)
|
|
||||||
return -2.*$s*$s*$dx;
|
|
||||||
elseif ($a_k == 1)
|
|
||||||
return -2.*$s*$s*$dy;
|
|
||||||
else
|
|
||||||
return 2.*$s*($dx*$dx + $dy*$dy);
|
|
||||||
} // function grad()
|
|
||||||
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @return array[] $a
|
|
||||||
*/
|
|
||||||
function initial() {
|
|
||||||
$a[0] = 0.05;
|
|
||||||
$a[1] = 0.1;
|
|
||||||
$a[2] = 1.0;
|
|
||||||
|
|
||||||
return $a;
|
|
||||||
} // function initial()
|
|
||||||
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @return Object[] $a
|
|
||||||
*/
|
|
||||||
function testdata() {
|
|
||||||
$npts = 25;
|
|
||||||
|
|
||||||
$a[0] = 0.;
|
|
||||||
$a[1] = 0.;
|
|
||||||
$a[2] = 0.9;
|
|
||||||
|
|
||||||
$i = 0;
|
|
||||||
|
|
||||||
for ($r = -2; $r <= 2; ++$r) {
|
|
||||||
for ($c = -2; $c <= 2; ++$c) {
|
|
||||||
$x[$i][0] = $c;
|
|
||||||
$x[$i][1] = $r;
|
|
||||||
$y[$i] = $this->val($x[$i], $a);
|
|
||||||
print("Quad ".$c.",".$r." -> ".$y[$i]."<br />");
|
|
||||||
$s[$i] = 1.;
|
|
||||||
++$i;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
print("quad x= ");
|
|
||||||
|
|
||||||
$qx = new Matrix($x);
|
|
||||||
$qx->print(10, 2);
|
|
||||||
|
|
||||||
print("quad y= ");
|
|
||||||
$qy = new Matrix($y, $npts);
|
|
||||||
$qy->print(10, 2);
|
|
||||||
|
|
||||||
$o[0] = $x;
|
|
||||||
$o[1] = $a;
|
|
||||||
$o[2] = $y;
|
|
||||||
$o[3] = $s;
|
|
||||||
|
|
||||||
return $o;
|
|
||||||
} // function testdata()
|
|
||||||
|
|
||||||
} // class LMQuadTest
|
|
||||||
|
|
@ -1,59 +0,0 @@
|
||||||
<?php
|
|
||||||
|
|
||||||
require_once "../Matrix.php";
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Given n points (x0,y0)...(xn-1,yn-1), the following methid computes
|
|
||||||
* the polynomial factors of the n-1't degree polynomial passing through
|
|
||||||
* the n points.
|
|
||||||
*
|
|
||||||
* Example: Passing in three points (2,3) (1,4) and (3,7) will produce
|
|
||||||
* the results [2.5, -8.5, 10] which means that the points are on the
|
|
||||||
* curve y = 2.5x² - 8.5x + 10.
|
|
||||||
*
|
|
||||||
* @see http://geosoft.no/software/lagrange/LagrangeInterpolation.java.html
|
|
||||||
* @author Jacob Dreyer
|
|
||||||
* @author Paul Meagher (port to PHP and minor changes)
|
|
||||||
*
|
|
||||||
* @param x[] float
|
|
||||||
* @param y[] float
|
|
||||||
*/
|
|
||||||
class LagrangeInterpolation {
|
|
||||||
|
|
||||||
public function findPolynomialFactors($x, $y) {
|
|
||||||
$n = count($x);
|
|
||||||
|
|
||||||
$data = array(); // double[n][n];
|
|
||||||
$rhs = array(); // double[n];
|
|
||||||
|
|
||||||
for ($i = 0; $i < $n; ++$i) {
|
|
||||||
$v = 1;
|
|
||||||
for ($j = 0; $j < $n; ++$j) {
|
|
||||||
$data[$i][$n-$j-1] = $v;
|
|
||||||
$v *= $x[$i];
|
|
||||||
}
|
|
||||||
$rhs[$i] = $y[$i];
|
|
||||||
}
|
|
||||||
|
|
||||||
// Solve m * s = b
|
|
||||||
$m = new Matrix($data);
|
|
||||||
$b = new Matrix($rhs, $n);
|
|
||||||
|
|
||||||
$s = $m->solve($b);
|
|
||||||
|
|
||||||
return $s->getRowPackedCopy();
|
|
||||||
} // function findPolynomialFactors()
|
|
||||||
|
|
||||||
} // class LagrangeInterpolation
|
|
||||||
|
|
||||||
|
|
||||||
$x = array(2.0, 1.0, 3.0);
|
|
||||||
$y = array(3.0, 4.0, 7.0);
|
|
||||||
|
|
||||||
$li = new LagrangeInterpolation;
|
|
||||||
$f = $li->findPolynomialFactors($x, $y);
|
|
||||||
|
|
||||||
|
|
||||||
for ($i = 0; $i < 3; ++$i) {
|
|
||||||
echo $f[$i]."<br />";
|
|
||||||
}
|
|
||||||
|
|
@ -1,59 +0,0 @@
|
||||||
<?php
|
|
||||||
|
|
||||||
require_once "../Matrix.php";
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Given n points (x0,y0)...(xn-1,yn-1), the following method computes
|
|
||||||
* the polynomial factors of the n-1't degree polynomial passing through
|
|
||||||
* the n points.
|
|
||||||
*
|
|
||||||
* Example: Passing in three points (2,3) (1,4) and (3,7) will produce
|
|
||||||
* the results [2.5, -8.5, 10] which means that the points are on the
|
|
||||||
* curve y = 2.5x² - 8.5x + 10.
|
|
||||||
*
|
|
||||||
* @see http://geosoft.no/software/lagrange/LagrangeInterpolation.java.html
|
|
||||||
* @see http://source.freehep.org/jcvsweb/ilc/LCSIM/wdview/lcsim/src/org/lcsim/fit/polynomial/PolynomialFitter.java
|
|
||||||
* @author Jacob Dreyer
|
|
||||||
* @author Paul Meagher (port to PHP and minor changes)
|
|
||||||
*
|
|
||||||
* @param x[] float
|
|
||||||
* @param y[] float
|
|
||||||
*/
|
|
||||||
class LagrangeInterpolation {
|
|
||||||
|
|
||||||
public function findPolynomialFactors($x, $y) {
|
|
||||||
$n = count($x);
|
|
||||||
|
|
||||||
$data = array(); // double[n][n];
|
|
||||||
$rhs = array(); // double[n];
|
|
||||||
|
|
||||||
for ($i = 0; $i < $n; ++$i) {
|
|
||||||
$v = 1;
|
|
||||||
for ($j = 0; $j < $n; ++$j) {
|
|
||||||
$data[$i][$n-$j-1] = $v;
|
|
||||||
$v *= $x[$i];
|
|
||||||
}
|
|
||||||
$rhs[$i] = $y[$i];
|
|
||||||
}
|
|
||||||
|
|
||||||
// Solve m * s = b
|
|
||||||
$m = new Matrix($data);
|
|
||||||
$b = new Matrix($rhs, $n);
|
|
||||||
|
|
||||||
$s = $m->solve($b);
|
|
||||||
|
|
||||||
return $s->getRowPackedCopy();
|
|
||||||
} // function findPolynomialFactors()
|
|
||||||
|
|
||||||
} // class LagrangeInterpolation
|
|
||||||
|
|
||||||
|
|
||||||
$x = array(2.0, 1.0, 3.0);
|
|
||||||
$y = array(3.0, 4.0, 7.0);
|
|
||||||
|
|
||||||
$li = new LagrangeInterpolation;
|
|
||||||
$f = $li->findPolynomialFactors($x, $y);
|
|
||||||
|
|
||||||
for ($i = 0; $i < 3; ++$i) {
|
|
||||||
echo $f[$i]."<br />";
|
|
||||||
}
|
|
||||||
|
|
@ -1,185 +0,0 @@
|
||||||
<?php
|
|
||||||
|
|
||||||
// Levenberg-Marquardt in PHP
|
|
||||||
|
|
||||||
// http://www.idiom.com/~zilla/Computer/Javanumeric/LM.java
|
|
||||||
|
|
||||||
class LevenbergMarquardt {
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Calculate the current sum-squared-error
|
|
||||||
*
|
|
||||||
* Chi-squared is the distribution of squared Gaussian errors,
|
|
||||||
* thus the name.
|
|
||||||
*
|
|
||||||
* @param double[][] $x
|
|
||||||
* @param double[] $a
|
|
||||||
* @param double[] $y,
|
|
||||||
* @param double[] $s,
|
|
||||||
* @param object $f
|
|
||||||
*/
|
|
||||||
function chiSquared($x, $a, $y, $s, $f) {
|
|
||||||
$npts = count($y);
|
|
||||||
$sum = 0.0;
|
|
||||||
|
|
||||||
for ($i = 0; $i < $npts; ++$i) {
|
|
||||||
$d = $y[$i] - $f->val($x[$i], $a);
|
|
||||||
$d = $d / $s[$i];
|
|
||||||
$sum = $sum + ($d*$d);
|
|
||||||
}
|
|
||||||
|
|
||||||
return $sum;
|
|
||||||
} // function chiSquared()
|
|
||||||
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Minimize E = sum {(y[k] - f(x[k],a)) / s[k]}^2
|
|
||||||
* The individual errors are optionally scaled by s[k].
|
|
||||||
* Note that LMfunc implements the value and gradient of f(x,a),
|
|
||||||
* NOT the value and gradient of E with respect to a!
|
|
||||||
*
|
|
||||||
* @param x array of domain points, each may be multidimensional
|
|
||||||
* @param y corresponding array of values
|
|
||||||
* @param a the parameters/state of the model
|
|
||||||
* @param vary false to indicate the corresponding a[k] is to be held fixed
|
|
||||||
* @param s2 sigma^2 for point i
|
|
||||||
* @param lambda blend between steepest descent (lambda high) and
|
|
||||||
* jump to bottom of quadratic (lambda zero).
|
|
||||||
* Start with 0.001.
|
|
||||||
* @param termepsilon termination accuracy (0.01)
|
|
||||||
* @param maxiter stop and return after this many iterations if not done
|
|
||||||
* @param verbose set to zero (no prints), 1, 2
|
|
||||||
*
|
|
||||||
* @return the new lambda for future iterations.
|
|
||||||
* Can use this and maxiter to interleave the LM descent with some other
|
|
||||||
* task, setting maxiter to something small.
|
|
||||||
*/
|
|
||||||
function solve($x, $a, $y, $s, $vary, $f, $lambda, $termepsilon, $maxiter, $verbose) {
|
|
||||||
$npts = count($y);
|
|
||||||
$nparm = count($a);
|
|
||||||
|
|
||||||
if ($verbose > 0) {
|
|
||||||
print("solve x[".count($x)."][".count($x[0])."]");
|
|
||||||
print(" a[".count($a)."]");
|
|
||||||
println(" y[".count(length)."]");
|
|
||||||
}
|
|
||||||
|
|
||||||
$e0 = $this->chiSquared($x, $a, $y, $s, $f);
|
|
||||||
|
|
||||||
//double lambda = 0.001;
|
|
||||||
$done = false;
|
|
||||||
|
|
||||||
// g = gradient, H = hessian, d = step to minimum
|
|
||||||
// H d = -g, solve for d
|
|
||||||
$H = array();
|
|
||||||
$g = array();
|
|
||||||
|
|
||||||
//double[] d = new double[nparm];
|
|
||||||
|
|
||||||
$oos2 = array();
|
|
||||||
|
|
||||||
for($i = 0; $i < $npts; ++$i) {
|
|
||||||
$oos2[$i] = 1./($s[$i]*$s[$i]);
|
|
||||||
}
|
|
||||||
$iter = 0;
|
|
||||||
$term = 0; // termination count test
|
|
||||||
|
|
||||||
do {
|
|
||||||
++$iter;
|
|
||||||
|
|
||||||
// hessian approximation
|
|
||||||
for( $r = 0; $r < $nparm; ++$r) {
|
|
||||||
for( $c = 0; $c < $nparm; ++$c) {
|
|
||||||
for( $i = 0; $i < $npts; ++$i) {
|
|
||||||
if ($i == 0) $H[$r][$c] = 0.;
|
|
||||||
$xi = $x[$i];
|
|
||||||
$H[$r][$c] += ($oos2[$i] * $f->grad($xi, $a, $r) * $f->grad($xi, $a, $c));
|
|
||||||
} //npts
|
|
||||||
} //c
|
|
||||||
} //r
|
|
||||||
|
|
||||||
// boost diagonal towards gradient descent
|
|
||||||
for( $r = 0; $r < $nparm; ++$r)
|
|
||||||
$H[$r][$r] *= (1. + $lambda);
|
|
||||||
|
|
||||||
// gradient
|
|
||||||
for( $r = 0; $r < $nparm; ++$r) {
|
|
||||||
for( $i = 0; $i < $npts; ++$i) {
|
|
||||||
if ($i == 0) $g[$r] = 0.;
|
|
||||||
$xi = $x[$i];
|
|
||||||
$g[$r] += ($oos2[$i] * ($y[$i]-$f->val($xi,$a)) * $f->grad($xi, $a, $r));
|
|
||||||
}
|
|
||||||
} //npts
|
|
||||||
|
|
||||||
// scale (for consistency with NR, not necessary)
|
|
||||||
if ($false) {
|
|
||||||
for( $r = 0; $r < $nparm; ++$r) {
|
|
||||||
$g[$r] = -0.5 * $g[$r];
|
|
||||||
for( $c = 0; $c < $nparm; ++$c) {
|
|
||||||
$H[$r][$c] *= 0.5;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// solve H d = -g, evaluate error at new location
|
|
||||||
//double[] d = DoubleMatrix.solve(H, g);
|
|
||||||
// double[] d = (new Matrix(H)).lu().solve(new Matrix(g, nparm)).getRowPackedCopy();
|
|
||||||
//double[] na = DoubleVector.add(a, d);
|
|
||||||
// double[] na = (new Matrix(a, nparm)).plus(new Matrix(d, nparm)).getRowPackedCopy();
|
|
||||||
// double e1 = chiSquared(x, na, y, s, f);
|
|
||||||
|
|
||||||
// if (verbose > 0) {
|
|
||||||
// System.out.println("\n\niteration "+iter+" lambda = "+lambda);
|
|
||||||
// System.out.print("a = ");
|
|
||||||
// (new Matrix(a, nparm)).print(10, 2);
|
|
||||||
// if (verbose > 1) {
|
|
||||||
// System.out.print("H = ");
|
|
||||||
// (new Matrix(H)).print(10, 2);
|
|
||||||
// System.out.print("g = ");
|
|
||||||
// (new Matrix(g, nparm)).print(10, 2);
|
|
||||||
// System.out.print("d = ");
|
|
||||||
// (new Matrix(d, nparm)).print(10, 2);
|
|
||||||
// }
|
|
||||||
// System.out.print("e0 = " + e0 + ": ");
|
|
||||||
// System.out.print("moved from ");
|
|
||||||
// (new Matrix(a, nparm)).print(10, 2);
|
|
||||||
// System.out.print("e1 = " + e1 + ": ");
|
|
||||||
// if (e1 < e0) {
|
|
||||||
// System.out.print("to ");
|
|
||||||
// (new Matrix(na, nparm)).print(10, 2);
|
|
||||||
// } else {
|
|
||||||
// System.out.println("move rejected");
|
|
||||||
// }
|
|
||||||
// }
|
|
||||||
|
|
||||||
// termination test (slightly different than NR)
|
|
||||||
// if (Math.abs(e1-e0) > termepsilon) {
|
|
||||||
// term = 0;
|
|
||||||
// } else {
|
|
||||||
// term++;
|
|
||||||
// if (term == 4) {
|
|
||||||
// System.out.println("terminating after " + iter + " iterations");
|
|
||||||
// done = true;
|
|
||||||
// }
|
|
||||||
// }
|
|
||||||
// if (iter >= maxiter) done = true;
|
|
||||||
|
|
||||||
// in the C++ version, found that changing this to e1 >= e0
|
|
||||||
// was not a good idea. See comment there.
|
|
||||||
//
|
|
||||||
// if (e1 > e0 || Double.isNaN(e1)) { // new location worse than before
|
|
||||||
// lambda *= 10.;
|
|
||||||
// } else { // new location better, accept new parameters
|
|
||||||
// lambda *= 0.1;
|
|
||||||
// e0 = e1;
|
|
||||||
// // simply assigning a = na will not get results copied back to caller
|
|
||||||
// for( int i = 0; i < nparm; i++ ) {
|
|
||||||
// if (vary[i]) a[i] = na[i];
|
|
||||||
// }
|
|
||||||
// }
|
|
||||||
} while(!$done);
|
|
||||||
|
|
||||||
return $lambda;
|
|
||||||
} // function solve()
|
|
||||||
|
|
||||||
} // class LevenbergMarquardt
|
|
||||||
|
|
@ -1,182 +0,0 @@
|
||||||
<?php
|
|
||||||
/**
|
|
||||||
* @package JAMA
|
|
||||||
*/
|
|
||||||
|
|
||||||
require_once "../Matrix.php";
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Example of use of Matrix Class, featuring magic squares.
|
|
||||||
*/
|
|
||||||
class MagicSquareExample {
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Generate magic square test matrix.
|
|
||||||
* @param int n dimension of matrix
|
|
||||||
*/
|
|
||||||
function magic($n) {
|
|
||||||
|
|
||||||
// Odd order
|
|
||||||
|
|
||||||
if (($n % 2) == 1) {
|
|
||||||
$a = ($n+1)/2;
|
|
||||||
$b = ($n+1);
|
|
||||||
for ($j = 0; $j < $n; ++$j)
|
|
||||||
for ($i = 0; $i < $n; ++$i)
|
|
||||||
$M[$i][$j] = $n*(($i+$j+$a) % $n) + (($i+2*$j+$b) % $n) + 1;
|
|
||||||
|
|
||||||
// Doubly Even Order
|
|
||||||
|
|
||||||
} else if (($n % 4) == 0) {
|
|
||||||
for ($j = 0; $j < $n; ++$j) {
|
|
||||||
for ($i = 0; $i < $n; ++$i) {
|
|
||||||
if ((($i+1)/2)%2 == (($j+1)/2)%2)
|
|
||||||
$M[$i][$j] = $n*$n-$n*$i-$j;
|
|
||||||
else
|
|
||||||
$M[$i][$j] = $n*$i+$j+1;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// Singly Even Order
|
|
||||||
|
|
||||||
} else {
|
|
||||||
|
|
||||||
$p = $n/2;
|
|
||||||
$k = ($n-2)/4;
|
|
||||||
$A = $this->magic($p);
|
|
||||||
$M = array();
|
|
||||||
for ($j = 0; $j < $p; ++$j) {
|
|
||||||
for ($i = 0; $i < $p; ++$i) {
|
|
||||||
$aij = $A->get($i,$j);
|
|
||||||
$M[$i][$j] = $aij;
|
|
||||||
$M[$i][$j+$p] = $aij + 2*$p*$p;
|
|
||||||
$M[$i+$p][$j] = $aij + 3*$p*$p;
|
|
||||||
$M[$i+$p][$j+$p] = $aij + $p*$p;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
for ($i = 0; $i < $p; ++$i) {
|
|
||||||
for ($j = 0; $j < $k; ++$j) {
|
|
||||||
$t = $M[$i][$j];
|
|
||||||
$M[$i][$j] = $M[$i+$p][$j];
|
|
||||||
$M[$i+$p][$j] = $t;
|
|
||||||
}
|
|
||||||
for ($j = $n-$k+1; $j < $n; ++$j) {
|
|
||||||
$t = $M[$i][$j];
|
|
||||||
$M[$i][$j] = $M[$i+$p][$j];
|
|
||||||
$M[$i+$p][$j] = $t;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
$t = $M[$k][0]; $M[$k][0] = $M[$k+$p][0]; $M[$k+$p][0] = $t;
|
|
||||||
$t = $M[$k][$k]; $M[$k][$k] = $M[$k+$p][$k]; $M[$k+$p][$k] = $t;
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
return new Matrix($M);
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Simple function to replicate PHP 5 behaviour
|
|
||||||
*/
|
|
||||||
function microtime_float() {
|
|
||||||
list($usec, $sec) = explode(" ", microtime());
|
|
||||||
return ((float)$usec + (float)$sec);
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Tests LU, QR, SVD and symmetric Eig decompositions.
|
|
||||||
*
|
|
||||||
* n = order of magic square.
|
|
||||||
* trace = diagonal sum, should be the magic sum, (n^3 + n)/2.
|
|
||||||
* max_eig = maximum eigenvalue of (A + A')/2, should equal trace.
|
|
||||||
* rank = linear algebraic rank, should equal n if n is odd,
|
|
||||||
* be less than n if n is even.
|
|
||||||
* cond = L_2 condition number, ratio of singular values.
|
|
||||||
* lu_res = test of LU factorization, norm1(L*U-A(p,:))/(n*eps).
|
|
||||||
* qr_res = test of QR factorization, norm1(Q*R-A)/(n*eps).
|
|
||||||
*/
|
|
||||||
function main() {
|
|
||||||
?>
|
|
||||||
<p>Test of Matrix Class, using magic squares.</p>
|
|
||||||
<p>See MagicSquareExample.main() for an explanation.</p>
|
|
||||||
<table border='1' cellspacing='0' cellpadding='4'>
|
|
||||||
<tr>
|
|
||||||
<th>n</th>
|
|
||||||
<th>trace</th>
|
|
||||||
<th>max_eig</th>
|
|
||||||
<th>rank</th>
|
|
||||||
<th>cond</th>
|
|
||||||
<th>lu_res</th>
|
|
||||||
<th>qr_res</th>
|
|
||||||
</tr>
|
|
||||||
<?php
|
|
||||||
$start_time = $this->microtime_float();
|
|
||||||
$eps = pow(2.0,-52.0);
|
|
||||||
for ($n = 3; $n <= 6; ++$n) {
|
|
||||||
echo "<tr>";
|
|
||||||
|
|
||||||
echo "<td align='right'>$n</td>";
|
|
||||||
|
|
||||||
$M = $this->magic($n);
|
|
||||||
$t = (int) $M->trace();
|
|
||||||
|
|
||||||
echo "<td align='right'>$t</td>";
|
|
||||||
|
|
||||||
$O = $M->plus($M->transpose());
|
|
||||||
$E = new EigenvalueDecomposition($O->times(0.5));
|
|
||||||
$d = $E->getRealEigenvalues();
|
|
||||||
|
|
||||||
echo "<td align='right'>".$d[$n-1]."</td>";
|
|
||||||
|
|
||||||
$r = $M->rank();
|
|
||||||
|
|
||||||
echo "<td align='right'>".$r."</td>";
|
|
||||||
|
|
||||||
$c = $M->cond();
|
|
||||||
|
|
||||||
if ($c < 1/$eps)
|
|
||||||
echo "<td align='right'>".sprintf("%.3f",$c)."</td>";
|
|
||||||
else
|
|
||||||
echo "<td align='right'>Inf</td>";
|
|
||||||
|
|
||||||
$LU = new LUDecomposition($M);
|
|
||||||
$L = $LU->getL();
|
|
||||||
$U = $LU->getU();
|
|
||||||
$p = $LU->getPivot();
|
|
||||||
// Java version: R = L.times(U).minus(M.getMatrix(p,0,n-1));
|
|
||||||
$S = $L->times($U);
|
|
||||||
$R = $S->minus($M->getMatrix($p,0,$n-1));
|
|
||||||
$res = $R->norm1()/($n*$eps);
|
|
||||||
|
|
||||||
echo "<td align='right'>".sprintf("%.3f",$res)."</td>";
|
|
||||||
|
|
||||||
$QR = new QRDecomposition($M);
|
|
||||||
$Q = $QR->getQ();
|
|
||||||
$R = $QR->getR();
|
|
||||||
$S = $Q->times($R);
|
|
||||||
$R = $S->minus($M);
|
|
||||||
$res = $R->norm1()/($n*$eps);
|
|
||||||
|
|
||||||
echo "<td align='right'>".sprintf("%.3f",$res)."</td>";
|
|
||||||
|
|
||||||
echo "</tr>";
|
|
||||||
|
|
||||||
}
|
|
||||||
echo "<table>";
|
|
||||||
echo "<br />";
|
|
||||||
|
|
||||||
$stop_time = $this->microtime_float();
|
|
||||||
$etime = $stop_time - $start_time;
|
|
||||||
|
|
||||||
echo "<p>Elapsed time is ". sprintf("%.4f",$etime) ." seconds.</p>";
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
$magic = new MagicSquareExample();
|
|
||||||
$magic->main();
|
|
||||||
|
|
||||||
?>
|
|
||||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,263 +0,0 @@
|
||||||
<?php
|
|
||||||
|
|
||||||
error_reporting(E_ALL);
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @package JAMA
|
|
||||||
*/
|
|
||||||
|
|
||||||
require_once '../Matrix.php';
|
|
||||||
require_once 'Stats.php';
|
|
||||||
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Example of use of Matrix Class, featuring magic squares.
|
|
||||||
*/
|
|
||||||
class Benchmark {
|
|
||||||
public $stat;
|
|
||||||
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Simple function to replicate PHP 5 behaviour
|
|
||||||
*/
|
|
||||||
function microtime_float() {
|
|
||||||
list($usec, $sec) = explode(" ", microtime());
|
|
||||||
|
|
||||||
return ((float)$usec + (float)$sec);
|
|
||||||
} // function microtime_float()
|
|
||||||
|
|
||||||
|
|
||||||
function displayStats($times = null) {
|
|
||||||
$this->stat->setData($times);
|
|
||||||
$stats = $this->stat->calcFull();
|
|
||||||
|
|
||||||
echo '<table style="margin-left:32px;">';
|
|
||||||
echo '<tr><td style="text-align:right;"><b>n:</b><td style="text-align:right;">' . $stats['count'] . ' </td></tr>';
|
|
||||||
echo '<tr><td style="text-align:right;"><b>Mean:</b><td style="text-align:right;">' . $stats['mean'] . ' </td></tr>';
|
|
||||||
echo '<tr><td style="text-align:right;"><b>Min.:</b><td style="text-align:right;">' . $stats['min'] . ' </td></tr>';
|
|
||||||
echo '<tr><td style="text-align:right;"><b>Max.:</b><td style="text-align:right;">' . $stats['max'] . ' </td></tr>';
|
|
||||||
echo '<tr><td style="text-align:right;"><b>σ:</b><td style="text-align:right;">' . $stats['stdev'] . ' </td></tr>';
|
|
||||||
echo '<tr><td style="text-align:right;"><b>Variance:</b><td style="text-align:right;">' . $stats['variance'] . ' </td></tr>';
|
|
||||||
echo '<tr><td style="text-align:right;"><b>Range:</b><td style="text-align:right;">' . $stats['range'] . ' </td></tr>';
|
|
||||||
echo '</table>';
|
|
||||||
|
|
||||||
return $stats;
|
|
||||||
} // function displayStats()
|
|
||||||
|
|
||||||
|
|
||||||
function runEig($n = 4, $t = 100) {
|
|
||||||
$times = array();
|
|
||||||
|
|
||||||
for ($i = 0; $i < $t; ++$i) {
|
|
||||||
$M = Matrix::random($n, $n);
|
|
||||||
$start_time = $this->microtime_float();
|
|
||||||
$E = new EigenvalueDecomposition($M);
|
|
||||||
$stop_time = $this->microtime_float();
|
|
||||||
$times[] = $stop_time - $start_time;
|
|
||||||
}
|
|
||||||
|
|
||||||
return $times;
|
|
||||||
} // function runEig()
|
|
||||||
|
|
||||||
|
|
||||||
function runLU($n = 4, $t = 100) {
|
|
||||||
$times = array();
|
|
||||||
|
|
||||||
for ($i = 0; $i < $t; ++$i) {
|
|
||||||
$M = Matrix::random($n, $n);
|
|
||||||
$start_time = $this->microtime_float();
|
|
||||||
$E = new LUDecomposition($M);
|
|
||||||
$stop_time = $this->microtime_float();
|
|
||||||
$times[] = $stop_time - $start_time;
|
|
||||||
}
|
|
||||||
|
|
||||||
return $times;
|
|
||||||
} // function runLU()
|
|
||||||
|
|
||||||
|
|
||||||
function runQR($n = 4, $t = 100) {
|
|
||||||
$times = array();
|
|
||||||
|
|
||||||
for ($i = 0; $i < $t; ++$i) {
|
|
||||||
$M = Matrix::random($n, $n);
|
|
||||||
$start_time = $this->microtime_float();
|
|
||||||
$E = new QRDecomposition($M);
|
|
||||||
$stop_time = $this->microtime_float();
|
|
||||||
$times[] = $stop_time - $start_time;
|
|
||||||
}
|
|
||||||
|
|
||||||
return $times;
|
|
||||||
} // function runQR()
|
|
||||||
|
|
||||||
|
|
||||||
function runCholesky($n = 4, $t = 100) {
|
|
||||||
$times = array();
|
|
||||||
|
|
||||||
for ($i = 0; $i < $t; ++$i) {
|
|
||||||
$M = Matrix::random($n, $n);
|
|
||||||
$start_time = $this->microtime_float();
|
|
||||||
$E = new CholeskyDecomposition($M);
|
|
||||||
$stop_time = $this->microtime_float();
|
|
||||||
$times[] = $stop_time - $start_time;
|
|
||||||
}
|
|
||||||
|
|
||||||
return $times;
|
|
||||||
} // function runCholesky()
|
|
||||||
|
|
||||||
|
|
||||||
function runSVD($n = 4, $t = 100) {
|
|
||||||
$times = array();
|
|
||||||
|
|
||||||
for ($i = 0; $i < $t; ++$i) {
|
|
||||||
$M = Matrix::random($n, $n);
|
|
||||||
$start_time = $this->microtime_float();
|
|
||||||
$E = new SingularValueDecomposition($M);
|
|
||||||
$stop_time = $this->microtime_float();
|
|
||||||
$times[] = $stop_time - $start_time;
|
|
||||||
}
|
|
||||||
|
|
||||||
return $times;
|
|
||||||
} // function runSVD()
|
|
||||||
|
|
||||||
|
|
||||||
function run() {
|
|
||||||
$n = 8;
|
|
||||||
$t = 16;
|
|
||||||
$sum = 0;
|
|
||||||
echo "<b>Cholesky decomposition: $t random {$n}x{$n} matrices</b><br />";
|
|
||||||
$r = $this->displayStats($this->runCholesky($n, $t));
|
|
||||||
$sum += $r['mean'] * $n;
|
|
||||||
|
|
||||||
echo '<hr />';
|
|
||||||
|
|
||||||
echo "<b>Eigenvalue decomposition: $t random {$n}x{$n} matrices</b><br />";
|
|
||||||
$r = $this->displayStats($this->runEig($n, $t));
|
|
||||||
$sum += $r['mean'] * $n;
|
|
||||||
|
|
||||||
echo '<hr />';
|
|
||||||
|
|
||||||
echo "<b>LU decomposition: $t random {$n}x{$n} matrices</b><br />";
|
|
||||||
$r = $this->displayStats($this->runLU($n, $t));
|
|
||||||
$sum += $r['mean'] * $n;
|
|
||||||
|
|
||||||
echo '<hr />';
|
|
||||||
|
|
||||||
echo "<b>QR decomposition: $t random {$n}x{$n} matrices</b><br />";
|
|
||||||
$r = $this->displayStats($this->runQR($n, $t));
|
|
||||||
$sum += $r['mean'] * $n;
|
|
||||||
|
|
||||||
echo '<hr />';
|
|
||||||
|
|
||||||
echo "<b>Singular Value decomposition: $t random {$n}x{$n} matrices</b><br />";
|
|
||||||
$r = $this->displayStats($this->runSVD($n, $t));
|
|
||||||
$sum += $r['mean'] * $n;
|
|
||||||
|
|
||||||
return $sum;
|
|
||||||
} // function run()
|
|
||||||
|
|
||||||
|
|
||||||
public function __construct() {
|
|
||||||
$this->stat = new Base();
|
|
||||||
} // function Benchmark()
|
|
||||||
|
|
||||||
} // class Benchmark (end MagicSquareExample)
|
|
||||||
|
|
||||||
|
|
||||||
$benchmark = new Benchmark();
|
|
||||||
|
|
||||||
switch($_REQUEST['decomposition']) {
|
|
||||||
case 'cholesky':
|
|
||||||
$m = array();
|
|
||||||
for ($i = 2; $i <= 8; $i *= 2) {
|
|
||||||
$t = 32 / $i;
|
|
||||||
echo "<b>Cholesky decomposition: $t random {$i}x{$i} matrices</b><br />";
|
|
||||||
$s = $benchmark->displayStats($benchmark->runCholesky($i, $t));
|
|
||||||
$m[$i] = $s['mean'];
|
|
||||||
echo "<br />";
|
|
||||||
}
|
|
||||||
echo '<pre>';
|
|
||||||
foreach($m as $x => $y) {
|
|
||||||
echo "$x\t" . 1000*$y . "\n";
|
|
||||||
}
|
|
||||||
echo '</pre>';
|
|
||||||
break;
|
|
||||||
case 'eigenvalue':
|
|
||||||
$m = array();
|
|
||||||
for ($i = 2; $i <= 8; $i *= 2) {
|
|
||||||
$t = 32 / $i;
|
|
||||||
echo "<b>Eigenvalue decomposition: $t random {$i}x{$i} matrices</b><br />";
|
|
||||||
$s = $benchmark->displayStats($benchmark->runEig($i, $t));
|
|
||||||
$m[$i] = $s['mean'];
|
|
||||||
echo "<br />";
|
|
||||||
}
|
|
||||||
echo '<pre>';
|
|
||||||
foreach($m as $x => $y) {
|
|
||||||
echo "$x\t" . 1000*$y . "\n";
|
|
||||||
}
|
|
||||||
echo '</pre>';
|
|
||||||
break;
|
|
||||||
case 'lu':
|
|
||||||
$m = array();
|
|
||||||
for ($i = 2; $i <= 8; $i *= 2) {
|
|
||||||
$t = 32 / $i;
|
|
||||||
echo "<b>LU decomposition: $t random {$i}x{$i} matrices</b><br />";
|
|
||||||
$s = $benchmark->displayStats($benchmark->runLU($i, $t));
|
|
||||||
$m[$i] = $s['mean'];
|
|
||||||
echo "<br />";
|
|
||||||
}
|
|
||||||
echo '<pre>';
|
|
||||||
foreach($m as $x => $y) {
|
|
||||||
echo "$x\t" . 1000*$y . "\n";
|
|
||||||
}
|
|
||||||
echo '</pre>';
|
|
||||||
break;
|
|
||||||
case 'qr':
|
|
||||||
$m = array();
|
|
||||||
for ($i = 2; $i <= 8; $i *= 2) {
|
|
||||||
$t = 32 / $i;
|
|
||||||
echo "<b>QR decomposition: $t random {$i}x{$i} matrices</b><br />";
|
|
||||||
$s = $benchmark->displayStats($benchmark->runQR($i, $t));
|
|
||||||
$m[$i] = $s['mean'];
|
|
||||||
echo "<br />";
|
|
||||||
}
|
|
||||||
echo '<pre>';
|
|
||||||
foreach($m as $x => $y) {
|
|
||||||
echo "$x\t" . 1000*$y . "\n";
|
|
||||||
}
|
|
||||||
echo '</pre>';
|
|
||||||
break;
|
|
||||||
case 'svd':
|
|
||||||
$m = array();
|
|
||||||
for($i = 2; $i <= 8; $i *= 2) {
|
|
||||||
$t = 32 / $i;
|
|
||||||
echo "<b>Singular value decomposition: $t random {$i}x{$i} matrices</b><br />";
|
|
||||||
$s = $benchmark->displayStats($benchmark->runSVD($i, $t));
|
|
||||||
$m[$i] = $s['mean'];
|
|
||||||
echo "<br />";
|
|
||||||
}
|
|
||||||
echo '<pre>';
|
|
||||||
foreach($m as $x => $y) {
|
|
||||||
echo "$x\t" . 1000*$y . "\n";
|
|
||||||
}
|
|
||||||
echo '</pre>';
|
|
||||||
break;
|
|
||||||
case 'all':
|
|
||||||
$s = $benchmark->run();
|
|
||||||
print("<br /><b>Total<b>: {$s}s<br />");
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
?>
|
|
||||||
<ul>
|
|
||||||
<li><a href="benchmark.php?decomposition=all">Complete Benchmark</a>
|
|
||||||
<ul>
|
|
||||||
<li><a href="benchmark.php?decomposition=cholesky">Cholesky</a></li>
|
|
||||||
<li><a href="benchmark.php?decomposition=eigenvalue">Eigenvalue</a></li>
|
|
||||||
<li><a href="benchmark.php?decomposition=lu">LU</a></li>
|
|
||||||
<li><a href="benchmark.php?decomposition=qr">QR</a></li>
|
|
||||||
<li><a href="benchmark.php?decomposition=svd">Singular Value</a></li>
|
|
||||||
</ul>
|
|
||||||
</li>
|
|
||||||
</ul>
|
|
||||||
<?php
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
@ -1,73 +0,0 @@
|
||||||
<?php
|
|
||||||
require_once "../Matrix.php";
|
|
||||||
/*
|
|
||||||
* @package JAMA
|
|
||||||
* @author Michael Bommarito
|
|
||||||
* @author Paul Meagher
|
|
||||||
* @version 0.1
|
|
||||||
*
|
|
||||||
* Function to fit an order n polynomial function through
|
|
||||||
* a series of x-y data points using least squares.
|
|
||||||
*
|
|
||||||
* @param $X array x values
|
|
||||||
* @param $Y array y values
|
|
||||||
* @param $n int order of polynomial to be used for fitting
|
|
||||||
* @returns array $coeffs of polynomial coefficients
|
|
||||||
* Pre-Conditions: the system is not underdetermined: sizeof($X) > $n+1
|
|
||||||
*/
|
|
||||||
function polyfit($X, $Y, $n) {
|
|
||||||
for ($i = 0; $i < sizeof($X); ++$i)
|
|
||||||
for ($j = 0; $j <= $n; ++$j)
|
|
||||||
$A[$i][$j] = pow($X[$i], $j);
|
|
||||||
for ($i=0; $i < sizeof($Y); ++$i)
|
|
||||||
$B[$i] = array($Y[$i]);
|
|
||||||
$matrixA = new Matrix($A);
|
|
||||||
$matrixB = new Matrix($B);
|
|
||||||
$C = $matrixA->solve($matrixB);
|
|
||||||
return $C->getMatrix(0, $n, 0, 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
function printpoly( $C = null ) {
|
|
||||||
for($i = $C->m - 1; $i >= 0; --$i) {
|
|
||||||
$r = $C->get($i, 0);
|
|
||||||
if ( abs($r) <= pow(10, -9) )
|
|
||||||
$r = 0;
|
|
||||||
if ($i == $C->m - 1)
|
|
||||||
echo $r . "x<sup>$i</sup>";
|
|
||||||
else if ($i < $C->m - 1)
|
|
||||||
echo " + " . $r . "x<sup>$i</sup>";
|
|
||||||
else if ($i == 0)
|
|
||||||
echo " + " . $r;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
$X = array(0,1,2,3,4,5);
|
|
||||||
$Y = array(4,3,12,67,228, 579);
|
|
||||||
$points = new Matrix(array($X, $Y));
|
|
||||||
$points->toHTML();
|
|
||||||
printpoly(polyfit($X, $Y, 4));
|
|
||||||
|
|
||||||
echo '<hr />';
|
|
||||||
|
|
||||||
$X = array(0,1,2,3,4,5);
|
|
||||||
$Y = array(1,2,5,10,17, 26);
|
|
||||||
$points = new Matrix(array($X, $Y));
|
|
||||||
$points->toHTML();
|
|
||||||
printpoly(polyfit($X, $Y, 2));
|
|
||||||
|
|
||||||
echo '<hr />';
|
|
||||||
|
|
||||||
$X = array(0,1,2,3,4,5,6);
|
|
||||||
$Y = array(-90,-104,-178,-252,-26, 1160, 4446);
|
|
||||||
$points = new Matrix(array($X, $Y));
|
|
||||||
$points->toHTML();
|
|
||||||
printpoly(polyfit($X, $Y, 5));
|
|
||||||
|
|
||||||
echo '<hr />';
|
|
||||||
|
|
||||||
$X = array(0,1,2,3,4);
|
|
||||||
$Y = array(mt_rand(0, 10), mt_rand(40, 80), mt_rand(240, 400), mt_rand(1800, 2215), mt_rand(8000, 9000));
|
|
||||||
$points = new Matrix(array($X, $Y));
|
|
||||||
$points->toHTML();
|
|
||||||
printpoly(polyfit($X, $Y, 3));
|
|
||||||
?>
|
|
||||||
|
|
@ -1,78 +0,0 @@
|
||||||
<?php
|
|
||||||
|
|
||||||
include "../Matrix.php";
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Tiling of matrix X in [rowWise by colWise] dimension. Tiling
|
|
||||||
* creates a larger matrix than the original data X. Example, if
|
|
||||||
* X is to be tiled in a [3 x 4] manner, then:
|
|
||||||
*
|
|
||||||
* / \
|
|
||||||
* | X X X X |
|
|
||||||
* C = | X X X X |
|
|
||||||
* | X X X X |
|
|
||||||
* \ /
|
|
||||||
*
|
|
||||||
* @param X Matrix
|
|
||||||
* @param rowWise int
|
|
||||||
* @param colWise int
|
|
||||||
* @return Matrix
|
|
||||||
*/
|
|
||||||
|
|
||||||
function tile(&$X, $rowWise, $colWise){
|
|
||||||
|
|
||||||
$xArray = $X->getArray();
|
|
||||||
print_r($xArray);
|
|
||||||
|
|
||||||
$countRow = 0;
|
|
||||||
$countColumn = 0;
|
|
||||||
|
|
||||||
$m = $X->getRowDimension();
|
|
||||||
$n = $X->getColumnDimension();
|
|
||||||
|
|
||||||
if( $rowWise<1 || $colWise<1 ){
|
|
||||||
die("tile : Array index is out-of-bound.");
|
|
||||||
}
|
|
||||||
|
|
||||||
$newRowDim = $m*$rowWise;
|
|
||||||
$newColDim = $n*$colWise;
|
|
||||||
|
|
||||||
$result = array();
|
|
||||||
|
|
||||||
for($i=0 ; $i<$newRowDim; ++$i) {
|
|
||||||
|
|
||||||
$holder = array();
|
|
||||||
|
|
||||||
for($j=0 ; $j<$newColDim ; ++$j) {
|
|
||||||
|
|
||||||
$holder[$j] = $xArray[$countRow][$countColumn++];
|
|
||||||
|
|
||||||
// reset the column-index to zero to avoid reference to out-of-bound index in xArray[][]
|
|
||||||
|
|
||||||
if($countColumn == $n) { $countColumn = 0; }
|
|
||||||
|
|
||||||
} // end for
|
|
||||||
|
|
||||||
++$countRow;
|
|
||||||
|
|
||||||
// reset the row-index to zero to avoid reference to out-of-bound index in xArray[][]
|
|
||||||
|
|
||||||
if($countRow == $m) { $countRow = 0; }
|
|
||||||
|
|
||||||
$result[$i] = $holder;
|
|
||||||
|
|
||||||
} // end for
|
|
||||||
|
|
||||||
return new Matrix($result);
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
$X =array(1,2,3,4,5,6,7,8,9);
|
|
||||||
$nRow = 3;
|
|
||||||
$nCol = 3;
|
|
||||||
$tiled_matrix = tile(new Matrix($X), $nRow, $nCol);
|
|
||||||
echo "<pre>";
|
|
||||||
print_r($tiled_matrix);
|
|
||||||
echo "</pre>";
|
|
||||||
?>
|
|
||||||
Loading…
Reference in New Issue