diff --git a/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/CHANGELOG.TXT b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/CHANGELOG.TXT new file mode 100644 index 0000000..1c18a5d --- /dev/null +++ b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/CHANGELOG.TXT @@ -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. + diff --git a/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/CholeskyDecomposition.php b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/CholeskyDecomposition.php new file mode 100644 index 0000000..cfbaa53 --- /dev/null +++ b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/CholeskyDecomposition.php @@ -0,0 +1,149 @@ +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 diff --git a/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php new file mode 100644 index 0000000..2a696d0 --- /dev/null +++ b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php @@ -0,0 +1,862 @@ +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 diff --git a/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/LUDecomposition.php b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/LUDecomposition.php new file mode 100644 index 0000000..08e500c --- /dev/null +++ b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/LUDecomposition.php @@ -0,0 +1,258 @@ += 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 diff --git a/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/Matrix.php b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/Matrix.php new file mode 100644 index 0000000..1e7c334 --- /dev/null +++ b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/Matrix.php @@ -0,0 +1,1057 @@ + 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 PHPExcel_Shared_JAMA_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 diff --git a/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/QRDecomposition.php b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/QRDecomposition.php new file mode 100644 index 0000000..7538462 --- /dev/null +++ b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/QRDecomposition.php @@ -0,0 +1,234 @@ += 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 diff --git a/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/SingularValueDecomposition.php b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/SingularValueDecomposition.php new file mode 100644 index 0000000..a4b096c --- /dev/null +++ b/includes/PHPExcel/Classes/PHPExcel/Shared/JAMA/SingularValueDecomposition.php @@ -0,0 +1,526 @@ += 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

= -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