Clean
This commit is contained in:
@ -1,16 +0,0 @@
|
||||
Mar 1, 2005 11:15 AST by PM
|
||||
|
||||
+ For consistency, renamed Math.php to Maths.java, utils to util,
|
||||
tests to test, docs to doc -
|
||||
|
||||
+ Removed conditional logic from top of Matrix class.
|
||||
|
||||
+ Switched to using hypo function in Maths.php for all php-hypot calls.
|
||||
NOTE TO SELF: Need to make sure that all decompositions have been
|
||||
switched over to using the bundled hypo.
|
||||
|
||||
Feb 25, 2005 at 10:00 AST by PM
|
||||
|
||||
+ Recommend using simpler Error.php instead of JAMA_Error.php but
|
||||
can be persuaded otherwise.
|
||||
|
||||
@ -1,149 +0,0 @@
|
||||
<?php
|
||||
/**
|
||||
* @package JAMA
|
||||
*
|
||||
* Cholesky decomposition class
|
||||
*
|
||||
* For a symmetric, positive definite matrix A, the Cholesky decomposition
|
||||
* is an lower triangular matrix L so that A = L*L'.
|
||||
*
|
||||
* If the matrix is not symmetric or positive definite, the constructor
|
||||
* returns a partial decomposition and sets an internal flag that may
|
||||
* be queried by the isSPD() method.
|
||||
*
|
||||
* @author Paul Meagher
|
||||
* @author Michael Bommarito
|
||||
* @version 1.2
|
||||
*/
|
||||
class CholeskyDecomposition {
|
||||
|
||||
/**
|
||||
* Decomposition storage
|
||||
* @var array
|
||||
* @access private
|
||||
*/
|
||||
private $L = array();
|
||||
|
||||
/**
|
||||
* Matrix row and column dimension
|
||||
* @var int
|
||||
* @access private
|
||||
*/
|
||||
private $m;
|
||||
|
||||
/**
|
||||
* Symmetric positive definite flag
|
||||
* @var boolean
|
||||
* @access private
|
||||
*/
|
||||
private $isspd = true;
|
||||
|
||||
|
||||
/**
|
||||
* CholeskyDecomposition
|
||||
*
|
||||
* Class constructor - decomposes symmetric positive definite matrix
|
||||
* @param mixed Matrix square symmetric positive definite matrix
|
||||
*/
|
||||
public function __construct($A = null) {
|
||||
if ($A instanceof Matrix) {
|
||||
$this->L = $A->getArray();
|
||||
$this->m = $A->getRowDimension();
|
||||
|
||||
for($i = 0; $i < $this->m; ++$i) {
|
||||
for($j = $i; $j < $this->m; ++$j) {
|
||||
for($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
|
||||
$sum -= $this->L[$i][$k] * $this->L[$j][$k];
|
||||
}
|
||||
if ($i == $j) {
|
||||
if ($sum >= 0) {
|
||||
$this->L[$i][$i] = sqrt($sum);
|
||||
} else {
|
||||
$this->isspd = false;
|
||||
}
|
||||
} else {
|
||||
if ($this->L[$i][$i] != 0) {
|
||||
$this->L[$j][$i] = $sum / $this->L[$i][$i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for ($k = $i+1; $k < $this->m; ++$k) {
|
||||
$this->L[$i][$k] = 0.0;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
throw new Exception(JAMAError(ArgumentTypeException));
|
||||
}
|
||||
} // function __construct()
|
||||
|
||||
|
||||
/**
|
||||
* Is the matrix symmetric and positive definite?
|
||||
*
|
||||
* @return boolean
|
||||
*/
|
||||
public function isSPD() {
|
||||
return $this->isspd;
|
||||
} // function isSPD()
|
||||
|
||||
|
||||
/**
|
||||
* getL
|
||||
*
|
||||
* Return triangular factor.
|
||||
* @return Matrix Lower triangular matrix
|
||||
*/
|
||||
public function getL() {
|
||||
return new Matrix($this->L);
|
||||
} // function getL()
|
||||
|
||||
|
||||
/**
|
||||
* Solve A*X = B
|
||||
*
|
||||
* @param $B Row-equal matrix
|
||||
* @return Matrix L * L' * X = B
|
||||
*/
|
||||
public function solve($B = null) {
|
||||
if ($B instanceof Matrix) {
|
||||
if ($B->getRowDimension() == $this->m) {
|
||||
if ($this->isspd) {
|
||||
$X = $B->getArrayCopy();
|
||||
$nx = $B->getColumnDimension();
|
||||
|
||||
for ($k = 0; $k < $this->m; ++$k) {
|
||||
for ($i = $k + 1; $i < $this->m; ++$i) {
|
||||
for ($j = 0; $j < $nx; ++$j) {
|
||||
$X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
|
||||
}
|
||||
}
|
||||
for ($j = 0; $j < $nx; ++$j) {
|
||||
$X[$k][$j] /= $this->L[$k][$k];
|
||||
}
|
||||
}
|
||||
|
||||
for ($k = $this->m - 1; $k >= 0; --$k) {
|
||||
for ($j = 0; $j < $nx; ++$j) {
|
||||
$X[$k][$j] /= $this->L[$k][$k];
|
||||
}
|
||||
for ($i = 0; $i < $k; ++$i) {
|
||||
for ($j = 0; $j < $nx; ++$j) {
|
||||
$X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return new Matrix($X, $this->m, $nx);
|
||||
} else {
|
||||
throw new Exception(JAMAError(MatrixSPDException));
|
||||
}
|
||||
} else {
|
||||
throw new Exception(JAMAError(MatrixDimensionException));
|
||||
}
|
||||
} else {
|
||||
throw new Exception(JAMAError(ArgumentTypeException));
|
||||
}
|
||||
} // function solve()
|
||||
|
||||
} // class CholeskyDecomposition
|
||||
@ -1,862 +0,0 @@
|
||||
<?php
|
||||
/**
|
||||
* @package JAMA
|
||||
*
|
||||
* Class to obtain eigenvalues and eigenvectors of a real matrix.
|
||||
*
|
||||
* If A is symmetric, then A = V*D*V' where the eigenvalue matrix D
|
||||
* is diagonal and the eigenvector matrix V is orthogonal (i.e.
|
||||
* A = V.times(D.times(V.transpose())) and V.times(V.transpose())
|
||||
* equals the identity matrix).
|
||||
*
|
||||
* If A is not symmetric, then the eigenvalue matrix D is block diagonal
|
||||
* with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
|
||||
* lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
|
||||
* columns of V represent the eigenvectors in the sense that A*V = V*D,
|
||||
* i.e. A.times(V) equals V.times(D). The matrix V may be badly
|
||||
* conditioned, or even singular, so the validity of the equation
|
||||
* A = V*D*inverse(V) depends upon V.cond().
|
||||
*
|
||||
* @author Paul Meagher
|
||||
* @license PHP v3.0
|
||||
* @version 1.1
|
||||
*/
|
||||
class EigenvalueDecomposition {
|
||||
|
||||
/**
|
||||
* Row and column dimension (square matrix).
|
||||
* @var int
|
||||
*/
|
||||
private $n;
|
||||
|
||||
/**
|
||||
* Internal symmetry flag.
|
||||
* @var int
|
||||
*/
|
||||
private $issymmetric;
|
||||
|
||||
/**
|
||||
* Arrays for internal storage of eigenvalues.
|
||||
* @var array
|
||||
*/
|
||||
private $d = array();
|
||||
private $e = array();
|
||||
|
||||
/**
|
||||
* Array for internal storage of eigenvectors.
|
||||
* @var array
|
||||
*/
|
||||
private $V = array();
|
||||
|
||||
/**
|
||||
* Array for internal storage of nonsymmetric Hessenberg form.
|
||||
* @var array
|
||||
*/
|
||||
private $H = array();
|
||||
|
||||
/**
|
||||
* Working storage for nonsymmetric algorithm.
|
||||
* @var array
|
||||
*/
|
||||
private $ort;
|
||||
|
||||
/**
|
||||
* Used for complex scalar division.
|
||||
* @var float
|
||||
*/
|
||||
private $cdivr;
|
||||
private $cdivi;
|
||||
|
||||
|
||||
/**
|
||||
* Symmetric Householder reduction to tridiagonal form.
|
||||
*
|
||||
* @access private
|
||||
*/
|
||||
private function tred2 () {
|
||||
// This is derived from the Algol procedures tred2 by
|
||||
// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
|
||||
// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
|
||||
// Fortran subroutine in EISPACK.
|
||||
$this->d = $this->V[$this->n-1];
|
||||
// Householder reduction to tridiagonal form.
|
||||
for ($i = $this->n-1; $i > 0; --$i) {
|
||||
$i_ = $i -1;
|
||||
// Scale to avoid under/overflow.
|
||||
$h = $scale = 0.0;
|
||||
$scale += array_sum(array_map(abs, $this->d));
|
||||
if ($scale == 0.0) {
|
||||
$this->e[$i] = $this->d[$i_];
|
||||
$this->d = array_slice($this->V[$i_], 0, $i_);
|
||||
for ($j = 0; $j < $i; ++$j) {
|
||||
$this->V[$j][$i] = $this->V[$i][$j] = 0.0;
|
||||
}
|
||||
} else {
|
||||
// Generate Householder vector.
|
||||
for ($k = 0; $k < $i; ++$k) {
|
||||
$this->d[$k] /= $scale;
|
||||
$h += pow($this->d[$k], 2);
|
||||
}
|
||||
$f = $this->d[$i_];
|
||||
$g = sqrt($h);
|
||||
if ($f > 0) {
|
||||
$g = -$g;
|
||||
}
|
||||
$this->e[$i] = $scale * $g;
|
||||
$h = $h - $f * $g;
|
||||
$this->d[$i_] = $f - $g;
|
||||
for ($j = 0; $j < $i; ++$j) {
|
||||
$this->e[$j] = 0.0;
|
||||
}
|
||||
// Apply similarity transformation to remaining columns.
|
||||
for ($j = 0; $j < $i; ++$j) {
|
||||
$f = $this->d[$j];
|
||||
$this->V[$j][$i] = $f;
|
||||
$g = $this->e[$j] + $this->V[$j][$j] * $f;
|
||||
for ($k = $j+1; $k <= $i_; ++$k) {
|
||||
$g += $this->V[$k][$j] * $this->d[$k];
|
||||
$this->e[$k] += $this->V[$k][$j] * $f;
|
||||
}
|
||||
$this->e[$j] = $g;
|
||||
}
|
||||
$f = 0.0;
|
||||
for ($j = 0; $j < $i; ++$j) {
|
||||
$this->e[$j] /= $h;
|
||||
$f += $this->e[$j] * $this->d[$j];
|
||||
}
|
||||
$hh = $f / (2 * $h);
|
||||
for ($j=0; $j < $i; ++$j) {
|
||||
$this->e[$j] -= $hh * $this->d[$j];
|
||||
}
|
||||
for ($j = 0; $j < $i; ++$j) {
|
||||
$f = $this->d[$j];
|
||||
$g = $this->e[$j];
|
||||
for ($k = $j; $k <= $i_; ++$k) {
|
||||
$this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
|
||||
}
|
||||
$this->d[$j] = $this->V[$i-1][$j];
|
||||
$this->V[$i][$j] = 0.0;
|
||||
}
|
||||
}
|
||||
$this->d[$i] = $h;
|
||||
}
|
||||
|
||||
// Accumulate transformations.
|
||||
for ($i = 0; $i < $this->n-1; ++$i) {
|
||||
$this->V[$this->n-1][$i] = $this->V[$i][$i];
|
||||
$this->V[$i][$i] = 1.0;
|
||||
$h = $this->d[$i+1];
|
||||
if ($h != 0.0) {
|
||||
for ($k = 0; $k <= $i; ++$k) {
|
||||
$this->d[$k] = $this->V[$k][$i+1] / $h;
|
||||
}
|
||||
for ($j = 0; $j <= $i; ++$j) {
|
||||
$g = 0.0;
|
||||
for ($k = 0; $k <= $i; ++$k) {
|
||||
$g += $this->V[$k][$i+1] * $this->V[$k][$j];
|
||||
}
|
||||
for ($k = 0; $k <= $i; ++$k) {
|
||||
$this->V[$k][$j] -= $g * $this->d[$k];
|
||||
}
|
||||
}
|
||||
}
|
||||
for ($k = 0; $k <= $i; ++$k) {
|
||||
$this->V[$k][$i+1] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
$this->d = $this->V[$this->n-1];
|
||||
$this->V[$this->n-1] = array_fill(0, $j, 0.0);
|
||||
$this->V[$this->n-1][$this->n-1] = 1.0;
|
||||
$this->e[0] = 0.0;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Symmetric tridiagonal QL algorithm.
|
||||
*
|
||||
* This is derived from the Algol procedures tql2, by
|
||||
* Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
|
||||
* Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
|
||||
* Fortran subroutine in EISPACK.
|
||||
*
|
||||
* @access private
|
||||
*/
|
||||
private function tql2() {
|
||||
for ($i = 1; $i < $this->n; ++$i) {
|
||||
$this->e[$i-1] = $this->e[$i];
|
||||
}
|
||||
$this->e[$this->n-1] = 0.0;
|
||||
$f = 0.0;
|
||||
$tst1 = 0.0;
|
||||
$eps = pow(2.0,-52.0);
|
||||
|
||||
for ($l = 0; $l < $this->n; ++$l) {
|
||||
// Find small subdiagonal element
|
||||
$tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
|
||||
$m = $l;
|
||||
while ($m < $this->n) {
|
||||
if (abs($this->e[$m]) <= $eps * $tst1)
|
||||
break;
|
||||
++$m;
|
||||
}
|
||||
// If m == l, $this->d[l] is an eigenvalue,
|
||||
// otherwise, iterate.
|
||||
if ($m > $l) {
|
||||
$iter = 0;
|
||||
do {
|
||||
// Could check iteration count here.
|
||||
$iter += 1;
|
||||
// Compute implicit shift
|
||||
$g = $this->d[$l];
|
||||
$p = ($this->d[$l+1] - $g) / (2.0 * $this->e[$l]);
|
||||
$r = hypo($p, 1.0);
|
||||
if ($p < 0)
|
||||
$r *= -1;
|
||||
$this->d[$l] = $this->e[$l] / ($p + $r);
|
||||
$this->d[$l+1] = $this->e[$l] * ($p + $r);
|
||||
$dl1 = $this->d[$l+1];
|
||||
$h = $g - $this->d[$l];
|
||||
for ($i = $l + 2; $i < $this->n; ++$i)
|
||||
$this->d[$i] -= $h;
|
||||
$f += $h;
|
||||
// Implicit QL transformation.
|
||||
$p = $this->d[$m];
|
||||
$c = 1.0;
|
||||
$c2 = $c3 = $c;
|
||||
$el1 = $this->e[$l + 1];
|
||||
$s = $s2 = 0.0;
|
||||
for ($i = $m-1; $i >= $l; --$i) {
|
||||
$c3 = $c2;
|
||||
$c2 = $c;
|
||||
$s2 = $s;
|
||||
$g = $c * $this->e[$i];
|
||||
$h = $c * $p;
|
||||
$r = hypo($p, $this->e[$i]);
|
||||
$this->e[$i+1] = $s * $r;
|
||||
$s = $this->e[$i] / $r;
|
||||
$c = $p / $r;
|
||||
$p = $c * $this->d[$i] - $s * $g;
|
||||
$this->d[$i+1] = $h + $s * ($c * $g + $s * $this->d[$i]);
|
||||
// Accumulate transformation.
|
||||
for ($k = 0; $k < $this->n; ++$k) {
|
||||
$h = $this->V[$k][$i+1];
|
||||
$this->V[$k][$i+1] = $s * $this->V[$k][$i] + $c * $h;
|
||||
$this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
|
||||
}
|
||||
}
|
||||
$p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
|
||||
$this->e[$l] = $s * $p;
|
||||
$this->d[$l] = $c * $p;
|
||||
// Check for convergence.
|
||||
} while (abs($this->e[$l]) > $eps * $tst1);
|
||||
}
|
||||
$this->d[$l] = $this->d[$l] + $f;
|
||||
$this->e[$l] = 0.0;
|
||||
}
|
||||
|
||||
// Sort eigenvalues and corresponding vectors.
|
||||
for ($i = 0; $i < $this->n - 1; ++$i) {
|
||||
$k = $i;
|
||||
$p = $this->d[$i];
|
||||
for ($j = $i+1; $j < $this->n; ++$j) {
|
||||
if ($this->d[$j] < $p) {
|
||||
$k = $j;
|
||||
$p = $this->d[$j];
|
||||
}
|
||||
}
|
||||
if ($k != $i) {
|
||||
$this->d[$k] = $this->d[$i];
|
||||
$this->d[$i] = $p;
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
$p = $this->V[$j][$i];
|
||||
$this->V[$j][$i] = $this->V[$j][$k];
|
||||
$this->V[$j][$k] = $p;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Nonsymmetric reduction to Hessenberg form.
|
||||
*
|
||||
* This is derived from the Algol procedures orthes and ortran,
|
||||
* by Martin and Wilkinson, Handbook for Auto. Comp.,
|
||||
* Vol.ii-Linear Algebra, and the corresponding
|
||||
* Fortran subroutines in EISPACK.
|
||||
*
|
||||
* @access private
|
||||
*/
|
||||
private function orthes () {
|
||||
$low = 0;
|
||||
$high = $this->n-1;
|
||||
|
||||
for ($m = $low+1; $m <= $high-1; ++$m) {
|
||||
// Scale column.
|
||||
$scale = 0.0;
|
||||
for ($i = $m; $i <= $high; ++$i) {
|
||||
$scale = $scale + abs($this->H[$i][$m-1]);
|
||||
}
|
||||
if ($scale != 0.0) {
|
||||
// Compute Householder transformation.
|
||||
$h = 0.0;
|
||||
for ($i = $high; $i >= $m; --$i) {
|
||||
$this->ort[$i] = $this->H[$i][$m-1] / $scale;
|
||||
$h += $this->ort[$i] * $this->ort[$i];
|
||||
}
|
||||
$g = sqrt($h);
|
||||
if ($this->ort[$m] > 0) {
|
||||
$g *= -1;
|
||||
}
|
||||
$h -= $this->ort[$m] * $g;
|
||||
$this->ort[$m] -= $g;
|
||||
// Apply Householder similarity transformation
|
||||
// H = (I -u * u' / h) * H * (I -u * u') / h)
|
||||
for ($j = $m; $j < $this->n; ++$j) {
|
||||
$f = 0.0;
|
||||
for ($i = $high; $i >= $m; --$i) {
|
||||
$f += $this->ort[$i] * $this->H[$i][$j];
|
||||
}
|
||||
$f /= $h;
|
||||
for ($i = $m; $i <= $high; ++$i) {
|
||||
$this->H[$i][$j] -= $f * $this->ort[$i];
|
||||
}
|
||||
}
|
||||
for ($i = 0; $i <= $high; ++$i) {
|
||||
$f = 0.0;
|
||||
for ($j = $high; $j >= $m; --$j) {
|
||||
$f += $this->ort[$j] * $this->H[$i][$j];
|
||||
}
|
||||
$f = $f / $h;
|
||||
for ($j = $m; $j <= $high; ++$j) {
|
||||
$this->H[$i][$j] -= $f * $this->ort[$j];
|
||||
}
|
||||
}
|
||||
$this->ort[$m] = $scale * $this->ort[$m];
|
||||
$this->H[$m][$m-1] = $scale * $g;
|
||||
}
|
||||
}
|
||||
|
||||
// Accumulate transformations (Algol's ortran).
|
||||
for ($i = 0; $i < $this->n; ++$i) {
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
$this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
|
||||
}
|
||||
}
|
||||
for ($m = $high-1; $m >= $low+1; --$m) {
|
||||
if ($this->H[$m][$m-1] != 0.0) {
|
||||
for ($i = $m+1; $i <= $high; ++$i) {
|
||||
$this->ort[$i] = $this->H[$i][$m-1];
|
||||
}
|
||||
for ($j = $m; $j <= $high; ++$j) {
|
||||
$g = 0.0;
|
||||
for ($i = $m; $i <= $high; ++$i) {
|
||||
$g += $this->ort[$i] * $this->V[$i][$j];
|
||||
}
|
||||
// Double division avoids possible underflow
|
||||
$g = ($g / $this->ort[$m]) / $this->H[$m][$m-1];
|
||||
for ($i = $m; $i <= $high; ++$i) {
|
||||
$this->V[$i][$j] += $g * $this->ort[$i];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Performs complex division.
|
||||
*
|
||||
* @access private
|
||||
*/
|
||||
private function cdiv($xr, $xi, $yr, $yi) {
|
||||
if (abs($yr) > abs($yi)) {
|
||||
$r = $yi / $yr;
|
||||
$d = $yr + $r * $yi;
|
||||
$this->cdivr = ($xr + $r * $xi) / $d;
|
||||
$this->cdivi = ($xi - $r * $xr) / $d;
|
||||
} else {
|
||||
$r = $yr / $yi;
|
||||
$d = $yi + $r * $yr;
|
||||
$this->cdivr = ($r * $xr + $xi) / $d;
|
||||
$this->cdivi = ($r * $xi - $xr) / $d;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Nonsymmetric reduction from Hessenberg to real Schur form.
|
||||
*
|
||||
* Code is derived from the Algol procedure hqr2,
|
||||
* by Martin and Wilkinson, Handbook for Auto. Comp.,
|
||||
* Vol.ii-Linear Algebra, and the corresponding
|
||||
* Fortran subroutine in EISPACK.
|
||||
*
|
||||
* @access private
|
||||
*/
|
||||
private function hqr2 () {
|
||||
// Initialize
|
||||
$nn = $this->n;
|
||||
$n = $nn - 1;
|
||||
$low = 0;
|
||||
$high = $nn - 1;
|
||||
$eps = pow(2.0, -52.0);
|
||||
$exshift = 0.0;
|
||||
$p = $q = $r = $s = $z = 0;
|
||||
// Store roots isolated by balanc and compute matrix norm
|
||||
$norm = 0.0;
|
||||
|
||||
for ($i = 0; $i < $nn; ++$i) {
|
||||
if (($i < $low) OR ($i > $high)) {
|
||||
$this->d[$i] = $this->H[$i][$i];
|
||||
$this->e[$i] = 0.0;
|
||||
}
|
||||
for ($j = max($i-1, 0); $j < $nn; ++$j) {
|
||||
$norm = $norm + abs($this->H[$i][$j]);
|
||||
}
|
||||
}
|
||||
|
||||
// Outer loop over eigenvalue index
|
||||
$iter = 0;
|
||||
while ($n >= $low) {
|
||||
// Look for single small sub-diagonal element
|
||||
$l = $n;
|
||||
while ($l > $low) {
|
||||
$s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]);
|
||||
if ($s == 0.0) {
|
||||
$s = $norm;
|
||||
}
|
||||
if (abs($this->H[$l][$l-1]) < $eps * $s) {
|
||||
break;
|
||||
}
|
||||
--$l;
|
||||
}
|
||||
// Check for convergence
|
||||
// One root found
|
||||
if ($l == $n) {
|
||||
$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
|
||||
$this->d[$n] = $this->H[$n][$n];
|
||||
$this->e[$n] = 0.0;
|
||||
--$n;
|
||||
$iter = 0;
|
||||
// Two roots found
|
||||
} else if ($l == $n-1) {
|
||||
$w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
|
||||
$p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0;
|
||||
$q = $p * $p + $w;
|
||||
$z = sqrt(abs($q));
|
||||
$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
|
||||
$this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift;
|
||||
$x = $this->H[$n][$n];
|
||||
// Real pair
|
||||
if ($q >= 0) {
|
||||
if ($p >= 0) {
|
||||
$z = $p + $z;
|
||||
} else {
|
||||
$z = $p - $z;
|
||||
}
|
||||
$this->d[$n-1] = $x + $z;
|
||||
$this->d[$n] = $this->d[$n-1];
|
||||
if ($z != 0.0) {
|
||||
$this->d[$n] = $x - $w / $z;
|
||||
}
|
||||
$this->e[$n-1] = 0.0;
|
||||
$this->e[$n] = 0.0;
|
||||
$x = $this->H[$n][$n-1];
|
||||
$s = abs($x) + abs($z);
|
||||
$p = $x / $s;
|
||||
$q = $z / $s;
|
||||
$r = sqrt($p * $p + $q * $q);
|
||||
$p = $p / $r;
|
||||
$q = $q / $r;
|
||||
// Row modification
|
||||
for ($j = $n-1; $j < $nn; ++$j) {
|
||||
$z = $this->H[$n-1][$j];
|
||||
$this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j];
|
||||
$this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
|
||||
}
|
||||
// Column modification
|
||||
for ($i = 0; $i <= n; ++$i) {
|
||||
$z = $this->H[$i][$n-1];
|
||||
$this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n];
|
||||
$this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
|
||||
}
|
||||
// Accumulate transformations
|
||||
for ($i = $low; $i <= $high; ++$i) {
|
||||
$z = $this->V[$i][$n-1];
|
||||
$this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n];
|
||||
$this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
|
||||
}
|
||||
// Complex pair
|
||||
} else {
|
||||
$this->d[$n-1] = $x + $p;
|
||||
$this->d[$n] = $x + $p;
|
||||
$this->e[$n-1] = $z;
|
||||
$this->e[$n] = -$z;
|
||||
}
|
||||
$n = $n - 2;
|
||||
$iter = 0;
|
||||
// No convergence yet
|
||||
} else {
|
||||
// Form shift
|
||||
$x = $this->H[$n][$n];
|
||||
$y = 0.0;
|
||||
$w = 0.0;
|
||||
if ($l < $n) {
|
||||
$y = $this->H[$n-1][$n-1];
|
||||
$w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
|
||||
}
|
||||
// Wilkinson's original ad hoc shift
|
||||
if ($iter == 10) {
|
||||
$exshift += $x;
|
||||
for ($i = $low; $i <= $n; ++$i) {
|
||||
$this->H[$i][$i] -= $x;
|
||||
}
|
||||
$s = abs($this->H[$n][$n-1]) + abs($this->H[$n-1][$n-2]);
|
||||
$x = $y = 0.75 * $s;
|
||||
$w = -0.4375 * $s * $s;
|
||||
}
|
||||
// MATLAB's new ad hoc shift
|
||||
if ($iter == 30) {
|
||||
$s = ($y - $x) / 2.0;
|
||||
$s = $s * $s + $w;
|
||||
if ($s > 0) {
|
||||
$s = sqrt($s);
|
||||
if ($y < $x) {
|
||||
$s = -$s;
|
||||
}
|
||||
$s = $x - $w / (($y - $x) / 2.0 + $s);
|
||||
for ($i = $low; $i <= $n; ++$i) {
|
||||
$this->H[$i][$i] -= $s;
|
||||
}
|
||||
$exshift += $s;
|
||||
$x = $y = $w = 0.964;
|
||||
}
|
||||
}
|
||||
// Could check iteration count here.
|
||||
$iter = $iter + 1;
|
||||
// Look for two consecutive small sub-diagonal elements
|
||||
$m = $n - 2;
|
||||
while ($m >= $l) {
|
||||
$z = $this->H[$m][$m];
|
||||
$r = $x - $z;
|
||||
$s = $y - $z;
|
||||
$p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1];
|
||||
$q = $this->H[$m+1][$m+1] - $z - $r - $s;
|
||||
$r = $this->H[$m+2][$m+1];
|
||||
$s = abs($p) + abs($q) + abs($r);
|
||||
$p = $p / $s;
|
||||
$q = $q / $s;
|
||||
$r = $r / $s;
|
||||
if ($m == $l) {
|
||||
break;
|
||||
}
|
||||
if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) <
|
||||
$eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) {
|
||||
break;
|
||||
}
|
||||
--$m;
|
||||
}
|
||||
for ($i = $m + 2; $i <= $n; ++$i) {
|
||||
$this->H[$i][$i-2] = 0.0;
|
||||
if ($i > $m+2) {
|
||||
$this->H[$i][$i-3] = 0.0;
|
||||
}
|
||||
}
|
||||
// Double QR step involving rows l:n and columns m:n
|
||||
for ($k = $m; $k <= $n-1; ++$k) {
|
||||
$notlast = ($k != $n-1);
|
||||
if ($k != $m) {
|
||||
$p = $this->H[$k][$k-1];
|
||||
$q = $this->H[$k+1][$k-1];
|
||||
$r = ($notlast ? $this->H[$k+2][$k-1] : 0.0);
|
||||
$x = abs($p) + abs($q) + abs($r);
|
||||
if ($x != 0.0) {
|
||||
$p = $p / $x;
|
||||
$q = $q / $x;
|
||||
$r = $r / $x;
|
||||
}
|
||||
}
|
||||
if ($x == 0.0) {
|
||||
break;
|
||||
}
|
||||
$s = sqrt($p * $p + $q * $q + $r * $r);
|
||||
if ($p < 0) {
|
||||
$s = -$s;
|
||||
}
|
||||
if ($s != 0) {
|
||||
if ($k != $m) {
|
||||
$this->H[$k][$k-1] = -$s * $x;
|
||||
} elseif ($l != $m) {
|
||||
$this->H[$k][$k-1] = -$this->H[$k][$k-1];
|
||||
}
|
||||
$p = $p + $s;
|
||||
$x = $p / $s;
|
||||
$y = $q / $s;
|
||||
$z = $r / $s;
|
||||
$q = $q / $p;
|
||||
$r = $r / $p;
|
||||
// Row modification
|
||||
for ($j = $k; $j < $nn; ++$j) {
|
||||
$p = $this->H[$k][$j] + $q * $this->H[$k+1][$j];
|
||||
if ($notlast) {
|
||||
$p = $p + $r * $this->H[$k+2][$j];
|
||||
$this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z;
|
||||
}
|
||||
$this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
|
||||
$this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y;
|
||||
}
|
||||
// Column modification
|
||||
for ($i = 0; $i <= min($n, $k+3); ++$i) {
|
||||
$p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1];
|
||||
if ($notlast) {
|
||||
$p = $p + $z * $this->H[$i][$k+2];
|
||||
$this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r;
|
||||
}
|
||||
$this->H[$i][$k] = $this->H[$i][$k] - $p;
|
||||
$this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q;
|
||||
}
|
||||
// Accumulate transformations
|
||||
for ($i = $low; $i <= $high; ++$i) {
|
||||
$p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1];
|
||||
if ($notlast) {
|
||||
$p = $p + $z * $this->V[$i][$k+2];
|
||||
$this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r;
|
||||
}
|
||||
$this->V[$i][$k] = $this->V[$i][$k] - $p;
|
||||
$this->V[$i][$k+1] = $this->V[$i][$k+1] - $p * $q;
|
||||
}
|
||||
} // ($s != 0)
|
||||
} // k loop
|
||||
} // check convergence
|
||||
} // while ($n >= $low)
|
||||
|
||||
// Backsubstitute to find vectors of upper triangular form
|
||||
if ($norm == 0.0) {
|
||||
return;
|
||||
}
|
||||
|
||||
for ($n = $nn-1; $n >= 0; --$n) {
|
||||
$p = $this->d[$n];
|
||||
$q = $this->e[$n];
|
||||
// Real vector
|
||||
if ($q == 0) {
|
||||
$l = $n;
|
||||
$this->H[$n][$n] = 1.0;
|
||||
for ($i = $n-1; $i >= 0; --$i) {
|
||||
$w = $this->H[$i][$i] - $p;
|
||||
$r = 0.0;
|
||||
for ($j = $l; $j <= $n; ++$j) {
|
||||
$r = $r + $this->H[$i][$j] * $this->H[$j][$n];
|
||||
}
|
||||
if ($this->e[$i] < 0.0) {
|
||||
$z = $w;
|
||||
$s = $r;
|
||||
} else {
|
||||
$l = $i;
|
||||
if ($this->e[$i] == 0.0) {
|
||||
if ($w != 0.0) {
|
||||
$this->H[$i][$n] = -$r / $w;
|
||||
} else {
|
||||
$this->H[$i][$n] = -$r / ($eps * $norm);
|
||||
}
|
||||
// Solve real equations
|
||||
} else {
|
||||
$x = $this->H[$i][$i+1];
|
||||
$y = $this->H[$i+1][$i];
|
||||
$q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
|
||||
$t = ($x * $s - $z * $r) / $q;
|
||||
$this->H[$i][$n] = $t;
|
||||
if (abs($x) > abs($z)) {
|
||||
$this->H[$i+1][$n] = (-$r - $w * $t) / $x;
|
||||
} else {
|
||||
$this->H[$i+1][$n] = (-$s - $y * $t) / $z;
|
||||
}
|
||||
}
|
||||
// Overflow control
|
||||
$t = abs($this->H[$i][$n]);
|
||||
if (($eps * $t) * $t > 1) {
|
||||
for ($j = $i; $j <= $n; ++$j) {
|
||||
$this->H[$j][$n] = $this->H[$j][$n] / $t;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Complex vector
|
||||
} else if ($q < 0) {
|
||||
$l = $n-1;
|
||||
// Last vector component imaginary so matrix is triangular
|
||||
if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) {
|
||||
$this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1];
|
||||
$this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1];
|
||||
} else {
|
||||
$this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q);
|
||||
$this->H[$n-1][$n-1] = $this->cdivr;
|
||||
$this->H[$n-1][$n] = $this->cdivi;
|
||||
}
|
||||
$this->H[$n][$n-1] = 0.0;
|
||||
$this->H[$n][$n] = 1.0;
|
||||
for ($i = $n-2; $i >= 0; --$i) {
|
||||
// double ra,sa,vr,vi;
|
||||
$ra = 0.0;
|
||||
$sa = 0.0;
|
||||
for ($j = $l; $j <= $n; ++$j) {
|
||||
$ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1];
|
||||
$sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
|
||||
}
|
||||
$w = $this->H[$i][$i] - $p;
|
||||
if ($this->e[$i] < 0.0) {
|
||||
$z = $w;
|
||||
$r = $ra;
|
||||
$s = $sa;
|
||||
} else {
|
||||
$l = $i;
|
||||
if ($this->e[$i] == 0) {
|
||||
$this->cdiv(-$ra, -$sa, $w, $q);
|
||||
$this->H[$i][$n-1] = $this->cdivr;
|
||||
$this->H[$i][$n] = $this->cdivi;
|
||||
} else {
|
||||
// Solve complex equations
|
||||
$x = $this->H[$i][$i+1];
|
||||
$y = $this->H[$i+1][$i];
|
||||
$vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
|
||||
$vi = ($this->d[$i] - $p) * 2.0 * $q;
|
||||
if ($vr == 0.0 & $vi == 0.0) {
|
||||
$vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
|
||||
}
|
||||
$this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
|
||||
$this->H[$i][$n-1] = $this->cdivr;
|
||||
$this->H[$i][$n] = $this->cdivi;
|
||||
if (abs($x) > (abs($z) + abs($q))) {
|
||||
$this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x;
|
||||
$this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x;
|
||||
} else {
|
||||
$this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q);
|
||||
$this->H[$i+1][$n-1] = $this->cdivr;
|
||||
$this->H[$i+1][$n] = $this->cdivi;
|
||||
}
|
||||
}
|
||||
// Overflow control
|
||||
$t = max(abs($this->H[$i][$n-1]),abs($this->H[$i][$n]));
|
||||
if (($eps * $t) * $t > 1) {
|
||||
for ($j = $i; $j <= $n; ++$j) {
|
||||
$this->H[$j][$n-1] = $this->H[$j][$n-1] / $t;
|
||||
$this->H[$j][$n] = $this->H[$j][$n] / $t;
|
||||
}
|
||||
}
|
||||
} // end else
|
||||
} // end for
|
||||
} // end else for complex case
|
||||
} // end for
|
||||
|
||||
// Vectors of isolated roots
|
||||
for ($i = 0; $i < $nn; ++$i) {
|
||||
if ($i < $low | $i > $high) {
|
||||
for ($j = $i; $j < $nn; ++$j) {
|
||||
$this->V[$i][$j] = $this->H[$i][$j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Back transformation to get eigenvectors of original matrix
|
||||
for ($j = $nn-1; $j >= $low; --$j) {
|
||||
for ($i = $low; $i <= $high; ++$i) {
|
||||
$z = 0.0;
|
||||
for ($k = $low; $k <= min($j,$high); ++$k) {
|
||||
$z = $z + $this->V[$i][$k] * $this->H[$k][$j];
|
||||
}
|
||||
$this->V[$i][$j] = $z;
|
||||
}
|
||||
}
|
||||
} // end hqr2
|
||||
|
||||
|
||||
/**
|
||||
* Constructor: Check for symmetry, then construct the eigenvalue decomposition
|
||||
*
|
||||
* @access public
|
||||
* @param A Square matrix
|
||||
* @return Structure to access D and V.
|
||||
*/
|
||||
public function __construct($Arg) {
|
||||
$this->A = $Arg->getArray();
|
||||
$this->n = $Arg->getColumnDimension();
|
||||
|
||||
$issymmetric = true;
|
||||
for ($j = 0; ($j < $this->n) & $issymmetric; ++$j) {
|
||||
for ($i = 0; ($i < $this->n) & $issymmetric; ++$i) {
|
||||
$issymmetric = ($this->A[$i][$j] == $this->A[$j][$i]);
|
||||
}
|
||||
}
|
||||
|
||||
if ($issymmetric) {
|
||||
$this->V = $this->A;
|
||||
// Tridiagonalize.
|
||||
$this->tred2();
|
||||
// Diagonalize.
|
||||
$this->tql2();
|
||||
} else {
|
||||
$this->H = $this->A;
|
||||
$this->ort = array();
|
||||
// Reduce to Hessenberg form.
|
||||
$this->orthes();
|
||||
// Reduce Hessenberg to real Schur form.
|
||||
$this->hqr2();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Return the eigenvector matrix
|
||||
*
|
||||
* @access public
|
||||
* @return V
|
||||
*/
|
||||
public function getV() {
|
||||
return new Matrix($this->V, $this->n, $this->n);
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Return the real parts of the eigenvalues
|
||||
*
|
||||
* @access public
|
||||
* @return real(diag(D))
|
||||
*/
|
||||
public function getRealEigenvalues() {
|
||||
return $this->d;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Return the imaginary parts of the eigenvalues
|
||||
*
|
||||
* @access public
|
||||
* @return imag(diag(D))
|
||||
*/
|
||||
public function getImagEigenvalues() {
|
||||
return $this->e;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Return the block diagonal eigenvalue matrix
|
||||
*
|
||||
* @access public
|
||||
* @return D
|
||||
*/
|
||||
public function getD() {
|
||||
for ($i = 0; $i < $this->n; ++$i) {
|
||||
$D[$i] = array_fill(0, $this->n, 0.0);
|
||||
$D[$i][$i] = $this->d[$i];
|
||||
if ($this->e[$i] == 0) {
|
||||
continue;
|
||||
}
|
||||
$o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
|
||||
$D[$i][$o] = $this->e[$i];
|
||||
}
|
||||
return new Matrix($D);
|
||||
}
|
||||
|
||||
} // class EigenvalueDecomposition
|
||||
@ -1,258 +0,0 @@
|
||||
<?php
|
||||
/**
|
||||
* @package JAMA
|
||||
*
|
||||
* For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
|
||||
* unit lower triangular matrix L, an n-by-n upper triangular matrix U,
|
||||
* and a permutation vector piv of length m so that A(piv,:) = L*U.
|
||||
* If m < n, then L is m-by-m and U is m-by-n.
|
||||
*
|
||||
* The LU decompostion with pivoting always exists, even if the matrix is
|
||||
* singular, so the constructor will never fail. The primary use of the
|
||||
* LU decomposition is in the solution of square systems of simultaneous
|
||||
* linear equations. This will fail if isNonsingular() returns false.
|
||||
*
|
||||
* @author Paul Meagher
|
||||
* @author Bartosz Matosiuk
|
||||
* @author Michael Bommarito
|
||||
* @version 1.1
|
||||
* @license PHP v3.0
|
||||
*/
|
||||
class PHPExcel_Shared_JAMA_LUDecomposition {
|
||||
|
||||
const MatrixSingularException = "Can only perform operation on singular matrix.";
|
||||
const MatrixSquareException = "Mismatched Row dimension";
|
||||
|
||||
/**
|
||||
* Decomposition storage
|
||||
* @var array
|
||||
*/
|
||||
private $LU = array();
|
||||
|
||||
/**
|
||||
* Row dimension.
|
||||
* @var int
|
||||
*/
|
||||
private $m;
|
||||
|
||||
/**
|
||||
* Column dimension.
|
||||
* @var int
|
||||
*/
|
||||
private $n;
|
||||
|
||||
/**
|
||||
* Pivot sign.
|
||||
* @var int
|
||||
*/
|
||||
private $pivsign;
|
||||
|
||||
/**
|
||||
* Internal storage of pivot vector.
|
||||
* @var array
|
||||
*/
|
||||
private $piv = array();
|
||||
|
||||
|
||||
/**
|
||||
* LU Decomposition constructor.
|
||||
*
|
||||
* @param $A Rectangular matrix
|
||||
* @return Structure to access L, U and piv.
|
||||
*/
|
||||
public function __construct($A) {
|
||||
if ($A instanceof PHPExcel_Shared_JAMA_Matrix) {
|
||||
// Use a "left-looking", dot-product, Crout/Doolittle algorithm.
|
||||
$this->LU = $A->getArray();
|
||||
$this->m = $A->getRowDimension();
|
||||
$this->n = $A->getColumnDimension();
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
$this->piv[$i] = $i;
|
||||
}
|
||||
$this->pivsign = 1;
|
||||
$LUrowi = $LUcolj = array();
|
||||
|
||||
// Outer loop.
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
// Make a copy of the j-th column to localize references.
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
$LUcolj[$i] = &$this->LU[$i][$j];
|
||||
}
|
||||
// Apply previous transformations.
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
$LUrowi = $this->LU[$i];
|
||||
// Most of the time is spent in the following dot product.
|
||||
$kmax = min($i,$j);
|
||||
$s = 0.0;
|
||||
for ($k = 0; $k < $kmax; ++$k) {
|
||||
$s += $LUrowi[$k] * $LUcolj[$k];
|
||||
}
|
||||
$LUrowi[$j] = $LUcolj[$i] -= $s;
|
||||
}
|
||||
// Find pivot and exchange if necessary.
|
||||
$p = $j;
|
||||
for ($i = $j+1; $i < $this->m; ++$i) {
|
||||
if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
|
||||
$p = $i;
|
||||
}
|
||||
}
|
||||
if ($p != $j) {
|
||||
for ($k = 0; $k < $this->n; ++$k) {
|
||||
$t = $this->LU[$p][$k];
|
||||
$this->LU[$p][$k] = $this->LU[$j][$k];
|
||||
$this->LU[$j][$k] = $t;
|
||||
}
|
||||
$k = $this->piv[$p];
|
||||
$this->piv[$p] = $this->piv[$j];
|
||||
$this->piv[$j] = $k;
|
||||
$this->pivsign = $this->pivsign * -1;
|
||||
}
|
||||
// Compute multipliers.
|
||||
if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
|
||||
for ($i = $j+1; $i < $this->m; ++$i) {
|
||||
$this->LU[$i][$j] /= $this->LU[$j][$j];
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
throw new Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
|
||||
}
|
||||
} // function __construct()
|
||||
|
||||
|
||||
/**
|
||||
* Get lower triangular factor.
|
||||
*
|
||||
* @return array Lower triangular factor
|
||||
*/
|
||||
public function getL() {
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
if ($i > $j) {
|
||||
$L[$i][$j] = $this->LU[$i][$j];
|
||||
} elseif ($i == $j) {
|
||||
$L[$i][$j] = 1.0;
|
||||
} else {
|
||||
$L[$i][$j] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
return new PHPExcel_Shared_JAMA_Matrix($L);
|
||||
} // function getL()
|
||||
|
||||
|
||||
/**
|
||||
* Get upper triangular factor.
|
||||
*
|
||||
* @return array Upper triangular factor
|
||||
*/
|
||||
public function getU() {
|
||||
for ($i = 0; $i < $this->n; ++$i) {
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
if ($i <= $j) {
|
||||
$U[$i][$j] = $this->LU[$i][$j];
|
||||
} else {
|
||||
$U[$i][$j] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
return new PHPExcel_Shared_JAMA_Matrix($U);
|
||||
} // function getU()
|
||||
|
||||
|
||||
/**
|
||||
* Return pivot permutation vector.
|
||||
*
|
||||
* @return array Pivot vector
|
||||
*/
|
||||
public function getPivot() {
|
||||
return $this->piv;
|
||||
} // function getPivot()
|
||||
|
||||
|
||||
/**
|
||||
* Alias for getPivot
|
||||
*
|
||||
* @see getPivot
|
||||
*/
|
||||
public function getDoublePivot() {
|
||||
return $this->getPivot();
|
||||
} // function getDoublePivot()
|
||||
|
||||
|
||||
/**
|
||||
* Is the matrix nonsingular?
|
||||
*
|
||||
* @return true if U, and hence A, is nonsingular.
|
||||
*/
|
||||
public function isNonsingular() {
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
if ($this->LU[$j][$j] == 0) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
} // function isNonsingular()
|
||||
|
||||
|
||||
/**
|
||||
* Count determinants
|
||||
*
|
||||
* @return array d matrix deterninat
|
||||
*/
|
||||
public function det() {
|
||||
if ($this->m == $this->n) {
|
||||
$d = $this->pivsign;
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
$d *= $this->LU[$j][$j];
|
||||
}
|
||||
return $d;
|
||||
} else {
|
||||
throw new Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
|
||||
}
|
||||
} // function det()
|
||||
|
||||
|
||||
/**
|
||||
* Solve A*X = B
|
||||
*
|
||||
* @param $B A Matrix with as many rows as A and any number of columns.
|
||||
* @return X so that L*U*X = B(piv,:)
|
||||
* @exception IllegalArgumentException Matrix row dimensions must agree.
|
||||
* @exception RuntimeException Matrix is singular.
|
||||
*/
|
||||
public function solve($B) {
|
||||
if ($B->getRowDimension() == $this->m) {
|
||||
if ($this->isNonsingular()) {
|
||||
// Copy right hand side with pivoting
|
||||
$nx = $B->getColumnDimension();
|
||||
$X = $B->getMatrix($this->piv, 0, $nx-1);
|
||||
// Solve L*Y = B(piv,:)
|
||||
for ($k = 0; $k < $this->n; ++$k) {
|
||||
for ($i = $k+1; $i < $this->n; ++$i) {
|
||||
for ($j = 0; $j < $nx; ++$j) {
|
||||
$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
|
||||
}
|
||||
}
|
||||
}
|
||||
// Solve U*X = Y;
|
||||
for ($k = $this->n-1; $k >= 0; --$k) {
|
||||
for ($j = 0; $j < $nx; ++$j) {
|
||||
$X->A[$k][$j] /= $this->LU[$k][$k];
|
||||
}
|
||||
for ($i = 0; $i < $k; ++$i) {
|
||||
for ($j = 0; $j < $nx; ++$j) {
|
||||
$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
|
||||
}
|
||||
}
|
||||
}
|
||||
return $X;
|
||||
} else {
|
||||
throw new Exception(self::MatrixSingularException);
|
||||
}
|
||||
} else {
|
||||
throw new Exception(self::MatrixSquareException);
|
||||
}
|
||||
} // function solve()
|
||||
|
||||
} // class PHPExcel_Shared_JAMA_LUDecomposition
|
||||
File diff suppressed because it is too large
Load Diff
@ -1,234 +0,0 @@
|
||||
<?php
|
||||
/**
|
||||
* @package JAMA
|
||||
*
|
||||
* For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
|
||||
* orthogonal matrix Q and an n-by-n upper triangular matrix R so that
|
||||
* A = Q*R.
|
||||
*
|
||||
* The QR decompostion always exists, even if the matrix does not have
|
||||
* full rank, so the constructor will never fail. The primary use of the
|
||||
* QR decomposition is in the least squares solution of nonsquare systems
|
||||
* of simultaneous linear equations. This will fail if isFullRank()
|
||||
* returns false.
|
||||
*
|
||||
* @author Paul Meagher
|
||||
* @license PHP v3.0
|
||||
* @version 1.1
|
||||
*/
|
||||
class PHPExcel_Shared_JAMA_QRDecomposition {
|
||||
|
||||
const MatrixRankException = "Can only perform operation on full-rank matrix.";
|
||||
|
||||
/**
|
||||
* Array for internal storage of decomposition.
|
||||
* @var array
|
||||
*/
|
||||
private $QR = array();
|
||||
|
||||
/**
|
||||
* Row dimension.
|
||||
* @var integer
|
||||
*/
|
||||
private $m;
|
||||
|
||||
/**
|
||||
* Column dimension.
|
||||
* @var integer
|
||||
*/
|
||||
private $n;
|
||||
|
||||
/**
|
||||
* Array for internal storage of diagonal of R.
|
||||
* @var array
|
||||
*/
|
||||
private $Rdiag = array();
|
||||
|
||||
|
||||
/**
|
||||
* QR Decomposition computed by Householder reflections.
|
||||
*
|
||||
* @param matrix $A Rectangular matrix
|
||||
* @return Structure to access R and the Householder vectors and compute Q.
|
||||
*/
|
||||
public function __construct($A) {
|
||||
if($A instanceof PHPExcel_Shared_JAMA_Matrix) {
|
||||
// Initialize.
|
||||
$this->QR = $A->getArrayCopy();
|
||||
$this->m = $A->getRowDimension();
|
||||
$this->n = $A->getColumnDimension();
|
||||
// Main loop.
|
||||
for ($k = 0; $k < $this->n; ++$k) {
|
||||
// Compute 2-norm of k-th column without under/overflow.
|
||||
$nrm = 0.0;
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$nrm = hypo($nrm, $this->QR[$i][$k]);
|
||||
}
|
||||
if ($nrm != 0.0) {
|
||||
// Form k-th Householder vector.
|
||||
if ($this->QR[$k][$k] < 0) {
|
||||
$nrm = -$nrm;
|
||||
}
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$this->QR[$i][$k] /= $nrm;
|
||||
}
|
||||
$this->QR[$k][$k] += 1.0;
|
||||
// Apply transformation to remaining columns.
|
||||
for ($j = $k+1; $j < $this->n; ++$j) {
|
||||
$s = 0.0;
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$s += $this->QR[$i][$k] * $this->QR[$i][$j];
|
||||
}
|
||||
$s = -$s/$this->QR[$k][$k];
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$this->QR[$i][$j] += $s * $this->QR[$i][$k];
|
||||
}
|
||||
}
|
||||
}
|
||||
$this->Rdiag[$k] = -$nrm;
|
||||
}
|
||||
} else {
|
||||
throw new Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
|
||||
}
|
||||
} // function __construct()
|
||||
|
||||
|
||||
/**
|
||||
* Is the matrix full rank?
|
||||
*
|
||||
* @return boolean true if R, and hence A, has full rank, else false.
|
||||
*/
|
||||
public function isFullRank() {
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
if ($this->Rdiag[$j] == 0) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
} // function isFullRank()
|
||||
|
||||
|
||||
/**
|
||||
* Return the Householder vectors
|
||||
*
|
||||
* @return Matrix Lower trapezoidal matrix whose columns define the reflections
|
||||
*/
|
||||
public function getH() {
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
if ($i >= $j) {
|
||||
$H[$i][$j] = $this->QR[$i][$j];
|
||||
} else {
|
||||
$H[$i][$j] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
return new PHPExcel_Shared_JAMA_Matrix($H);
|
||||
} // function getH()
|
||||
|
||||
|
||||
/**
|
||||
* Return the upper triangular factor
|
||||
*
|
||||
* @return Matrix upper triangular factor
|
||||
*/
|
||||
public function getR() {
|
||||
for ($i = 0; $i < $this->n; ++$i) {
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
if ($i < $j) {
|
||||
$R[$i][$j] = $this->QR[$i][$j];
|
||||
} elseif ($i == $j) {
|
||||
$R[$i][$j] = $this->Rdiag[$i];
|
||||
} else {
|
||||
$R[$i][$j] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
return new PHPExcel_Shared_JAMA_Matrix($R);
|
||||
} // function getR()
|
||||
|
||||
|
||||
/**
|
||||
* Generate and return the (economy-sized) orthogonal factor
|
||||
*
|
||||
* @return Matrix orthogonal factor
|
||||
*/
|
||||
public function getQ() {
|
||||
for ($k = $this->n-1; $k >= 0; --$k) {
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
$Q[$i][$k] = 0.0;
|
||||
}
|
||||
$Q[$k][$k] = 1.0;
|
||||
for ($j = $k; $j < $this->n; ++$j) {
|
||||
if ($this->QR[$k][$k] != 0) {
|
||||
$s = 0.0;
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$s += $this->QR[$i][$k] * $Q[$i][$j];
|
||||
}
|
||||
$s = -$s/$this->QR[$k][$k];
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$Q[$i][$j] += $s * $this->QR[$i][$k];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
/*
|
||||
for($i = 0; $i < count($Q); ++$i) {
|
||||
for($j = 0; $j < count($Q); ++$j) {
|
||||
if(! isset($Q[$i][$j]) ) {
|
||||
$Q[$i][$j] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
*/
|
||||
return new PHPExcel_Shared_JAMA_Matrix($Q);
|
||||
} // function getQ()
|
||||
|
||||
|
||||
/**
|
||||
* Least squares solution of A*X = B
|
||||
*
|
||||
* @param Matrix $B A Matrix with as many rows as A and any number of columns.
|
||||
* @return Matrix Matrix that minimizes the two norm of Q*R*X-B.
|
||||
*/
|
||||
public function solve($B) {
|
||||
if ($B->getRowDimension() == $this->m) {
|
||||
if ($this->isFullRank()) {
|
||||
// Copy right hand side
|
||||
$nx = $B->getColumnDimension();
|
||||
$X = $B->getArrayCopy();
|
||||
// Compute Y = transpose(Q)*B
|
||||
for ($k = 0; $k < $this->n; ++$k) {
|
||||
for ($j = 0; $j < $nx; ++$j) {
|
||||
$s = 0.0;
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$s += $this->QR[$i][$k] * $X[$i][$j];
|
||||
}
|
||||
$s = -$s/$this->QR[$k][$k];
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$X[$i][$j] += $s * $this->QR[$i][$k];
|
||||
}
|
||||
}
|
||||
}
|
||||
// Solve R*X = Y;
|
||||
for ($k = $this->n-1; $k >= 0; --$k) {
|
||||
for ($j = 0; $j < $nx; ++$j) {
|
||||
$X[$k][$j] /= $this->Rdiag[$k];
|
||||
}
|
||||
for ($i = 0; $i < $k; ++$i) {
|
||||
for ($j = 0; $j < $nx; ++$j) {
|
||||
$X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];
|
||||
}
|
||||
}
|
||||
}
|
||||
$X = new PHPExcel_Shared_JAMA_Matrix($X);
|
||||
return ($X->getMatrix(0, $this->n-1, 0, $nx));
|
||||
} else {
|
||||
throw new Exception(self::MatrixRankException);
|
||||
}
|
||||
} else {
|
||||
throw new Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
|
||||
}
|
||||
} // function solve()
|
||||
|
||||
} // PHPExcel_Shared_JAMA_class QRDecomposition
|
||||
@ -1,526 +0,0 @@
|
||||
<?php
|
||||
/**
|
||||
* @package JAMA
|
||||
*
|
||||
* For an m-by-n matrix A with m >= n, the singular value decomposition is
|
||||
* an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
|
||||
* an n-by-n orthogonal matrix V so that A = U*S*V'.
|
||||
*
|
||||
* The singular values, sigma[$k] = S[$k][$k], are ordered so that
|
||||
* sigma[0] >= sigma[1] >= ... >= sigma[n-1].
|
||||
*
|
||||
* The singular value decompostion always exists, so the constructor will
|
||||
* never fail. The matrix condition number and the effective numerical
|
||||
* rank can be computed from this decomposition.
|
||||
*
|
||||
* @author Paul Meagher
|
||||
* @license PHP v3.0
|
||||
* @version 1.1
|
||||
*/
|
||||
class SingularValueDecomposition {
|
||||
|
||||
/**
|
||||
* Internal storage of U.
|
||||
* @var array
|
||||
*/
|
||||
private $U = array();
|
||||
|
||||
/**
|
||||
* Internal storage of V.
|
||||
* @var array
|
||||
*/
|
||||
private $V = array();
|
||||
|
||||
/**
|
||||
* Internal storage of singular values.
|
||||
* @var array
|
||||
*/
|
||||
private $s = array();
|
||||
|
||||
/**
|
||||
* Row dimension.
|
||||
* @var int
|
||||
*/
|
||||
private $m;
|
||||
|
||||
/**
|
||||
* Column dimension.
|
||||
* @var int
|
||||
*/
|
||||
private $n;
|
||||
|
||||
|
||||
/**
|
||||
* Construct the singular value decomposition
|
||||
*
|
||||
* Derived from LINPACK code.
|
||||
*
|
||||
* @param $A Rectangular matrix
|
||||
* @return Structure to access U, S and V.
|
||||
*/
|
||||
public function __construct($Arg) {
|
||||
|
||||
// Initialize.
|
||||
$A = $Arg->getArrayCopy();
|
||||
$this->m = $Arg->getRowDimension();
|
||||
$this->n = $Arg->getColumnDimension();
|
||||
$nu = min($this->m, $this->n);
|
||||
$e = array();
|
||||
$work = array();
|
||||
$wantu = true;
|
||||
$wantv = true;
|
||||
$nct = min($this->m - 1, $this->n);
|
||||
$nrt = max(0, min($this->n - 2, $this->m));
|
||||
|
||||
// Reduce A to bidiagonal form, storing the diagonal elements
|
||||
// in s and the super-diagonal elements in e.
|
||||
for ($k = 0; $k < max($nct,$nrt); ++$k) {
|
||||
|
||||
if ($k < $nct) {
|
||||
// Compute the transformation for the k-th column and
|
||||
// place the k-th diagonal in s[$k].
|
||||
// Compute 2-norm of k-th column without under/overflow.
|
||||
$this->s[$k] = 0;
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
|
||||
}
|
||||
if ($this->s[$k] != 0.0) {
|
||||
if ($A[$k][$k] < 0.0) {
|
||||
$this->s[$k] = -$this->s[$k];
|
||||
}
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$A[$i][$k] /= $this->s[$k];
|
||||
}
|
||||
$A[$k][$k] += 1.0;
|
||||
}
|
||||
$this->s[$k] = -$this->s[$k];
|
||||
}
|
||||
|
||||
for ($j = $k + 1; $j < $this->n; ++$j) {
|
||||
if (($k < $nct) & ($this->s[$k] != 0.0)) {
|
||||
// Apply the transformation.
|
||||
$t = 0;
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$t += $A[$i][$k] * $A[$i][$j];
|
||||
}
|
||||
$t = -$t / $A[$k][$k];
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$A[$i][$j] += $t * $A[$i][$k];
|
||||
}
|
||||
// Place the k-th row of A into e for the
|
||||
// subsequent calculation of the row transformation.
|
||||
$e[$j] = $A[$k][$j];
|
||||
}
|
||||
}
|
||||
|
||||
if ($wantu AND ($k < $nct)) {
|
||||
// Place the transformation in U for subsequent back
|
||||
// multiplication.
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$this->U[$i][$k] = $A[$i][$k];
|
||||
}
|
||||
}
|
||||
|
||||
if ($k < $nrt) {
|
||||
// Compute the k-th row transformation and place the
|
||||
// k-th super-diagonal in e[$k].
|
||||
// Compute 2-norm without under/overflow.
|
||||
$e[$k] = 0;
|
||||
for ($i = $k + 1; $i < $this->n; ++$i) {
|
||||
$e[$k] = hypo($e[$k], $e[$i]);
|
||||
}
|
||||
if ($e[$k] != 0.0) {
|
||||
if ($e[$k+1] < 0.0) {
|
||||
$e[$k] = -$e[$k];
|
||||
}
|
||||
for ($i = $k + 1; $i < $this->n; ++$i) {
|
||||
$e[$i] /= $e[$k];
|
||||
}
|
||||
$e[$k+1] += 1.0;
|
||||
}
|
||||
$e[$k] = -$e[$k];
|
||||
if (($k+1 < $this->m) AND ($e[$k] != 0.0)) {
|
||||
// Apply the transformation.
|
||||
for ($i = $k+1; $i < $this->m; ++$i) {
|
||||
$work[$i] = 0.0;
|
||||
}
|
||||
for ($j = $k+1; $j < $this->n; ++$j) {
|
||||
for ($i = $k+1; $i < $this->m; ++$i) {
|
||||
$work[$i] += $e[$j] * $A[$i][$j];
|
||||
}
|
||||
}
|
||||
for ($j = $k + 1; $j < $this->n; ++$j) {
|
||||
$t = -$e[$j] / $e[$k+1];
|
||||
for ($i = $k + 1; $i < $this->m; ++$i) {
|
||||
$A[$i][$j] += $t * $work[$i];
|
||||
}
|
||||
}
|
||||
}
|
||||
if ($wantv) {
|
||||
// Place the transformation in V for subsequent
|
||||
// back multiplication.
|
||||
for ($i = $k + 1; $i < $this->n; ++$i) {
|
||||
$this->V[$i][$k] = $e[$i];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Set up the final bidiagonal matrix or order p.
|
||||
$p = min($this->n, $this->m + 1);
|
||||
if ($nct < $this->n) {
|
||||
$this->s[$nct] = $A[$nct][$nct];
|
||||
}
|
||||
if ($this->m < $p) {
|
||||
$this->s[$p-1] = 0.0;
|
||||
}
|
||||
if ($nrt + 1 < $p) {
|
||||
$e[$nrt] = $A[$nrt][$p-1];
|
||||
}
|
||||
$e[$p-1] = 0.0;
|
||||
// If required, generate U.
|
||||
if ($wantu) {
|
||||
for ($j = $nct; $j < $nu; ++$j) {
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
$this->U[$i][$j] = 0.0;
|
||||
}
|
||||
$this->U[$j][$j] = 1.0;
|
||||
}
|
||||
for ($k = $nct - 1; $k >= 0; --$k) {
|
||||
if ($this->s[$k] != 0.0) {
|
||||
for ($j = $k + 1; $j < $nu; ++$j) {
|
||||
$t = 0;
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$t += $this->U[$i][$k] * $this->U[$i][$j];
|
||||
}
|
||||
$t = -$t / $this->U[$k][$k];
|
||||
for ($i = $k; $i < $this->m; ++$i) {
|
||||
$this->U[$i][$j] += $t * $this->U[$i][$k];
|
||||
}
|
||||
}
|
||||
for ($i = $k; $i < $this->m; ++$i ) {
|
||||
$this->U[$i][$k] = -$this->U[$i][$k];
|
||||
}
|
||||
$this->U[$k][$k] = 1.0 + $this->U[$k][$k];
|
||||
for ($i = 0; $i < $k - 1; ++$i) {
|
||||
$this->U[$i][$k] = 0.0;
|
||||
}
|
||||
} else {
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
$this->U[$i][$k] = 0.0;
|
||||
}
|
||||
$this->U[$k][$k] = 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// If required, generate V.
|
||||
if ($wantv) {
|
||||
for ($k = $this->n - 1; $k >= 0; --$k) {
|
||||
if (($k < $nrt) AND ($e[$k] != 0.0)) {
|
||||
for ($j = $k + 1; $j < $nu; ++$j) {
|
||||
$t = 0;
|
||||
for ($i = $k + 1; $i < $this->n; ++$i) {
|
||||
$t += $this->V[$i][$k]* $this->V[$i][$j];
|
||||
}
|
||||
$t = -$t / $this->V[$k+1][$k];
|
||||
for ($i = $k + 1; $i < $this->n; ++$i) {
|
||||
$this->V[$i][$j] += $t * $this->V[$i][$k];
|
||||
}
|
||||
}
|
||||
}
|
||||
for ($i = 0; $i < $this->n; ++$i) {
|
||||
$this->V[$i][$k] = 0.0;
|
||||
}
|
||||
$this->V[$k][$k] = 1.0;
|
||||
}
|
||||
}
|
||||
|
||||
// Main iteration loop for the singular values.
|
||||
$pp = $p - 1;
|
||||
$iter = 0;
|
||||
$eps = pow(2.0, -52.0);
|
||||
|
||||
while ($p > 0) {
|
||||
// Here is where a test for too many iterations would go.
|
||||
// This section of the program inspects for negligible
|
||||
// elements in the s and e arrays. On completion the
|
||||
// variables kase and k are set as follows:
|
||||
// kase = 1 if s(p) and e[k-1] are negligible and k<p
|
||||
// kase = 2 if s(k) is negligible and k<p
|
||||
// kase = 3 if e[k-1] is negligible, k<p, and
|
||||
// s(k), ..., s(p) are not negligible (qr step).
|
||||
// kase = 4 if e(p-1) is negligible (convergence).
|
||||
for ($k = $p - 2; $k >= -1; --$k) {
|
||||
if ($k == -1) {
|
||||
break;
|
||||
}
|
||||
if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) {
|
||||
$e[$k] = 0.0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if ($k == $p - 2) {
|
||||
$kase = 4;
|
||||
} else {
|
||||
for ($ks = $p - 1; $ks >= $k; --$ks) {
|
||||
if ($ks == $k) {
|
||||
break;
|
||||
}
|
||||
$t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.);
|
||||
if (abs($this->s[$ks]) <= $eps * $t) {
|
||||
$this->s[$ks] = 0.0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if ($ks == $k) {
|
||||
$kase = 3;
|
||||
} else if ($ks == $p-1) {
|
||||
$kase = 1;
|
||||
} else {
|
||||
$kase = 2;
|
||||
$k = $ks;
|
||||
}
|
||||
}
|
||||
++$k;
|
||||
|
||||
// Perform the task indicated by kase.
|
||||
switch ($kase) {
|
||||
// Deflate negligible s(p).
|
||||
case 1:
|
||||
$f = $e[$p-2];
|
||||
$e[$p-2] = 0.0;
|
||||
for ($j = $p - 2; $j >= $k; --$j) {
|
||||
$t = hypo($this->s[$j],$f);
|
||||
$cs = $this->s[$j] / $t;
|
||||
$sn = $f / $t;
|
||||
$this->s[$j] = $t;
|
||||
if ($j != $k) {
|
||||
$f = -$sn * $e[$j-1];
|
||||
$e[$j-1] = $cs * $e[$j-1];
|
||||
}
|
||||
if ($wantv) {
|
||||
for ($i = 0; $i < $this->n; ++$i) {
|
||||
$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1];
|
||||
$this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1];
|
||||
$this->V[$i][$j] = $t;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
// Split at negligible s(k).
|
||||
case 2:
|
||||
$f = $e[$k-1];
|
||||
$e[$k-1] = 0.0;
|
||||
for ($j = $k; $j < $p; ++$j) {
|
||||
$t = hypo($this->s[$j], $f);
|
||||
$cs = $this->s[$j] / $t;
|
||||
$sn = $f / $t;
|
||||
$this->s[$j] = $t;
|
||||
$f = -$sn * $e[$j];
|
||||
$e[$j] = $cs * $e[$j];
|
||||
if ($wantu) {
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1];
|
||||
$this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1];
|
||||
$this->U[$i][$j] = $t;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
// Perform one qr step.
|
||||
case 3:
|
||||
// Calculate the shift.
|
||||
$scale = max(max(max(max(
|
||||
abs($this->s[$p-1]),abs($this->s[$p-2])),abs($e[$p-2])),
|
||||
abs($this->s[$k])), abs($e[$k]));
|
||||
$sp = $this->s[$p-1] / $scale;
|
||||
$spm1 = $this->s[$p-2] / $scale;
|
||||
$epm1 = $e[$p-2] / $scale;
|
||||
$sk = $this->s[$k] / $scale;
|
||||
$ek = $e[$k] / $scale;
|
||||
$b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
|
||||
$c = ($sp * $epm1) * ($sp * $epm1);
|
||||
$shift = 0.0;
|
||||
if (($b != 0.0) || ($c != 0.0)) {
|
||||
$shift = sqrt($b * $b + $c);
|
||||
if ($b < 0.0) {
|
||||
$shift = -$shift;
|
||||
}
|
||||
$shift = $c / ($b + $shift);
|
||||
}
|
||||
$f = ($sk + $sp) * ($sk - $sp) + $shift;
|
||||
$g = $sk * $ek;
|
||||
// Chase zeros.
|
||||
for ($j = $k; $j < $p-1; ++$j) {
|
||||
$t = hypo($f,$g);
|
||||
$cs = $f/$t;
|
||||
$sn = $g/$t;
|
||||
if ($j != $k) {
|
||||
$e[$j-1] = $t;
|
||||
}
|
||||
$f = $cs * $this->s[$j] + $sn * $e[$j];
|
||||
$e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
|
||||
$g = $sn * $this->s[$j+1];
|
||||
$this->s[$j+1] = $cs * $this->s[$j+1];
|
||||
if ($wantv) {
|
||||
for ($i = 0; $i < $this->n; ++$i) {
|
||||
$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1];
|
||||
$this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1];
|
||||
$this->V[$i][$j] = $t;
|
||||
}
|
||||
}
|
||||
$t = hypo($f,$g);
|
||||
$cs = $f/$t;
|
||||
$sn = $g/$t;
|
||||
$this->s[$j] = $t;
|
||||
$f = $cs * $e[$j] + $sn * $this->s[$j+1];
|
||||
$this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1];
|
||||
$g = $sn * $e[$j+1];
|
||||
$e[$j+1] = $cs * $e[$j+1];
|
||||
if ($wantu && ($j < $this->m - 1)) {
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1];
|
||||
$this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1];
|
||||
$this->U[$i][$j] = $t;
|
||||
}
|
||||
}
|
||||
}
|
||||
$e[$p-2] = $f;
|
||||
$iter = $iter + 1;
|
||||
break;
|
||||
// Convergence.
|
||||
case 4:
|
||||
// Make the singular values positive.
|
||||
if ($this->s[$k] <= 0.0) {
|
||||
$this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);
|
||||
if ($wantv) {
|
||||
for ($i = 0; $i <= $pp; ++$i) {
|
||||
$this->V[$i][$k] = -$this->V[$i][$k];
|
||||
}
|
||||
}
|
||||
}
|
||||
// Order the singular values.
|
||||
while ($k < $pp) {
|
||||
if ($this->s[$k] >= $this->s[$k+1]) {
|
||||
break;
|
||||
}
|
||||
$t = $this->s[$k];
|
||||
$this->s[$k] = $this->s[$k+1];
|
||||
$this->s[$k+1] = $t;
|
||||
if ($wantv AND ($k < $this->n - 1)) {
|
||||
for ($i = 0; $i < $this->n; ++$i) {
|
||||
$t = $this->V[$i][$k+1];
|
||||
$this->V[$i][$k+1] = $this->V[$i][$k];
|
||||
$this->V[$i][$k] = $t;
|
||||
}
|
||||
}
|
||||
if ($wantu AND ($k < $this->m-1)) {
|
||||
for ($i = 0; $i < $this->m; ++$i) {
|
||||
$t = $this->U[$i][$k+1];
|
||||
$this->U[$i][$k+1] = $this->U[$i][$k];
|
||||
$this->U[$i][$k] = $t;
|
||||
}
|
||||
}
|
||||
++$k;
|
||||
}
|
||||
$iter = 0;
|
||||
--$p;
|
||||
break;
|
||||
} // end switch
|
||||
} // end while
|
||||
|
||||
} // end constructor
|
||||
|
||||
|
||||
/**
|
||||
* Return the left singular vectors
|
||||
*
|
||||
* @access public
|
||||
* @return U
|
||||
*/
|
||||
public function getU() {
|
||||
return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Return the right singular vectors
|
||||
*
|
||||
* @access public
|
||||
* @return V
|
||||
*/
|
||||
public function getV() {
|
||||
return new Matrix($this->V);
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Return the one-dimensional array of singular values
|
||||
*
|
||||
* @access public
|
||||
* @return diagonal of S.
|
||||
*/
|
||||
public function getSingularValues() {
|
||||
return $this->s;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Return the diagonal matrix of singular values
|
||||
*
|
||||
* @access public
|
||||
* @return S
|
||||
*/
|
||||
public function getS() {
|
||||
for ($i = 0; $i < $this->n; ++$i) {
|
||||
for ($j = 0; $j < $this->n; ++$j) {
|
||||
$S[$i][$j] = 0.0;
|
||||
}
|
||||
$S[$i][$i] = $this->s[$i];
|
||||
}
|
||||
return new Matrix($S);
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Two norm
|
||||
*
|
||||
* @access public
|
||||
* @return max(S)
|
||||
*/
|
||||
public function norm2() {
|
||||
return $this->s[0];
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Two norm condition number
|
||||
*
|
||||
* @access public
|
||||
* @return max(S)/min(S)
|
||||
*/
|
||||
public function cond() {
|
||||
return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Effective numerical matrix rank
|
||||
*
|
||||
* @access public
|
||||
* @return Number of nonnegligible singular values.
|
||||
*/
|
||||
public function rank() {
|
||||
$eps = pow(2.0, -52.0);
|
||||
$tol = max($this->m, $this->n) * $this->s[0] * $eps;
|
||||
$r = 0;
|
||||
for ($i = 0; $i < count($this->s); ++$i) {
|
||||
if ($this->s[$i] > $tol) {
|
||||
++$r;
|
||||
}
|
||||
}
|
||||
return $r;
|
||||
}
|
||||
|
||||
} // class SingularValueDecomposition
|
||||
@ -1,6 +0,0 @@
|
||||
<?php
|
||||
require_once "includes/header.php";
|
||||
require_once "includes/navbar.php";
|
||||
require_once "sections/Home.php";
|
||||
require_once "includes/footer.php";
|
||||
?>
|
||||
@ -1,65 +0,0 @@
|
||||
<?php
|
||||
/**
|
||||
* Script to create REGRESS package for download
|
||||
*
|
||||
* @author Mike Bommarito
|
||||
* @author Paul Meagher
|
||||
* @version 0.3
|
||||
* @modified Apr 2, 2006
|
||||
*
|
||||
* Note: Script requires the PEAR Archive_Tar package be installed:
|
||||
*
|
||||
* @see http://pear.php.net/package/Archive_Tar
|
||||
*/
|
||||
|
||||
// name and directory of package
|
||||
$pkgName = "JAMA";
|
||||
|
||||
// root of PHP/Math build directory
|
||||
$buildDir = substr(dirname(__FILE__), 0, -5 - strlen($pkgName));
|
||||
|
||||
// switch to PHP/Math build directory
|
||||
chdir($buildDir);
|
||||
|
||||
$tarName = "$pkgName.tar.gz";
|
||||
|
||||
$tarPath = $buildDir.$pkgName."/downloads/".$tarName;
|
||||
|
||||
if($_GET['op'] == "download") {
|
||||
|
||||
require_once('Archive/Tar.php');
|
||||
|
||||
$tar = new Archive_Tar($tarPath);
|
||||
|
||||
// create $pkgName archive under $pkgName folder
|
||||
$files = glob("$pkgName/*.php");
|
||||
$files = array_merge($files, glob("$pkgName/*.TXT"));
|
||||
$files = array_merge($files, glob("$pkgName/docs/*.php"));
|
||||
$files = array_merge($files, glob("$pkgName/docs/includes/*.php"));
|
||||
$files = array_merge($files, glob("$pkgName/examples/*.php"));
|
||||
$files = array_merge($files, glob("$pkgName/tests/*.php"));
|
||||
$files = array_merge($files, glob("$pkgName/utils/*.php"));
|
||||
|
||||
$tar->create($files);
|
||||
|
||||
// create the download url
|
||||
$webDir = substr($_SERVER['PHP_SELF'], 0, -18);
|
||||
$urlPath = "http://".$_SERVER['HTTP_HOST'].$webDir."/downloads";
|
||||
|
||||
// redirect to download url
|
||||
header("Location: $urlPath/$tarName");
|
||||
|
||||
}
|
||||
|
||||
include_once "includes/header.php";
|
||||
include_once "includes/navbar.php";
|
||||
?>
|
||||
<p>
|
||||
Download current version:
|
||||
</p>
|
||||
<ul>
|
||||
<li><a href='<?php echo $_SERVER['PHP_SELF']."?op=download"; ?>'><?php echo $tarName ?></a></li>
|
||||
</ul>
|
||||
<?php
|
||||
include_once "includes/footer.php";
|
||||
?>
|
||||
@ -1,166 +0,0 @@
|
||||
<?php
|
||||
include_once "includes/header.php";
|
||||
include_once "includes/navbar.php";
|
||||
?>
|
||||
<h2>Magic Square Example</h2>
|
||||
<p>
|
||||
The Jama distribution comes with a magic square example that is used to
|
||||
test and benchmark the LU, QR, SVD and symmetric Eig decompositions.
|
||||
The example outputs a multi-column table with these column headings:
|
||||
</p>
|
||||
|
||||
<table border='1' cellpadding='5' cellspacing='0' align='center'>
|
||||
<tr>
|
||||
<td><b>n</b></td>
|
||||
<td>Order of magic square.</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><b>trace</b></td>
|
||||
<td>Diagonal sum, should be the magic sum, (n^3 + n)/2.</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><b>max_eig</b></td>
|
||||
<td>Maximum eigenvalue of (A + A')/2, should equal trace.</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><b>rank</b></td>
|
||||
<td>Linear algebraic rank, should equal n if n is odd, be less than n if n is even.</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><b>cond</b></td>
|
||||
<td>L_2 condition number, ratio of singular values.</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><b>lu_res</b></td>
|
||||
<td>test of LU factorization, norm1(L*U-A(p,:))/(n*eps).</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><b>qr_res</b></td>
|
||||
<td>test of QR factorization, norm1(Q*R-A)/(n*eps).</td>
|
||||
</tr>
|
||||
</table>
|
||||
<p>
|
||||
Running the Java-based version of the matix square example produces these results:
|
||||
</p>
|
||||
|
||||
<table border='1' align='center'>
|
||||
<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>
|
||||
<tr>
|
||||
<td>3</td><td>15</td><td>15.000</td><td>3</td><td>4.330</td><td>0.000</td><td>11.333</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>4</td><td>34</td><td>34.000</td><td>3</td><td> Inf</td><td>0.000</td><td>13.500</td>
|
||||
<tr>
|
||||
<td>5</td><td>65</td><td>65.000</td><td>5</td><td>5.462</td><td>0.000</td><td>14.400</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>6</td><td>111</td><td>111.000</td><td>5</td><td> Inf</td><td>5.333</td><td>16.000</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>7</td><td>175</td><td>175.000</td><td>7</td><td>7.111</td><td>2.286</td><td>37.714</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>8</td><td>260</td><td>260.000</td><td>3</td><td> Inf</td><td>0.000</td><td>59.000</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>9</td><td>369</td><td>369.000</td><td>9</td><td>9.102</td><td>7.111</td><td>53.333</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>10</td><td>505</td><td>505.000</td><td>7</td><td> Inf</td><td>3.200</td><td>159.200</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>11</td><td>671</td><td>671.000</td><td>11</td><td>11.102</td><td>2.909</td><td>215.273</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>12</td><td>870</td><td>870.000</td><td>3</td><td> Inf</td><td>0.000</td><td>185.333</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>13</td><td>1105</td><td>1105.000</td><td>13</td><td>13.060</td><td>4.923</td><td>313.846</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>14</td><td>1379</td><td>1379.000</td><td>9</td><td> Inf</td><td>4.571</td><td>540.571</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>15</td><td>1695</td><td>1695.000</td><td>15</td><td>15.062</td><td>4.267</td><td>242.133</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>16</td><td>2056</td><td>2056.000</td><td>3</td><td> Inf</td><td>0.000</td><td>488.500</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>17</td><td>2465</td><td>2465.000</td><td>17</td><td>17.042</td><td>7.529</td><td>267.294</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>18</td><td>2925</td><td>2925.000</td><td>11</td><td> Inf</td><td>7.111</td><td>520.889</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>19</td><td>3439</td><td>3439.000</td><td>19</td><td>19.048</td><td>16.842</td><td>387.368</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>20</td><td>4010</td><td>4010.000</td><td>3</td><td> Inf</td><td>14.400</td><td>584.800</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>21</td><td>4641</td><td>4641.000</td><td>21</td><td>21.035</td><td>6.095</td><td>1158.095</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>22</td><td>5335</td><td>5335.000</td><td>13</td><td> Inf</td><td>6.545</td><td>1132.364</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>23</td><td>6095</td><td>6095.000</td><td>23</td><td>23.037</td><td>11.130</td><td>1268.870</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>24</td><td>6924</td><td>6924.000</td><td>3</td><td> Inf</td><td>10.667</td><td>827.500</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>25</td><td>7825</td><td>7825.000</td><td>25</td><td>25.029</td><td>35.840</td><td>1190.400</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>26</td><td>8801</td><td>8801.000</td><td>15</td><td> Inf</td><td>4.923</td><td>1859.077</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>27</td><td>9855</td><td>9855.000</td><td>27</td><td>27.032</td><td>37.926</td><td>1365.333</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>28</td><td>10990</td><td>10990.000</td><td>3</td><td> Inf</td><td>34.286</td><td>1365.714</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>29</td><td>12209</td><td>12209.000</td><td>29</td><td>29.025</td><td>30.897</td><td>1647.448</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>30</td><td>13515</td><td>13515.000</td><td>17</td><td> Inf</td><td>8.533</td><td>2571.733</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>31</td><td>14911</td><td>14911.000</td><td>31</td><td>31.027</td><td>33.032</td><td>1426.581</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>32</td><td>16400</td><td>16400.000</td><td>3</td><td> Inf</td><td>0.000</td><td>1600.125</td>
|
||||
</tr>
|
||||
</table>
|
||||
<center>Elapsed Time = 0.710 seconds</center>
|
||||
|
||||
<p>
|
||||
The magic square example does not fare well when <a href='../examples/MagicSquareExample.php'>run as a PHP script</a>. For a 32x32 matrix array
|
||||
it takes around a second to complete just the last row of computations in the above table.
|
||||
Hopefully this result will spur PHP developers to find optimizations and better attuned algorithms
|
||||
to speed things up. Matrix algebra is a great testing ground for ideas about time and memory
|
||||
performance optimation. Keep in perspective that PHP JAMA scripts are still plenty fast for use as
|
||||
a tool for learning about matrix algebra and quickly extending your knowledge with new scripts
|
||||
to apply knowledge.
|
||||
</p>
|
||||
|
||||
<p>
|
||||
To learn more about the subject of magic squares you can visit the <a href='http://mathforum.org/alejandre/magic.square.html'>Drexel Math Forum on Magic Squares</a>.
|
||||
You can also learn more by carefully examining the <code>MagicSquareExample.php</code> source code below.
|
||||
</p>
|
||||
|
||||
<?php
|
||||
highlight_file("../examples/MagicSquareExample.php");
|
||||
include_once "includes/footer.php";
|
||||
?>
|
||||
@ -1,14 +0,0 @@
|
||||
<div id="credits">
|
||||
<p>
|
||||
Brought to you by:
|
||||
</p>
|
||||
<ul>
|
||||
<li><a href="http://math.nist.gov/">National Institute of Standards and Technology</a></li>
|
||||
<li><a href="http://math.nist.gov/">MathWorks</a></li>
|
||||
<li><a href="http://math.nist.gov/javanumerics/jama/">JAMA : A Java Matrix Package</a></li>
|
||||
<li>Paul Meagher</li>
|
||||
<li>Michael Bommarito</li>
|
||||
<li>Lukasz Karapuda</li>
|
||||
<li>Bartek Matosiuk</li>
|
||||
</ul>
|
||||
</div>
|
||||
@ -1,2 +0,0 @@
|
||||
</body>
|
||||
</html>
|
||||
@ -1,11 +0,0 @@
|
||||
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
|
||||
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
|
||||
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
|
||||
<head>
|
||||
<title>JAMA v1.0.1</title>
|
||||
<meta name="description" content="JAMA v1.0.1 - Port of the proposed standard Java Matrix Library to PHP" />
|
||||
<meta name="robots" content="index, follow" />
|
||||
<meta name="keywords" content="php, matrix, matrix library, cholesky decomposition, eigenvalue decomposition, eigenvector, lu decomposition, qr decomposition, singular value decomposition" />
|
||||
<link rel="stylesheet" type="text/css" href="style.css" />
|
||||
</head>
|
||||
<body>
|
||||
@ -1,5 +0,0 @@
|
||||
<center>
|
||||
<hr />
|
||||
[ <a href="index.php">index.php</a> ] [ <a href="docs.php">docs.php</a> ] [ <a href="package.php">package.php</a> ] [ <a href="test.php">test.php</a> ] [ <a href="example.php">example.php</a> ] [ <a href="download.php">download.php</a> ]
|
||||
<hr />
|
||||
</center>
|
||||
@ -1,30 +0,0 @@
|
||||
<?php
|
||||
include_once "includes/header.php";
|
||||
include_once "includes/navbar.php";
|
||||
?>
|
||||
|
||||
<p>
|
||||
JAMA is a proposed standard matrix class for Java. The JAMA introduction
|
||||
describes "JAMA : A Java Matrix Package" in this way:
|
||||
</p>
|
||||
|
||||
<blockquote>
|
||||
JAMA is a basic linear algebra package for Java. It provides user-level classes for
|
||||
constructing and manipulating real, dense matrices. It is meant to provide sufficient
|
||||
functionality for routine problems, packaged in a way that is natural and understandable
|
||||
to non-experts. It is intended to serve as the standard matrix class for Java, and
|
||||
will be proposed as such to the Java Grande Forum and then to Sun. A straightforward
|
||||
public-domain reference implementation has been developed by the MathWorks and NIST as
|
||||
a strawman for such a class. We are releasing this version in order to obtain public
|
||||
comment. There is no guarantee that future versions of JAMA will be compatible with this one.
|
||||
</blockquote>
|
||||
|
||||
<p>
|
||||
The development team below has successfully ported the JAMA API to PHP. You can explore
|
||||
this site to learn more about this project and it's current development status.
|
||||
</p>
|
||||
|
||||
<?php
|
||||
include_once "includes/credits.php";
|
||||
include_once "includes/footer.php";
|
||||
?>
|
||||
@ -1,37 +0,0 @@
|
||||
<?php
|
||||
include_once "includes/header.php";
|
||||
include_once "includes/navbar.php";
|
||||
?>
|
||||
<p>
|
||||
Source Listing:
|
||||
</p>
|
||||
<ul>
|
||||
<?php
|
||||
chdir("../");
|
||||
$files = glob("*.php");
|
||||
$files = array_merge($files, glob("util/*.php"));
|
||||
foreach ($files as $fileName) {
|
||||
?>
|
||||
<li><a href="package.php?view=<?php echo sha1($fileName);?>"><?php echo $fileName;?></a> - <?php echo date ("F d Y - g:i a", filemtime($fileName));?></li>
|
||||
<?php
|
||||
}
|
||||
?>
|
||||
</ul>
|
||||
<?php
|
||||
if( isset($_REQUEST['view']) ) {
|
||||
$hash = $_REQUEST['view'];
|
||||
$n = array_search($hash, array_map(sha1, $files));
|
||||
$fileName = $files[$n];
|
||||
?>
|
||||
<hr />
|
||||
Viewing: <?php echo $fileName;?>
|
||||
<hr />
|
||||
<?php
|
||||
highlight_file($fileName);
|
||||
?>
|
||||
<hr />
|
||||
<?php
|
||||
}
|
||||
include_once "includes/footer.php";
|
||||
?>
|
||||
|
||||
@ -1,28 +0,0 @@
|
||||
<?php
|
||||
include_once "includes/header.php";
|
||||
include_once "includes/navbar.php";
|
||||
?>
|
||||
<p>
|
||||
The first script your should run when you install Jama is the TestMatrix.php script.
|
||||
</p>
|
||||
<p>
|
||||
This will run the unit tests for methods in the <code>Matrix.php</code> class. Because
|
||||
the Matrix.php class can be used to invoke all the decomposition methods the <code>TestMatrix.php</code>
|
||||
script is a test suite for the whole Jama package.
|
||||
</p>
|
||||
<p>
|
||||
The original <code>TestMatrix.java</code> code uses try/catch error handling. We will
|
||||
eventually create a build of JAMA that will take advantage of PHP5's new try/catch error
|
||||
handling capabilities. This will improve our ability to replicate all the unit tests that
|
||||
appeared in the original (except for some print methods that may not be worth porting).
|
||||
</p>
|
||||
<p>
|
||||
You can <a href='../test/TestMatrix.php'>run the TestMatrix.php script</a> to see what
|
||||
unit tests are currently implemented. The source of the <code>TestMatrix.php</code> script
|
||||
is provided below. It is worth studying carefully for an example of how to do matrix algebra
|
||||
programming with Jama.
|
||||
</p>
|
||||
<?php
|
||||
highlight_file("../test/TestMatrix.php");
|
||||
include_once "includes/footer.php";
|
||||
?>
|
||||
@ -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>";
|
||||
?>
|
||||
@ -1,415 +0,0 @@
|
||||
<?php
|
||||
|
||||
require_once "../Matrix.php";
|
||||
|
||||
class TestMatrix {
|
||||
|
||||
function TestMatrix() {
|
||||
|
||||
// define test variables
|
||||
|
||||
$errorCount = 0;
|
||||
$warningCount = 0;
|
||||
$columnwise = array(1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.);
|
||||
$rowwise = array(1.,4.,7.,10.,2.,5.,8.,11.,3.,6.,9.,12.);
|
||||
$avals = array(array(1.,4.,7.,10.),array(2.,5.,8.,11.),array(3.,6.,9.,12.));
|
||||
$rankdef = $avals;
|
||||
$tvals = array(array(1.,2.,3.),array(4.,5.,6.),array(7.,8.,9.),array(10.,11.,12.));
|
||||
$subavals = array(array(5.,8.,11.),array(6.,9.,12.));
|
||||
$rvals = array(array(1.,4.,7.),array(2.,5.,8.,11.),array(3.,6.,9.,12.));
|
||||
$pvals = array(array(1.,1.,1.),array(1.,2.,3.),array(1.,3.,6.));
|
||||
$ivals = array(array(1.,0.,0.,0.),array(0.,1.,0.,0.),array(0.,0.,1.,0.));
|
||||
$evals = array(array(0.,1.,0.,0.),array(1.,0.,2.e-7,0.),array(0.,-2.e-7,0.,1.),array(0.,0.,1.,0.));
|
||||
$square = array(array(166.,188.,210.),array(188.,214.,240.),array(210.,240.,270.));
|
||||
$sqSolution = array(array(13.),array(15.));
|
||||
$condmat = array(array(1.,3.),array(7.,9.));
|
||||
$rows = 3;
|
||||
$cols = 4;
|
||||
$invalidID = 5; /* should trigger bad shape for construction with val */
|
||||
$raggedr = 0; /* (raggedr,raggedc) should be out of bounds in ragged array */
|
||||
$raggedc = 4;
|
||||
$validID = 3; /* leading dimension of intended test Matrices */
|
||||
$nonconformld = 4; /* leading dimension which is valid, but nonconforming */
|
||||
$ib = 1; /* index ranges for sub Matrix */
|
||||
$ie = 2;
|
||||
$jb = 1;
|
||||
$je = 3;
|
||||
$rowindexset = array(1,2);
|
||||
$badrowindexset = array(1,3);
|
||||
$columnindexset = array(1,2,3);
|
||||
$badcolumnindexset = array(1,2,4);
|
||||
$columnsummax = 33.;
|
||||
$rowsummax = 30.;
|
||||
$sumofdiagonals = 15;
|
||||
$sumofsquares = 650;
|
||||
|
||||
/**
|
||||
* Test matrix methods
|
||||
*/
|
||||
|
||||
/**
|
||||
* Constructors and constructor-like methods:
|
||||
*
|
||||
* Matrix(double[], int)
|
||||
* Matrix(double[][])
|
||||
* Matrix(int, int)
|
||||
* Matrix(int, int, double)
|
||||
* Matrix(int, int, double[][])
|
||||
* constructWithCopy(double[][])
|
||||
* random(int,int)
|
||||
* identity(int)
|
||||
*/
|
||||
echo "<p>Testing constructors and constructor-like methods...</p>";
|
||||
|
||||
$A = new Matrix($columnwise, 3);
|
||||
if($A instanceof Matrix) {
|
||||
$this->try_success("Column-packed constructor...");
|
||||
} else
|
||||
$errorCount = $this->try_failure($errorCount, "Column-packed constructor...", "Unable to construct Matrix");
|
||||
|
||||
$T = new Matrix($tvals);
|
||||
if($T instanceof Matrix)
|
||||
$this->try_success("2D array constructor...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "2D array constructor...", "Unable to construct Matrix");
|
||||
|
||||
$A = new Matrix($columnwise, $validID);
|
||||
$B = new Matrix($avals);
|
||||
$tmp = $B->get(0,0);
|
||||
$avals[0][0] = 0.0;
|
||||
$C = $B->minus($A);
|
||||
$avals[0][0] = $tmp;
|
||||
$B = Matrix::constructWithCopy($avals);
|
||||
$tmp = $B->get(0,0);
|
||||
$avals[0][0] = 0.0;
|
||||
/** check that constructWithCopy behaves properly **/
|
||||
if ( ( $tmp - $B->get(0,0) ) != 0.0 )
|
||||
$errorCount = $this->try_failure($errorCount,"constructWithCopy... ","copy not effected... data visible outside");
|
||||
else
|
||||
$this->try_success("constructWithCopy... ","");
|
||||
|
||||
$I = new Matrix($ivals);
|
||||
if ( $this->checkMatrices($I,Matrix::identity(3,4)) )
|
||||
$this->try_success("identity... ","");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount,"identity... ","identity Matrix not successfully created");
|
||||
|
||||
/**
|
||||
* Access Methods:
|
||||
*
|
||||
* getColumnDimension()
|
||||
* getRowDimension()
|
||||
* getArray()
|
||||
* getArrayCopy()
|
||||
* getColumnPackedCopy()
|
||||
* getRowPackedCopy()
|
||||
* get(int,int)
|
||||
* getMatrix(int,int,int,int)
|
||||
* getMatrix(int,int,int[])
|
||||
* getMatrix(int[],int,int)
|
||||
* getMatrix(int[],int[])
|
||||
* set(int,int,double)
|
||||
* setMatrix(int,int,int,int,Matrix)
|
||||
* setMatrix(int,int,int[],Matrix)
|
||||
* setMatrix(int[],int,int,Matrix)
|
||||
* setMatrix(int[],int[],Matrix)
|
||||
*/
|
||||
print "<p>Testing access methods...</p>";
|
||||
|
||||
$B = new Matrix($avals);
|
||||
if($B->getRowDimension() == $rows)
|
||||
$this->try_success("getRowDimension...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "getRowDimension...");
|
||||
|
||||
if($B->getColumnDimension() == $cols)
|
||||
$this->try_success("getColumnDimension...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "getColumnDimension...");
|
||||
|
||||
$barray = $B->getArray();
|
||||
if($this->checkArrays($barray, $avals))
|
||||
$this->try_success("getArray...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "getArray...");
|
||||
|
||||
$bpacked = $B->getColumnPackedCopy();
|
||||
if($this->checkArrays($bpacked, $columnwise))
|
||||
$this->try_success("getColumnPackedCopy...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "getColumnPackedCopy...");
|
||||
|
||||
$bpacked = $B->getRowPackedCopy();
|
||||
if($this->checkArrays($bpacked, $rowwise))
|
||||
$this->try_success("getRowPackedCopy...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "getRowPackedCopy...");
|
||||
|
||||
/**
|
||||
* Array-like methods:
|
||||
* minus
|
||||
* minusEquals
|
||||
* plus
|
||||
* plusEquals
|
||||
* arrayLeftDivide
|
||||
* arrayLeftDivideEquals
|
||||
* arrayRightDivide
|
||||
* arrayRightDivideEquals
|
||||
* arrayTimes
|
||||
* arrayTimesEquals
|
||||
* uminus
|
||||
*/
|
||||
print "<p>Testing array-like methods...</p>";
|
||||
|
||||
/**
|
||||
* I/O methods:
|
||||
* read
|
||||
* print
|
||||
* serializable:
|
||||
* writeObject
|
||||
* readObject
|
||||
*/
|
||||
print "<p>Testing I/O methods...</p>";
|
||||
|
||||
/**
|
||||
* Test linear algebra methods
|
||||
*/
|
||||
echo "<p>Testing linear algebra methods...<p>";
|
||||
|
||||
$A = new Matrix($columnwise, 3);
|
||||
if( $this->checkMatrices($A->transpose(), $T) )
|
||||
$this->try_success("Transpose check...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "Transpose check...", "Matrices are not equal");
|
||||
|
||||
if($this->checkScalars($A->norm1(), $columnsummax))
|
||||
$this->try_success("Maximum column sum...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "Maximum column sum...", "Incorrect: " . $A->norm1() . " != " . $columnsummax);
|
||||
|
||||
if($this->checkScalars($A->normInf(), $rowsummax))
|
||||
$this->try_success("Maximum row sum...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "Maximum row sum...", "Incorrect: " . $A->normInf() . " != " . $rowsummax );
|
||||
|
||||
if($this->checkScalars($A->normF(), sqrt($sumofsquares)))
|
||||
$this->try_success("Frobenius norm...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "Frobenius norm...", "Incorrect:" . $A->normF() . " != " . sqrt($sumofsquares));
|
||||
|
||||
if($this->checkScalars($A->trace(), $sumofdiagonals))
|
||||
$this->try_success("Matrix trace...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "Matrix trace...", "Incorrect: " . $A->trace() . " != " . $sumofdiagonals);
|
||||
|
||||
$B = $A->getMatrix(0, $A->getRowDimension(), 0, $A->getRowDimension());
|
||||
if( $B->det() == 0 )
|
||||
$this->try_success("Matrix determinant...");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "Matrix determinant...", "Incorrect: " . $B->det() . " != " . 0);
|
||||
|
||||
$A = new Matrix($columnwise,3);
|
||||
$SQ = new Matrix($square);
|
||||
if ($this->checkMatrices($SQ, $A->times($A->transpose())))
|
||||
$this->try_success("times(Matrix)...");
|
||||
else {
|
||||
$errorCount = $this->try_failure($errorCount, "times(Matrix)...", "Unable to multiply matrices");
|
||||
$SQ->toHTML();
|
||||
$AT->toHTML();
|
||||
}
|
||||
|
||||
$A = new Matrix($columnwise, 4);
|
||||
|
||||
$QR = $A->qr();
|
||||
$R = $QR->getR();
|
||||
$Q = $QR->getQ();
|
||||
if($this->checkMatrices($A, $Q->times($R)))
|
||||
$this->try_success("QRDecomposition...","");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount,"QRDecomposition...","incorrect qr decomposition calculation");
|
||||
|
||||
$A = new Matrix($columnwise, 4);
|
||||
$SVD = $A->svd();
|
||||
$U = $SVD->getU();
|
||||
$S = $SVD->getS();
|
||||
$V = $SVD->getV();
|
||||
if ($this->checkMatrices($A, $U->times($S->times($V->transpose()))))
|
||||
$this->try_success("SingularValueDecomposition...","");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount,"SingularValueDecomposition...","incorrect singular value decomposition calculation");
|
||||
|
||||
$n = $A->getColumnDimension();
|
||||
$A = $A->getMatrix(0,$n-1,0,$n-1);
|
||||
$A->set(0,0,0.);
|
||||
|
||||
$LU = $A->lu();
|
||||
$L = $LU->getL();
|
||||
if ( $this->checkMatrices($A->getMatrix($LU->getPivot(),0,$n-1), $L->times($LU->getU())) )
|
||||
$this->try_success("LUDecomposition...","");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount,"LUDecomposition...","incorrect LU decomposition calculation");
|
||||
|
||||
$X = $A->inverse();
|
||||
if ( $this->checkMatrices($A->times($X),Matrix::identity(3,3)) )
|
||||
$this->try_success("inverse()...","");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount, "inverse()...","incorrect inverse calculation");
|
||||
|
||||
$DEF = new Matrix($rankdef);
|
||||
if($this->checkScalars($DEF->rank(), min($DEF->getRowDimension(), $DEF->getColumnDimension())-1))
|
||||
$this->try_success("Rank...");
|
||||
else
|
||||
$this->try_failure("Rank...", "incorrect rank calculation");
|
||||
|
||||
$B = new Matrix($condmat);
|
||||
$SVD = $B->svd();
|
||||
$singularvalues = $SVD->getSingularValues();
|
||||
if($this->checkScalars($B->cond(), $singularvalues[0]/$singularvalues[min($B->getRowDimension(), $B->getColumnDimension())-1]))
|
||||
$this->try_success("Condition number...");
|
||||
else
|
||||
$this->try_failure("Condition number...", "incorrect condition number calculation");
|
||||
|
||||
$SUB = new Matrix($subavals);
|
||||
$O = new Matrix($SUB->getRowDimension(),1,1.0);
|
||||
$SOL = new Matrix($sqSolution);
|
||||
$SQ = $SUB->getMatrix(0,$SUB->getRowDimension()-1,0,$SUB->getRowDimension()-1);
|
||||
if ( $this->checkMatrices($SQ->solve($SOL),$O) )
|
||||
$this->try_success("solve()...","");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount,"solve()...","incorrect lu solve calculation");
|
||||
|
||||
$A = new Matrix($pvals);
|
||||
$Chol = $A->chol();
|
||||
$L = $Chol->getL();
|
||||
if ( $this->checkMatrices($A, $L->times($L->transpose())) )
|
||||
$this->try_success("CholeskyDecomposition...","");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount,"CholeskyDecomposition...","incorrect Cholesky decomposition calculation");
|
||||
|
||||
$X = $Chol->solve(Matrix::identity(3,3));
|
||||
if ( $this->checkMatrices($A->times($X), Matrix::identity(3,3)) )
|
||||
$this->try_success("CholeskyDecomposition solve()...","");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount,"CholeskyDecomposition solve()...","incorrect Choleskydecomposition solve calculation");
|
||||
|
||||
$Eig = $A->eig();
|
||||
$D = $Eig->getD();
|
||||
$V = $Eig->getV();
|
||||
if( $this->checkMatrices($A->times($V),$V->times($D)) )
|
||||
$this->try_success("EigenvalueDecomposition (symmetric)...","");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount,"EigenvalueDecomposition (symmetric)...","incorrect symmetric Eigenvalue decomposition calculation");
|
||||
|
||||
$A = new Matrix($evals);
|
||||
$Eig = $A->eig();
|
||||
$D = $Eig->getD();
|
||||
$V = $Eig->getV();
|
||||
if ( $this->checkMatrices($A->times($V),$V->times($D)) )
|
||||
$this->try_success("EigenvalueDecomposition (nonsymmetric)...","");
|
||||
else
|
||||
$errorCount = $this->try_failure($errorCount,"EigenvalueDecomposition (nonsymmetric)...","incorrect nonsymmetric Eigenvalue decomposition calculation");
|
||||
|
||||
print("<b>{$errorCount} total errors</b>.");
|
||||
}
|
||||
|
||||
/**
|
||||
* Print appropriate messages for successful outcome try
|
||||
* @param string $s
|
||||
* @param string $e
|
||||
*/
|
||||
function try_success($s, $e = "") {
|
||||
print "> ". $s ."success<br />";
|
||||
if ($e != "")
|
||||
print "> Message: ". $e ."<br />";
|
||||
}
|
||||
|
||||
/**
|
||||
* Print appropriate messages for unsuccessful outcome try
|
||||
* @param int $count
|
||||
* @param string $s
|
||||
* @param string $e
|
||||
* @return int incremented counter
|
||||
*/
|
||||
function try_failure($count, $s, $e="") {
|
||||
print "> ". $s ."*** failure ***<br />> Message: ". $e ."<br />";
|
||||
return ++$count;
|
||||
}
|
||||
|
||||
/**
|
||||
* Print appropriate messages for unsuccessful outcome try
|
||||
* @param int $count
|
||||
* @param string $s
|
||||
* @param string $e
|
||||
* @return int incremented counter
|
||||
*/
|
||||
function try_warning($count, $s, $e="") {
|
||||
print "> ". $s ."*** warning ***<br />> Message: ". $e ."<br />";
|
||||
return ++$count;
|
||||
}
|
||||
|
||||
/**
|
||||
* Check magnitude of difference of "scalars".
|
||||
* @param float $x
|
||||
* @param float $y
|
||||
*/
|
||||
function checkScalars($x, $y) {
|
||||
$eps = pow(2.0,-52.0);
|
||||
if ($x == 0 & abs($y) < 10*$eps) return;
|
||||
if ($y == 0 & abs($x) < 10*$eps) return;
|
||||
if (abs($x-$y) > 10 * $eps * max(abs($x),abs($y)))
|
||||
return false;
|
||||
else
|
||||
return true;
|
||||
}
|
||||
|
||||
/**
|
||||
* Check norm of difference of "vectors".
|
||||
* @param float $x[]
|
||||
* @param float $y[]
|
||||
*/
|
||||
function checkVectors($x, $y) {
|
||||
$nx = count($x);
|
||||
$ny = count($y);
|
||||
if ($nx == $ny)
|
||||
for($i=0; $i < $nx; ++$i)
|
||||
$this->checkScalars($x[$i],$y[$i]);
|
||||
else
|
||||
die("Attempt to compare vectors of different lengths");
|
||||
}
|
||||
|
||||
/**
|
||||
* Check norm of difference of "arrays".
|
||||
* @param float $x[][]
|
||||
* @param float $y[][]
|
||||
*/
|
||||
function checkArrays($x, $y) {
|
||||
$A = new Matrix($x);
|
||||
$B = new Matrix($y);
|
||||
return $this->checkMatrices($A,$B);
|
||||
}
|
||||
|
||||
/**
|
||||
* Check norm of difference of "matrices".
|
||||
* @param matrix $X
|
||||
* @param matrix $Y
|
||||
*/
|
||||
function checkMatrices($X = null, $Y = null) {
|
||||
if( $X == null || $Y == null )
|
||||
return false;
|
||||
|
||||
$eps = pow(2.0,-52.0);
|
||||
if ($X->norm1() == 0. & $Y->norm1() < 10*$eps) return true;
|
||||
if ($Y->norm1() == 0. & $X->norm1() < 10*$eps) return true;
|
||||
|
||||
$A = $X->minus($Y);
|
||||
|
||||
if ($A->norm1() > 1000 * $eps * max($X->norm1(),$Y->norm1()))
|
||||
die("The norm of (X-Y) is too large: ".$A->norm1());
|
||||
else
|
||||
return true;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
$test = new TestMatrix;
|
||||
?>
|
||||
@ -1,82 +0,0 @@
|
||||
<?php
|
||||
/**
|
||||
* @package JAMA
|
||||
*
|
||||
* Error handling
|
||||
* @author Michael Bommarito
|
||||
* @version 01292005
|
||||
*/
|
||||
|
||||
//Language constant
|
||||
define('JAMALANG', 'EN');
|
||||
|
||||
|
||||
//All errors may be defined by the following format:
|
||||
//define('ExceptionName', N);
|
||||
//$error['lang'][ExceptionName] = 'Error message';
|
||||
$error = array();
|
||||
|
||||
/*
|
||||
I've used Babelfish and a little poor knowledge of Romance/Germanic languages for the translations here.
|
||||
Feel free to correct anything that looks amiss to you.
|
||||
*/
|
||||
|
||||
define('PolymorphicArgumentException', -1);
|
||||
$error['EN'][PolymorphicArgumentException] = "Invalid argument pattern for polymorphic function.";
|
||||
$error['FR'][PolymorphicArgumentException] = "Modèle inadmissible d'argument pour la fonction polymorphe.".
|
||||
$error['DE'][PolymorphicArgumentException] = "Unzulässiges Argumentmuster für polymorphe Funktion.";
|
||||
|
||||
define('ArgumentTypeException', -2);
|
||||
$error['EN'][ArgumentTypeException] = "Invalid argument type.";
|
||||
$error['FR'][ArgumentTypeException] = "Type inadmissible d'argument.";
|
||||
$error['DE'][ArgumentTypeException] = "Unzulässige Argumentart.";
|
||||
|
||||
define('ArgumentBoundsException', -3);
|
||||
$error['EN'][ArgumentBoundsException] = "Invalid argument range.";
|
||||
$error['FR'][ArgumentBoundsException] = "Gamme inadmissible d'argument.";
|
||||
$error['DE'][ArgumentBoundsException] = "Unzulässige Argumentstrecke.";
|
||||
|
||||
define('MatrixDimensionException', -4);
|
||||
$error['EN'][MatrixDimensionException] = "Matrix dimensions are not equal.";
|
||||
$error['FR'][MatrixDimensionException] = "Les dimensions de Matrix ne sont pas égales.";
|
||||
$error['DE'][MatrixDimensionException] = "Matrixmaße sind nicht gleich.";
|
||||
|
||||
define('PrecisionLossException', -5);
|
||||
$error['EN'][PrecisionLossException] = "Significant precision loss detected.";
|
||||
$error['FR'][PrecisionLossException] = "Perte significative de précision détectée.";
|
||||
$error['DE'][PrecisionLossException] = "Bedeutender Präzision Verlust ermittelte.";
|
||||
|
||||
define('MatrixSPDException', -6);
|
||||
$error['EN'][MatrixSPDException] = "Can only perform operation on symmetric positive definite matrix.";
|
||||
$error['FR'][MatrixSPDException] = "Perte significative de précision détectée.";
|
||||
$error['DE'][MatrixSPDException] = "Bedeutender Präzision Verlust ermittelte.";
|
||||
|
||||
define('MatrixSingularException', -7);
|
||||
$error['EN'][MatrixSingularException] = "Can only perform operation on singular matrix.";
|
||||
|
||||
define('MatrixRankException', -8);
|
||||
$error['EN'][MatrixRankException] = "Can only perform operation on full-rank matrix.";
|
||||
|
||||
define('ArrayLengthException', -9);
|
||||
$error['EN'][ArrayLengthException] = "Array length must be a multiple of m.";
|
||||
|
||||
define('RowLengthException', -10);
|
||||
$error['EN'][RowLengthException] = "All rows must have the same length.";
|
||||
|
||||
/**
|
||||
* Custom error handler
|
||||
* @param int $num Error number
|
||||
*/
|
||||
function JAMAError($errorNumber = null) {
|
||||
global $error;
|
||||
|
||||
if (isset($errorNumber)) {
|
||||
if (isset($error[JAMALANG][$errorNumber])) {
|
||||
return $error[JAMALANG][$errorNumber];
|
||||
} else {
|
||||
return $error['EN'][$errorNumber];
|
||||
}
|
||||
} else {
|
||||
return ("Invalid argument to JAMAError()");
|
||||
}
|
||||
}
|
||||
@ -1,43 +0,0 @@
|
||||
<?php
|
||||
/**
|
||||
* @package JAMA
|
||||
*
|
||||
* Pythagorean Theorem:
|
||||
*
|
||||
* a = 3
|
||||
* b = 4
|
||||
* r = sqrt(square(a) + square(b))
|
||||
* r = 5
|
||||
*
|
||||
* r = sqrt(a^2 + b^2) without under/overflow.
|
||||
*/
|
||||
function hypo($a, $b) {
|
||||
if (abs($a) > abs($b)) {
|
||||
$r = $b / $a;
|
||||
$r = abs($a) * sqrt(1 + $r * $r);
|
||||
} elseif ($b != 0) {
|
||||
$r = $a / $b;
|
||||
$r = abs($b) * sqrt(1 + $r * $r);
|
||||
} else {
|
||||
$r = 0.0;
|
||||
}
|
||||
return $r;
|
||||
} // function hypo()
|
||||
|
||||
|
||||
/**
|
||||
* Mike Bommarito's version.
|
||||
* Compute n-dimensional hyotheneuse.
|
||||
*
|
||||
function hypot() {
|
||||
$s = 0;
|
||||
foreach (func_get_args() as $d) {
|
||||
if (is_numeric($d)) {
|
||||
$s += pow($d, 2);
|
||||
} else {
|
||||
throw new Exception(JAMAError(ArgumentTypeException));
|
||||
}
|
||||
}
|
||||
return sqrt($s);
|
||||
}
|
||||
*/
|
||||
Reference in New Issue
Block a user