/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php |
---|
New file |
0,0 → 1,862 |
<?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 |
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/Matrix.php |
---|
New file |
0,0 → 1,1059 |
<?php |
/** |
* @package JAMA |
*/ |
/** PHPExcel root directory */ |
if (!defined('PHPEXCEL_ROOT')) { |
/** |
* @ignore |
*/ |
define('PHPEXCEL_ROOT', dirname(__FILE__) . '/../../../'); |
require(PHPEXCEL_ROOT . 'PHPExcel/Autoloader.php'); |
} |
/* |
* Matrix class |
* |
* @author Paul Meagher |
* @author Michael Bommarito |
* @author Lukasz Karapuda |
* @author Bartek Matosiuk |
* @version 1.8 |
* @license PHP v3.0 |
* @see http://math.nist.gov/javanumerics/jama/ |
*/ |
class PHPExcel_Shared_JAMA_Matrix { |
const PolymorphicArgumentException = "Invalid argument pattern for polymorphic function."; |
const ArgumentTypeException = "Invalid argument type."; |
const ArgumentBoundsException = "Invalid argument range."; |
const MatrixDimensionException = "Matrix dimensions are not equal."; |
const ArrayLengthException = "Array length must be a multiple of m."; |
/** |
* Matrix storage |
* |
* @var array |
* @access public |
*/ |
public $A = array(); |
/** |
* Matrix row dimension |
* |
* @var int |
* @access private |
*/ |
private $m; |
/** |
* Matrix column dimension |
* |
* @var int |
* @access private |
*/ |
private $n; |
/** |
* Polymorphic constructor |
* |
* As PHP has no support for polymorphic constructors, we hack our own sort of polymorphism using func_num_args, func_get_arg, and gettype. In essence, we're just implementing a simple RTTI filter and calling the appropriate constructor. |
*/ |
public function __construct() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
//Rectangular matrix - m x n initialized from 2D array |
case 'array': |
$this->m = count($args[0]); |
$this->n = count($args[0][0]); |
$this->A = $args[0]; |
break; |
//Square matrix - n x n |
case 'integer': |
$this->m = $args[0]; |
$this->n = $args[0]; |
$this->A = array_fill(0, $this->m, array_fill(0, $this->n, 0)); |
break; |
//Rectangular matrix - m x n |
case 'integer,integer': |
$this->m = $args[0]; |
$this->n = $args[1]; |
$this->A = array_fill(0, $this->m, array_fill(0, $this->n, 0)); |
break; |
//Rectangular matrix - m x n initialized from packed array |
case 'array,integer': |
$this->m = $args[1]; |
if ($this->m != 0) { |
$this->n = count($args[0]) / $this->m; |
} else { |
$this->n = 0; |
} |
if (($this->m * $this->n) == count($args[0])) { |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$this->A[$i][$j] = $args[0][$i + $j * $this->m]; |
} |
} |
} else { |
throw new PHPExcel_Calculation_Exception(self::ArrayLengthException); |
} |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function __construct() |
/** |
* getArray |
* |
* @return array Matrix array |
*/ |
public function getArray() { |
return $this->A; |
} // function getArray() |
/** |
* getRowDimension |
* |
* @return int Row dimension |
*/ |
public function getRowDimension() { |
return $this->m; |
} // function getRowDimension() |
/** |
* getColumnDimension |
* |
* @return int Column dimension |
*/ |
public function getColumnDimension() { |
return $this->n; |
} // function getColumnDimension() |
/** |
* get |
* |
* Get the i,j-th element of the matrix. |
* @param int $i Row position |
* @param int $j Column position |
* @return mixed Element (int/float/double) |
*/ |
public function get($i = null, $j = null) { |
return $this->A[$i][$j]; |
} // function get() |
/** |
* getMatrix |
* |
* Get a submatrix |
* @param int $i0 Initial row index |
* @param int $iF Final row index |
* @param int $j0 Initial column index |
* @param int $jF Final column index |
* @return Matrix Submatrix |
*/ |
public function getMatrix() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
//A($i0...; $j0...) |
case 'integer,integer': |
list($i0, $j0) = $args; |
if ($i0 >= 0) { $m = $this->m - $i0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
if ($j0 >= 0) { $n = $this->n - $j0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n); |
for($i = $i0; $i < $this->m; ++$i) { |
for($j = $j0; $j < $this->n; ++$j) { |
$R->set($i, $j, $this->A[$i][$j]); |
} |
} |
return $R; |
break; |
//A($i0...$iF; $j0...$jF) |
case 'integer,integer,integer,integer': |
list($i0, $iF, $j0, $jF) = $args; |
if (($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0)) { $m = $iF - $i0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
if (($jF > $j0) && ($this->n >= $jF) && ($j0 >= 0)) { $n = $jF - $j0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
$R = new PHPExcel_Shared_JAMA_Matrix($m+1, $n+1); |
for($i = $i0; $i <= $iF; ++$i) { |
for($j = $j0; $j <= $jF; ++$j) { |
$R->set($i - $i0, $j - $j0, $this->A[$i][$j]); |
} |
} |
return $R; |
break; |
//$R = array of row indices; $C = array of column indices |
case 'array,array': |
list($RL, $CL) = $args; |
if (count($RL) > 0) { $m = count($RL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
if (count($CL) > 0) { $n = count($CL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n); |
for($i = 0; $i < $m; ++$i) { |
for($j = 0; $j < $n; ++$j) { |
$R->set($i - $i0, $j - $j0, $this->A[$RL[$i]][$CL[$j]]); |
} |
} |
return $R; |
break; |
//$RL = array of row indices; $CL = array of column indices |
case 'array,array': |
list($RL, $CL) = $args; |
if (count($RL) > 0) { $m = count($RL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
if (count($CL) > 0) { $n = count($CL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n); |
for($i = 0; $i < $m; ++$i) { |
for($j = 0; $j < $n; ++$j) { |
$R->set($i, $j, $this->A[$RL[$i]][$CL[$j]]); |
} |
} |
return $R; |
break; |
//A($i0...$iF); $CL = array of column indices |
case 'integer,integer,array': |
list($i0, $iF, $CL) = $args; |
if (($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0)) { $m = $iF - $i0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
if (count($CL) > 0) { $n = count($CL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n); |
for($i = $i0; $i < $iF; ++$i) { |
for($j = 0; $j < $n; ++$j) { |
$R->set($i - $i0, $j, $this->A[$RL[$i]][$j]); |
} |
} |
return $R; |
break; |
//$RL = array of row indices |
case 'array,integer,integer': |
list($RL, $j0, $jF) = $args; |
if (count($RL) > 0) { $m = count($RL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
if (($jF >= $j0) && ($this->n >= $jF) && ($j0 >= 0)) { $n = $jF - $j0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); } |
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n+1); |
for($i = 0; $i < $m; ++$i) { |
for($j = $j0; $j <= $jF; ++$j) { |
$R->set($i, $j - $j0, $this->A[$RL[$i]][$j]); |
} |
} |
return $R; |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function getMatrix() |
/** |
* checkMatrixDimensions |
* |
* Is matrix B the same size? |
* @param Matrix $B Matrix B |
* @return boolean |
*/ |
public function checkMatrixDimensions($B = null) { |
if ($B instanceof PHPExcel_Shared_JAMA_Matrix) { |
if (($this->m == $B->getRowDimension()) && ($this->n == $B->getColumnDimension())) { |
return true; |
} else { |
throw new PHPExcel_Calculation_Exception(self::MatrixDimensionException); |
} |
} else { |
throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); |
} |
} // function checkMatrixDimensions() |
/** |
* set |
* |
* Set the i,j-th element of the matrix. |
* @param int $i Row position |
* @param int $j Column position |
* @param mixed $c Int/float/double value |
* @return mixed Element (int/float/double) |
*/ |
public function set($i = null, $j = null, $c = null) { |
// Optimized set version just has this |
$this->A[$i][$j] = $c; |
} // function set() |
/** |
* identity |
* |
* Generate an identity matrix. |
* @param int $m Row dimension |
* @param int $n Column dimension |
* @return Matrix Identity matrix |
*/ |
public function identity($m = null, $n = null) { |
return $this->diagonal($m, $n, 1); |
} // function identity() |
/** |
* diagonal |
* |
* Generate a diagonal matrix |
* @param int $m Row dimension |
* @param int $n Column dimension |
* @param mixed $c Diagonal value |
* @return Matrix Diagonal matrix |
*/ |
public function diagonal($m = null, $n = null, $c = 1) { |
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n); |
for($i = 0; $i < $m; ++$i) { |
$R->set($i, $i, $c); |
} |
return $R; |
} // function diagonal() |
/** |
* getMatrixByRow |
* |
* Get a submatrix by row index/range |
* @param int $i0 Initial row index |
* @param int $iF Final row index |
* @return Matrix Submatrix |
*/ |
public function getMatrixByRow($i0 = null, $iF = null) { |
if (is_int($i0)) { |
if (is_int($iF)) { |
return $this->getMatrix($i0, 0, $iF + 1, $this->n); |
} else { |
return $this->getMatrix($i0, 0, $i0 + 1, $this->n); |
} |
} else { |
throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); |
} |
} // function getMatrixByRow() |
/** |
* getMatrixByCol |
* |
* Get a submatrix by column index/range |
* @param int $i0 Initial column index |
* @param int $iF Final column index |
* @return Matrix Submatrix |
*/ |
public function getMatrixByCol($j0 = null, $jF = null) { |
if (is_int($j0)) { |
if (is_int($jF)) { |
return $this->getMatrix(0, $j0, $this->m, $jF + 1); |
} else { |
return $this->getMatrix(0, $j0, $this->m, $j0 + 1); |
} |
} else { |
throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); |
} |
} // function getMatrixByCol() |
/** |
* transpose |
* |
* Tranpose matrix |
* @return Matrix Transposed matrix |
*/ |
public function transpose() { |
$R = new PHPExcel_Shared_JAMA_Matrix($this->n, $this->m); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$R->set($j, $i, $this->A[$i][$j]); |
} |
} |
return $R; |
} // function transpose() |
/** |
* trace |
* |
* Sum of diagonal elements |
* @return float Sum of diagonal elements |
*/ |
public function trace() { |
$s = 0; |
$n = min($this->m, $this->n); |
for($i = 0; $i < $n; ++$i) { |
$s += $this->A[$i][$i]; |
} |
return $s; |
} // function trace() |
/** |
* uminus |
* |
* Unary minus matrix -A |
* @return Matrix Unary minus matrix |
*/ |
public function uminus() { |
} // function uminus() |
/** |
* plus |
* |
* A + B |
* @param mixed $B Matrix/Array |
* @return Matrix Sum |
*/ |
public function plus() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$M->set($i, $j, $M->get($i, $j) + $this->A[$i][$j]); |
} |
} |
return $M; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function plus() |
/** |
* plusEquals |
* |
* A = A + B |
* @param mixed $B Matrix/Array |
* @return Matrix Sum |
*/ |
public function plusEquals() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$validValues = True; |
$value = $M->get($i, $j); |
if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) { |
$this->A[$i][$j] = trim($this->A[$i][$j],'"'); |
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($this->A[$i][$j]); |
} |
if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) { |
$value = trim($value,'"'); |
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($value); |
} |
if ($validValues) { |
$this->A[$i][$j] += $value; |
} else { |
$this->A[$i][$j] = PHPExcel_Calculation_Functions::NaN(); |
} |
} |
} |
return $this; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function plusEquals() |
/** |
* minus |
* |
* A - B |
* @param mixed $B Matrix/Array |
* @return Matrix Sum |
*/ |
public function minus() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$M->set($i, $j, $M->get($i, $j) - $this->A[$i][$j]); |
} |
} |
return $M; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function minus() |
/** |
* minusEquals |
* |
* A = A - B |
* @param mixed $B Matrix/Array |
* @return Matrix Sum |
*/ |
public function minusEquals() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$validValues = True; |
$value = $M->get($i, $j); |
if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) { |
$this->A[$i][$j] = trim($this->A[$i][$j],'"'); |
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($this->A[$i][$j]); |
} |
if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) { |
$value = trim($value,'"'); |
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($value); |
} |
if ($validValues) { |
$this->A[$i][$j] -= $value; |
} else { |
$this->A[$i][$j] = PHPExcel_Calculation_Functions::NaN(); |
} |
} |
} |
return $this; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function minusEquals() |
/** |
* arrayTimes |
* |
* Element-by-element multiplication |
* Cij = Aij * Bij |
* @param mixed $B Matrix/Array |
* @return Matrix Matrix Cij |
*/ |
public function arrayTimes() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$M->set($i, $j, $M->get($i, $j) * $this->A[$i][$j]); |
} |
} |
return $M; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function arrayTimes() |
/** |
* arrayTimesEquals |
* |
* Element-by-element multiplication |
* Aij = Aij * Bij |
* @param mixed $B Matrix/Array |
* @return Matrix Matrix Aij |
*/ |
public function arrayTimesEquals() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$validValues = True; |
$value = $M->get($i, $j); |
if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) { |
$this->A[$i][$j] = trim($this->A[$i][$j],'"'); |
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($this->A[$i][$j]); |
} |
if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) { |
$value = trim($value,'"'); |
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($value); |
} |
if ($validValues) { |
$this->A[$i][$j] *= $value; |
} else { |
$this->A[$i][$j] = PHPExcel_Calculation_Functions::NaN(); |
} |
} |
} |
return $this; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function arrayTimesEquals() |
/** |
* arrayRightDivide |
* |
* Element-by-element right division |
* A / B |
* @param Matrix $B Matrix B |
* @return Matrix Division result |
*/ |
public function arrayRightDivide() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$validValues = True; |
$value = $M->get($i, $j); |
if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) { |
$this->A[$i][$j] = trim($this->A[$i][$j],'"'); |
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($this->A[$i][$j]); |
} |
if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) { |
$value = trim($value,'"'); |
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($value); |
} |
if ($validValues) { |
if ($value == 0) { |
// Trap for Divide by Zero error |
$M->set($i, $j, '#DIV/0!'); |
} else { |
$M->set($i, $j, $this->A[$i][$j] / $value); |
} |
} else { |
$M->set($i, $j, PHPExcel_Calculation_Functions::NaN()); |
} |
} |
} |
return $M; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function arrayRightDivide() |
/** |
* arrayRightDivideEquals |
* |
* Element-by-element right division |
* Aij = Aij / Bij |
* @param mixed $B Matrix/Array |
* @return Matrix Matrix Aij |
*/ |
public function arrayRightDivideEquals() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$this->A[$i][$j] = $this->A[$i][$j] / $M->get($i, $j); |
} |
} |
return $M; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function arrayRightDivideEquals() |
/** |
* arrayLeftDivide |
* |
* Element-by-element Left division |
* A / B |
* @param Matrix $B Matrix B |
* @return Matrix Division result |
*/ |
public function arrayLeftDivide() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$M->set($i, $j, $M->get($i, $j) / $this->A[$i][$j]); |
} |
} |
return $M; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function arrayLeftDivide() |
/** |
* arrayLeftDivideEquals |
* |
* Element-by-element Left division |
* Aij = Aij / Bij |
* @param mixed $B Matrix/Array |
* @return Matrix Matrix Aij |
*/ |
public function arrayLeftDivideEquals() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$this->A[$i][$j] = $M->get($i, $j) / $this->A[$i][$j]; |
} |
} |
return $M; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function arrayLeftDivideEquals() |
/** |
* times |
* |
* Matrix multiplication |
* @param mixed $n Matrix/Array/Scalar |
* @return Matrix Product |
*/ |
public function times() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $B = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
if ($this->n == $B->m) { |
$C = new PHPExcel_Shared_JAMA_Matrix($this->m, $B->n); |
for($j = 0; $j < $B->n; ++$j) { |
for ($k = 0; $k < $this->n; ++$k) { |
$Bcolj[$k] = $B->A[$k][$j]; |
} |
for($i = 0; $i < $this->m; ++$i) { |
$Arowi = $this->A[$i]; |
$s = 0; |
for($k = 0; $k < $this->n; ++$k) { |
$s += $Arowi[$k] * $Bcolj[$k]; |
} |
$C->A[$i][$j] = $s; |
} |
} |
return $C; |
} else { |
throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionMismatch)); |
} |
break; |
case 'array': |
$B = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
if ($this->n == $B->m) { |
$C = new PHPExcel_Shared_JAMA_Matrix($this->m, $B->n); |
for($i = 0; $i < $C->m; ++$i) { |
for($j = 0; $j < $C->n; ++$j) { |
$s = "0"; |
for($k = 0; $k < $C->n; ++$k) { |
$s += $this->A[$i][$k] * $B->A[$k][$j]; |
} |
$C->A[$i][$j] = $s; |
} |
} |
return $C; |
} else { |
throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionMismatch)); |
} |
return $M; |
break; |
case 'integer': |
$C = new PHPExcel_Shared_JAMA_Matrix($this->A); |
for($i = 0; $i < $C->m; ++$i) { |
for($j = 0; $j < $C->n; ++$j) { |
$C->A[$i][$j] *= $args[0]; |
} |
} |
return $C; |
break; |
case 'double': |
$C = new PHPExcel_Shared_JAMA_Matrix($this->m, $this->n); |
for($i = 0; $i < $C->m; ++$i) { |
for($j = 0; $j < $C->n; ++$j) { |
$C->A[$i][$j] = $args[0] * $this->A[$i][$j]; |
} |
} |
return $C; |
break; |
case 'float': |
$C = new PHPExcel_Shared_JAMA_Matrix($this->A); |
for($i = 0; $i < $C->m; ++$i) { |
for($j = 0; $j < $C->n; ++$j) { |
$C->A[$i][$j] *= $args[0]; |
} |
} |
return $C; |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function times() |
/** |
* power |
* |
* A = A ^ B |
* @param mixed $B Matrix/Array |
* @return Matrix Sum |
*/ |
public function power() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
break; |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$validValues = True; |
$value = $M->get($i, $j); |
if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) { |
$this->A[$i][$j] = trim($this->A[$i][$j],'"'); |
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($this->A[$i][$j]); |
} |
if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) { |
$value = trim($value,'"'); |
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($value); |
} |
if ($validValues) { |
$this->A[$i][$j] = pow($this->A[$i][$j],$value); |
} else { |
$this->A[$i][$j] = PHPExcel_Calculation_Functions::NaN(); |
} |
} |
} |
return $this; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function power() |
/** |
* concat |
* |
* A = A & B |
* @param mixed $B Matrix/Array |
* @return Matrix Sum |
*/ |
public function concat() { |
if (func_num_args() > 0) { |
$args = func_get_args(); |
$match = implode(",", array_map('gettype', $args)); |
switch($match) { |
case 'object': |
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); } |
case 'array': |
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]); |
break; |
default: |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
break; |
} |
$this->checkMatrixDimensions($M); |
for($i = 0; $i < $this->m; ++$i) { |
for($j = 0; $j < $this->n; ++$j) { |
$this->A[$i][$j] = trim($this->A[$i][$j],'"').trim($M->get($i, $j),'"'); |
} |
} |
return $this; |
} else { |
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException); |
} |
} // function concat() |
/** |
* Solve A*X = B. |
* |
* @param Matrix $B Right hand side |
* @return Matrix ... Solution if A is square, least squares solution otherwise |
*/ |
public function solve($B) { |
if ($this->m == $this->n) { |
$LU = new PHPExcel_Shared_JAMA_LUDecomposition($this); |
return $LU->solve($B); |
} else { |
$QR = new QRDecomposition($this); |
return $QR->solve($B); |
} |
} // function solve() |
/** |
* Matrix inverse or pseudoinverse. |
* |
* @return Matrix ... Inverse(A) if A is square, pseudoinverse otherwise. |
*/ |
public function inverse() { |
return $this->solve($this->identity($this->m, $this->m)); |
} // function inverse() |
/** |
* det |
* |
* Calculate determinant |
* @return float Determinant |
*/ |
public function det() { |
$L = new PHPExcel_Shared_JAMA_LUDecomposition($this); |
return $L->det(); |
} // function det() |
} // class PHPExcel_Shared_JAMA_Matrix |
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/utils/Error.php |
---|
New file |
0,0 → 1,82 |
<?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()"); |
} |
} |
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/utils/Maths.php |
---|
New file |
0,0 → 1,43 |
<?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 PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException)); |
} |
} |
return sqrt($s); |
} |
*/ |
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/LUDecomposition.php |
---|
New file |
0,0 → 1,258 |
<?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 PHPExcel_Calculation_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 PHPExcel_Calculation_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,:) |
* @PHPExcel_Calculation_Exception IllegalArgumentException Matrix row dimensions must agree. |
* @PHPExcel_Calculation_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 PHPExcel_Calculation_Exception(self::MatrixSingularException); |
} |
} else { |
throw new PHPExcel_Calculation_Exception(self::MatrixSquareException); |
} |
} // function solve() |
} // class PHPExcel_Shared_JAMA_LUDecomposition |
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/CholeskyDecomposition.php |
---|
New file |
0,0 → 1,149 |
<?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 PHPExcel_Calculation_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 PHPExcel_Calculation_Exception(JAMAError(MatrixSPDException)); |
} |
} else { |
throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionException)); |
} |
} else { |
throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException)); |
} |
} // function solve() |
} // class CholeskyDecomposition |
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/SingularValueDecomposition.php |
---|
New file |
0,0 → 1,526 |
<?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 |
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/CHANGELOG.TXT |
---|
New file |
0,0 → 1,16 |
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. |
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/QRDecomposition.php |
---|
New file |
0,0 → 1,234 |
<?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 PHPExcel_Calculation_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 PHPExcel_Calculation_Exception(self::MatrixRankException); |
} |
} else { |
throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException); |
} |
} // function solve() |
} // PHPExcel_Shared_JAMA_class QRDecomposition |