diff --git a/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer.php b/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer.php index 588cae727..1084d94f6 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer.php +++ b/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer.php @@ -22,7 +22,7 @@ * @package PHPExcel_Shared_Escher * @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel) * @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL - * @version 1.6.7, 2009-04-22 + * @version 1.7.0, 2009-08-10 */ /** diff --git a/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer/SpgrContainer.php b/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer/SpgrContainer.php index 3114f9c62..846b56990 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer/SpgrContainer.php +++ b/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer/SpgrContainer.php @@ -22,7 +22,7 @@ * @package PHPExcel_Shared_Escher * @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel) * @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL - * @version 1.6.7, 2009-04-22 + * @version 1.7.0, 2009-08-10 */ /** diff --git a/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer/SpgrContainer/SpContainer.php b/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer/SpgrContainer/SpContainer.php index e05188b74..481e7db55 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer/SpgrContainer/SpContainer.php +++ b/libraries/PHPExcel/PHPExcel/Shared/Escher/DgContainer/SpgrContainer/SpContainer.php @@ -22,7 +22,7 @@ * @package PHPExcel_Shared_Escher * @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel) * @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL - * @version 1.6.7, 2009-04-22 + * @version 1.7.0, 2009-08-10 */ /** diff --git a/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer.php b/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer.php index 7b1741c82..cdbdde01f 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer.php +++ b/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer.php @@ -22,7 +22,7 @@ * @package PHPExcel_Shared_Escher * @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel) * @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL - * @version 1.6.7, 2009-04-22 + * @version 1.7.0, 2009-08-10 */ /** diff --git a/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer.php b/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer.php index 4c26bb384..4209b424c 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer.php +++ b/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer.php @@ -22,7 +22,7 @@ * @package PHPExcel_Shared_Escher * @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel) * @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL - * @version 1.6.7, 2009-04-22 + * @version 1.7.0, 2009-08-10 */ /** diff --git a/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer/BSE.php b/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer/BSE.php index ae0f246d2..e472b218c 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer/BSE.php +++ b/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer/BSE.php @@ -22,7 +22,7 @@ * @package PHPExcel_Shared_Escher * @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel) * @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL - * @version 1.6.7, 2009-04-22 + * @version 1.7.0, 2009-08-10 */ /** diff --git a/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer/BSE/Blip.php b/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer/BSE/Blip.php index 41d095657..d835a0687 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer/BSE/Blip.php +++ b/libraries/PHPExcel/PHPExcel/Shared/Escher/DggContainer/BstoreContainer/BSE/Blip.php @@ -22,7 +22,7 @@ * @package PHPExcel_Shared_Escher * @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel) * @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL - * @version 1.6.7, 2009-04-22 + * @version 1.7.0, 2009-08-10 */ /** diff --git a/libraries/PHPExcel/PHPExcel/Shared/JAMA/CholeskyDecomposition.php b/libraries/PHPExcel/PHPExcel/Shared/JAMA/CholeskyDecomposition.php index b401b8491..9d064f9e6 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/JAMA/CholeskyDecomposition.php +++ b/libraries/PHPExcel/PHPExcel/Shared/JAMA/CholeskyDecomposition.php @@ -1,133 +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 { - trigger_error(ArgumentTypeException, ERROR); - } - } - - /** - * Is the matrix symmetric and positive definite? - * @return boolean - */ - function isSPD () { - return $this->isspd; - } - - /** - * getL - * Return triangular factor. - * @return Matrix Lower triangular matrix - */ - function getL () { - return new Matrix($this->L); - } - - /** - * Solve A*X = B - * @param $B Row-equal matrix - * @return Matrix L * L' * X = B - */ - function solve ( $B = null ) { - if( is_a($B, '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 { - trigger_error(MatrixSPDException, ERROR); - } - } else { - trigger_error(MatrixDimensionException, ERROR); - } - } else { - trigger_error(ArgumentTypeException, ERROR); - } - } -} + /** + * Decomposition storage + * @var array + * @access private + */ + private $L = array(); + + /** + * Matrix row and column dimension + * @var int + * @access private + */ + private $m; + + /** + * Symmetric positive definite flag + * @var boolean + * @access private + */ + private $isspd = true; + + + /** + * CholeskyDecomposition + * + * Class constructor - decomposes symmetric positive definite matrix + * @param mixed Matrix square symmetric positive definite matrix + */ + public function __construct($A = null) { + if ($A instanceof Matrix) { + $this->L = $A->getArray(); + $this->m = $A->getRowDimension(); + + for($i = 0; $i < $this->m; ++$i) { + for($j = $i; $j < $this->m; ++$j) { + for($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) { + $sum -= $this->L[$i][$k] * $this->L[$j][$k]; + } + if ($i == $j) { + if ($sum >= 0) { + $this->L[$i][$i] = sqrt($sum); + } else { + $this->isspd = false; + } + } else { + if ($this->L[$i][$i] != 0) { + $this->L[$j][$i] = $sum / $this->L[$i][$i]; + } + } + } + + for ($k = $i+1; $k < $this->m; ++$k) { + $this->L[$i][$k] = 0.0; + } + } + } else { + throw new Exception(JAMAError(ArgumentTypeException)); + } + } // function __construct() + + + /** + * Is the matrix symmetric and positive definite? + * + * @return boolean + */ + public function isSPD() { + return $this->isspd; + } // function isSPD() + + + /** + * getL + * + * Return triangular factor. + * @return Matrix Lower triangular matrix + */ + public function getL() { + return new Matrix($this->L); + } // function getL() + + + /** + * Solve A*X = B + * + * @param $B Row-equal matrix + * @return Matrix L * L' * X = B + */ + public function solve($B = null) { + if ($B instanceof Matrix) { + if ($B->getRowDimension() == $this->m) { + if ($this->isspd) { + $X = $B->getArrayCopy(); + $nx = $B->getColumnDimension(); + + for ($k = 0; $k < $this->m; ++$k) { + for ($i = $k + 1; $i < $this->m; ++$i) { + for ($j = 0; $j < $nx; ++$j) { + $X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k]; + } + } + for ($j = 0; $j < $nx; ++$j) { + $X[$k][$j] /= $this->L[$k][$k]; + } + } + + for ($k = $this->m - 1; $k >= 0; --$k) { + for ($j = 0; $j < $nx; ++$j) { + $X[$k][$j] /= $this->L[$k][$k]; + } + for ($i = 0; $i < $k; ++$i) { + for ($j = 0; $j < $nx; ++$j) { + $X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i]; + } + } + } + + return new Matrix($X, $this->m, $nx); + } else { + throw new Exception(JAMAError(MatrixSPDException)); + } + } else { + throw new Exception(JAMAError(MatrixDimensionException)); + } + } else { + throw new Exception(JAMAError(ArgumentTypeException)); + } + } // function solve() + +} // class CholeskyDecomposition diff --git a/libraries/PHPExcel/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php b/libraries/PHPExcel/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php index a010884b7..2a696d00f 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php +++ b/libraries/PHPExcel/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php @@ -1,789 +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; - } + /** + * 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; + } - $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; - } + // 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; + } + } - /** - * 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 - */ - 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; - } - } - } - } + $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; + } - /** - * 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 - */ - 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 - */ - 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 - */ - 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; - else if ($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. - */ - function EigenvalueDecomposition($Arg) { - $this->A = $Arg->getArray(); - $this->n = $Arg->getColumnDimension(); - $this->V = array(); - $this->d = array(); - $this->e = array(); - $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(); - } - } + /** + * 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); - /** - * Return the eigenvector matrix - * @access public - * @return V - */ - function getV() { - return new Matrix($this->V, $this->n, $this->n); - } + 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; + } - /** - * Return the real parts of the eigenvalues - * @access public - * @return real(diag(D)) - */ - function getRealEigenvalues() { - return $this->d; - } + // 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; + } + } + } + } - /** - * Return the imaginary parts of the eigenvalues - * @access public - * @return imag(diag(D)) - */ - function getImagEigenvalues() { - return $this->e; - } - /** - * Return the block diagonal eigenvalue matrix - * @access public - * @return D - */ - 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); - } -} + /** + * 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/libraries/PHPExcel/PHPExcel/Shared/JAMA/LUDecomposition.php b/libraries/PHPExcel/PHPExcel/Shared/JAMA/LUDecomposition.php index de177a5f6..4fd43f91f 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/JAMA/LUDecomposition.php +++ b/libraries/PHPExcel/PHPExcel/Shared/JAMA/LUDecomposition.php @@ -1,222 +1,255 @@ = 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 -*/ + * @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 LUDecomposition { - /** - * Decomposition storage - * @var array - */ - var $LU = array(); - - /** - * Row dimension. - * @var int - */ - var $m; - /** - * Column dimension. - * @var int - */ - var $n; - - /** - * Pivot sign. - * @var int - */ - var $pivsign; + /** + * Decomposition storage + * @var array + */ + private $LU = array(); - /** - * Internal storage of pivot vector. - * @var array - */ - var $piv = array(); - - /** - * LU Decomposition constructor. - * @param $A Rectangular matrix - * @return Structure to access L, U and piv. - */ - function LUDecomposition ($A) { - if( is_a($A, 'Matrix') ) { - // Use a "left-looking", dot-product, Crout/Doolittle algorithm. - $this->LU = $A->getArrayCopy(); - $this->m = $A->getRowDimension(); - $this->n = $A->getColumnDimension(); - for ($i = 0; $i < $this->m; $i++) - $this->piv[$i] = $i; - $this->pivsign = 1; - $LUrowi = array(); - $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) AND ($this->LU[$j][$j] != 0.0) ) { - for ($i = $j+1; $i < $this->m; $i++) - $this->LU[$i][$j] /= $this->LU[$j][$j]; - } - } - } else { - trigger_error(ArgumentTypeException, ERROR); - } - } - - /** - * Get lower triangular factor. - * @return array Lower triangular factor - */ - 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]; - else if($i == $j) - $L[$i][$j] = 1.0; - else - $L[$i][$j] = 0.0; - } - } - return new Matrix($L); - } + /** + * Row dimension. + * @var int + */ + private $m; - /** - * Get upper triangular factor. - * @return array Upper triangular factor - */ - 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 Matrix($U); - } - - /** - * Return pivot permutation vector. - * @return array Pivot vector - */ - function getPivot () { - return $this->piv; - } - - /** - * Alias for getPivot - * @see getPivot - */ - function getDoublePivot () { - return $this->getPivot(); - } + /** + * Column dimension. + * @var int + */ + private $n; - /** - * Is the matrix nonsingular? - * @return true if U, and hence A, is nonsingular. - */ - function isNonsingular () { - for ($j = 0; $j < $this->n; $j++) { - if ($this->LU[$j][$j] == 0) - return false; - } - return true; - } + /** + * Pivot sign. + * @var int + */ + private $pivsign; - /** - * Count determinants - * @return array d matrix deterninat - */ - 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 { - trigger_error(MatrixDimensionException, ERROR); - } - } - - /** - * Solve A*X = B - * @param $B A Matrix with as many rows as A and any number of columns. - * @return X so that L*U*X = B(piv,:) - * @exception IllegalArgumentException Matrix row dimensions must agree. - * @exception RuntimeException Matrix is singular. - */ - 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 { - trigger_error(MatrixSingularException, ERROR); - } - } else { - trigger_error(MatrixSquareException, ERROR); - } - } -} + /** + * 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 Matrix) { + // Use a "left-looking", dot-product, Crout/Doolittle algorithm. + $this->LU = $A->getArrayCopy(); + $this->m = $A->getRowDimension(); + $this->n = $A->getColumnDimension(); + for ($i = 0; $i < $this->m; ++$i) { + $this->piv[$i] = $i; + } + $this->pivsign = 1; + $LUrowi = $LUcolj = array(); + + // Outer loop. + for ($j = 0; $j < $this->n; ++$j) { + // Make a copy of the j-th column to localize references. + for ($i = 0; $i < $this->m; ++$i) { + $LUcolj[$i] = &$this->LU[$i][$j]; + } + // Apply previous transformations. + for ($i = 0; $i < $this->m; ++$i) { + $LUrowi = $this->LU[$i]; + // Most of the time is spent in the following dot product. + $kmax = min($i,$j); + $s = 0.0; + for ($k = 0; $k < $kmax; ++$k) { + $s += $LUrowi[$k] * $LUcolj[$k]; + } + $LUrowi[$j] = $LUcolj[$i] -= $s; + } + // Find pivot and exchange if necessary. + $p = $j; + for ($i = $j+1; $i < $this->m; ++$i) { + if (abs($LUcolj[$i]) > abs($LUcolj[$p])) { + $p = $i; + } + } + if ($p != $j) { + for ($k = 0; $k < $this->n; ++$k) { + $t = $this->LU[$p][$k]; + $this->LU[$p][$k] = $this->LU[$j][$k]; + $this->LU[$j][$k] = $t; + } + $k = $this->piv[$p]; + $this->piv[$p] = $this->piv[$j]; + $this->piv[$j] = $k; + $this->pivsign = $this->pivsign * -1; + } + // Compute multipliers. + if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) { + for ($i = $j+1; $i < $this->m; ++$i) { + $this->LU[$i][$j] /= $this->LU[$j][$j]; + } + } + } + } else { + throw new Exception(JAMAError(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 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 Matrix($U); + } // function getU() + + + /** + * Return pivot permutation vector. + * + * @return array Pivot vector + */ + public function getPivot() { + return $this->piv; + } // function getPivot() + + + /** + * Alias for getPivot + * + * @see getPivot + */ + public function getDoublePivot() { + return $this->getPivot(); + } // function getDoublePivot() + + + /** + * Is the matrix nonsingular? + * + * @return true if U, and hence A, is nonsingular. + */ + public function isNonsingular() { + for ($j = 0; $j < $this->n; ++$j) { + if ($this->LU[$j][$j] == 0) { + return false; + } + } + return true; + } // function isNonsingular() + + + /** + * Count determinants + * + * @return array d matrix deterninat + */ + public function det() { + if ($this->m == $this->n) { + $d = $this->pivsign; + for ($j = 0; $j < $this->n; ++$j) { + $d *= $this->LU[$j][$j]; + } + return $d; + } else { + throw new Exception(JAMAError(MatrixDimensionException)); + } + } // function det() + + + /** + * Solve A*X = B + * + * @param $B A Matrix with as many rows as A and any number of columns. + * @return X so that L*U*X = B(piv,:) + * @exception IllegalArgumentException Matrix row dimensions must agree. + * @exception RuntimeException Matrix is singular. + */ + public function solve($B) { + if ($B->getRowDimension() == $this->m) { + if ($this->isNonsingular()) { + // Copy right hand side with pivoting + $nx = $B->getColumnDimension(); + $X = $B->getMatrix($this->piv, 0, $nx-1); + // Solve L*Y = B(piv,:) + for ($k = 0; $k < $this->n; ++$k) { + for ($i = $k+1; $i < $this->n; ++$i) { + for ($j = 0; $j < $nx; ++$j) { + $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k]; + } + } + } + // Solve U*X = Y; + for ($k = $this->n-1; $k >= 0; --$k) { + for ($j = 0; $j < $nx; ++$j) { + $X->A[$k][$j] /= $this->LU[$k][$k]; + } + for ($i = 0; $i < $k; ++$i) { + for ($j = 0; $j < $nx; ++$j) { + $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k]; + } + } + } + return $X; + } else { + throw new Exception(JAMAError(MatrixSingularException)); + } + } else { + throw new Exception(JAMAError(MatrixSquareException)); + } + } // function solve() + +} // class LUDecomposition diff --git a/libraries/PHPExcel/PHPExcel/Shared/JAMA/Matrix.php b/libraries/PHPExcel/PHPExcel/Shared/JAMA/Matrix.php index c1b79ae04..efababbc4 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/JAMA/Matrix.php +++ b/libraries/PHPExcel/PHPExcel/Shared/JAMA/Matrix.php @@ -1,1331 +1,1361 @@ 0 ) { - - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - - //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 constant-filled - m x n filled with c - case 'integer,integer,integer': - $this->m = $args[0]; - $this->n = $args[1]; - $this->A = array_fill(0, $this->m, array_fill(0, $this->n, $args[2])); - break; - - //Rectangular matrix constant-filled - m x n filled with c - case 'integer,integer,double': - $this->m = $args[0]; - $this->n = $args[1]; - $this->A = array_fill(0, $this->m, array_fill(0, $this->n, $args[2])); - break; - - //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; - - //Rectangular matrix - m x n initialized from 2D array - case 'array,integer,integer': - $this->m = $args[1]; - $this->n = $args[2]; - $this->A = $args[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 - trigger_error(ArrayLengthException, ERROR); - - break; - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else - trigger_error(PolymorphicArgumentException, ERROR); - } - - /** - * getArray - * @return array Matrix array - */ - function &getArray() { - return $this->A; - } - - /** - * getArrayCopy - * @return array Matrix array copy - */ - function getArrayCopy() { - return $this->A; - } - - /** Construct a matrix from a copy of a 2-D array. - * @param double A[][] Two-dimensional array of doubles. - * @exception IllegalArgumentException All rows must have the same length - */ - function constructWithCopy($A) { - $this->m = count($A); - $this->n = count($A[0]); - $X = new Matrix($this->m, $this->n); - for ($i = 0; $i < $this->m; $i++) { - if (count($A[$i]) != $this->n) - trigger_error(RowLengthException, ERROR); - for ($j = 0; $j < $this->n; $j++) - $X->A[$i][$j] = $A[$i][$j]; - } - return $X; - } - - /** - * getColumnPacked - * Get a column-packed array - * @return array Column-packed matrix array - */ - function getColumnPackedCopy() { - $P = array(); - for($i = 0; $i < $this->m; $i++) { - for($j = 0; $j < $this->n; $j++) { - array_push($P, $this->A[$j][$i]); - } - } - return $P; - } - - /** - * getRowPacked - * Get a row-packed array - * @return array Row-packed matrix array - */ - function getRowPackedCopy() { - $P = array(); - for($i = 0; $i < $this->m; $i++) { - for($j = 0; $j < $this->n; $j++) { - array_push($P, $this->A[$i][$j]); - } - } - return $P; - } - - /** - * getRowDimension - * @return int Row dimension - */ - function getRowDimension() { - return $this->m; - } - - /** - * getColumnDimension - * @return int Column dimension - */ - function getColumnDimension() { - return $this->n; - } - - /** - * 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) - */ - function get( $i = null, $j = null ) { - return $this->A[$i][$j]; - } - - /** - * 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 - */ - 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; - $m = $i0 >= 0 ? $this->m - $i0 : trigger_error(ArgumentBoundsException, ERROR); - $n = $j0 >= 0 ? $this->n - $j0 : trigger_error(ArgumentBoundsException, ERROR); - $R = new 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; - $m = ( ($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0) ) ? $iF - $i0 : trigger_error(ArgumentBoundsException, ERROR); - $n = ( ($jF > $j0) && ($this->n >= $jF) && ($j0 >= 0) ) ? $jF - $j0 : trigger_error(ArgumentBoundsException, ERROR); - $R = new 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; - $m = count($RL) > 0 ? count($RL) : trigger_error(ArgumentBoundsException, ERROR); - $n = count($CL) > 0 ? count($CL) : trigger_error(ArgumentBoundsException, ERROR); - $R = new 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; - $m = count($RL) > 0 ? count($RL) : trigger_error(ArgumentBoundsException, ERROR); - $n = count($CL) > 0 ? count($CL) : trigger_error(ArgumentBoundsException, ERROR); - $R = new 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; - $m = ( ($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0) ) ? $iF - $i0 : trigger_error(ArgumentBoundsException, ERROR); - $n = count($CL) > 0 ? count($CL) : trigger_error(ArgumentBoundsException, ERROR); - $R = new 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; - $m = count($RL) > 0 ? count($RL) : trigger_error(ArgumentBoundsException, ERROR); - $n = ( ($jF >= $j0) && ($this->n >= $jF) && ($j0 >= 0) ) ? $jF - $j0 : trigger_error(ArgumentBoundsException, ERROR); - $R = new 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: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * setMatrix - * Set a submatrix - * @param int $i0 Initial row index - * @param int $j0 Initial column index - * @param mixed $S Matrix/Array submatrix - * ($i0, $j0, $S) $S = Matrix - * ($i0, $j0, $S) $S = Array - */ - function setMatrix( ) { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'integer,integer,object': - $M = is_a($args[2], 'Matrix') ? $args[2] : trigger_error(ArgumentTypeException, ERROR); - $i0 = ( ($args[0] + $M->m) <= $this->m ) ? $args[0] : trigger_error(ArgumentBoundsException, ERROR); - $j0 = ( ($args[1] + $M->n) <= $this->n ) ? $args[1] : trigger_error(ArgumentBoundsException, ERROR); - - for($i = $i0; $i < $i0 + $M->m; $i++) { - for($j = $j0; $j < $j0 + $M->n; $j++) { - $this->A[$i][$j] = $M->get($i - $i0, $j - $j0); - } - } - - break; - - case 'integer,integer,array': - $M = new Matrix($args[2]); - $i0 = ( ($args[0] + $M->m) <= $this->m ) ? $args[0] : trigger_error(ArgumentBoundsException, ERROR); - $j0 = ( ($args[1] + $M->n) <= $this->n ) ? $args[1] : trigger_error(ArgumentBoundsException, ERROR); - - for($i = $i0; $i < $i0 + $M->m; $i++) { - for($j = $j0; $j < $j0 + $M->n; $j++) { - $this->A[$i][$j] = $M->get($i - $i0, $j - $j0); - } - } - - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * checkMatrixDimensions - * Is matrix B the same size? - * @param Matrix $B Matrix B - * @return boolean - */ - function checkMatrixDimensions( $B = null ) { - if( is_a($B, 'Matrix') ) - if( ($this->m == $B->m) && ($this->n == $B->n) ) - return true; - else - trigger_error(MatrixDimensionException, ERROR); - - else - trigger_error(ArgumentTypeException, ERROR); - } - - - - /** - * 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) - */ - function set( $i = null, $j = null, $c = null ) { - // Optimized set version just has this - $this->A[$i][$j] = $c; - /* - if( is_int($i) && is_int($j) && is_numeric($c) ) { - if( ( $i < $this->m ) && ( $j < $this->n ) ) { - $this->A[$i][$j] = $c; - } else { - echo "A[$i][$j] = $c
"; - trigger_error(ArgumentBoundsException, WARNING); - } - } else { - trigger_error(ArgumentTypeException, WARNING); - } - */ - } - - /** - * identity - * Generate an identity matrix. - * @param int $m Row dimension - * @param int $n Column dimension - * @return Matrix Identity matrix - */ - function &identity( $m = null, $n = null ) { - return Matrix::diagonal($m, $n, 1); - } - - /** - * diagonal - * Generate a diagonal matrix - * @param int $m Row dimension - * @param int $n Column dimension - * @param mixed $c Diagonal value - * @return Matrix Diagonal matrix - */ - function &diagonal( $m = null, $n = null, $c = 1 ) { - $R = new Matrix($m, $n); - for($i = 0; $i < $m; $i++) - $R->set($i, $i, $c); - return $R; - } - - /** - * filled - * Generate a filled matrix - * @param int $m Row dimension - * @param int $n Column dimension - * @param int $c Fill constant - * @return Matrix Filled matrix - */ - function &filled( $m = null, $n = null, $c = 0 ) { - if( is_int($m) && is_int($n) && is_numeric($c) ) { - $R = new Matrix($m, $n, $c); - return $R; - } else { - trigger_error(ArgumentTypeException, ERROR); - } - } - - /** - * random - * Generate a random matrix - * @param int $m Row dimension - * @param int $n Column dimension - * @return Matrix Random matrix - */ - function &random( $m = null, $n = null, $a = RAND_MIN, $b = RAND_MAX ) { - if( is_int($m) && is_int($n) && is_numeric($a) && is_numeric($b) ) { - $R = new Matrix($m, $n); - - for($i = 0; $i < $m; $i++) - for($j = 0; $j < $n; $j++) - $R->set($i, $j, mt_rand($a, $b)); - - return $R; - } else { - trigger_error(ArgumentTypeException, ERROR); - } - } - - /** - * packed - * Alias for getRowPacked - * @return array Packed array - */ - function &packed() { - return $this->getRowPacked(); - } - - - /** - * getMatrixByRow - * Get a submatrix by row index/range - * @param int $i0 Initial row index - * @param int $iF Final row index - * @return Matrix Submatrix - */ - 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 - trigger_error(ArgumentTypeException, ERROR); - } - - /** - * getMatrixByCol - * Get a submatrix by column index/range - * @param int $i0 Initial column index - * @param int $iF Final column index - * @return Matrix Submatrix - */ - 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 - trigger_error(ArgumentTypeException, ERROR); - } - - /** - * transpose - * Tranpose matrix - * @return Matrix Transposed matrix - */ - function transpose() { - $R = new 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; - } - -/* - public Matrix transpose () { - Matrix X = new Matrix(n,m); - double[][] C = X.getArray(); - for (int i = 0; i < m; i++) { - for (int j = 0; j < n; j++) { - C[j][i] = A[i][j]; - } - } - return X; - } -*/ - - /** - * norm1 - * One norm - * @return float Maximum column sum - */ - function norm1() { - $r = 0; - - for($j = 0; $j < $this->n; $j++) { - $s = 0; - - for($i = 0; $i < $this->m; $i++) { - $s += abs($this->A[$i][$j]); - } - - $r = ( $r > $s ) ? $r : $s; - } - - return $r; - } - - - /** - * norm2 - * Maximum singular value - * @return float Maximum singular value - */ - function norm2() { - - } - - /** - * normInf - * Infinite norm - * @return float Maximum row sum - */ - function normInf() { - $r = 0; - - for($i = 0; $i < $this->m; $i++) { - $s = 0; - - for($j = 0; $j < $this->n; $j++) { - $s += abs($this->A[$i][$j]); - } - - $r = ( $r > $s ) ? $r : $s; - } - - return $r; - } - - /** - * normF - * Frobenius norm - * @return float Square root of the sum of all elements squared - */ - function normF() { - $f = 0; - for ($i = 0; $i < $this->m; $i++) - for ($j = 0; $j < $this->n; $j++) - $f = hypo($f,$this->A[$i][$j]); - return $f; - } - - /** - * Matrix rank - * @return effective numerical rank, obtained from SVD. - */ - function rank () { - $svd = new SingularValueDecomposition($this); - return $svd->rank(); - } - - /** - * Matrix condition (2 norm) - * @return ratio of largest to smallest singular value. - */ - function cond () { - $svd = new SingularValueDecomposition($this); - return $svd->cond(); - } - - /** - * trace - * Sum of diagonal elements - * @return float Sum of diagonal elements - */ - function trace() { - $s = 0; - $n = min($this->m, $this->n); - - for($i = 0; $i < $n; $i++) - $s += $this->A[$i][$i]; - - return $s; - } - - - /** - * uminus - * Unary minus matrix -A - * @return Matrix Unary minus matrix - */ - function uminus() { - - } - - /** - * plus - * A + B - * @param mixed $B Matrix/Array - * @return Matrix Sum - */ - function plus() { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'object': - $M = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - //$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; - break; - - case 'array': - $M = new Matrix($args[0]); - //$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; - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * plusEquals - * A = A + B - * @param mixed $B Matrix/Array - * @return Matrix Sum - */ - function &plusEquals() { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'object': - $M = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - $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); - } - } - - return $this; - break; - - case 'array': - $M = new Matrix($args[0]); - $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); - } - } - - return $this; - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * minus - * A - B - * @param mixed $B Matrix/Array - * @return Matrix Sum - */ - function minus() { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'object': - $M = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - $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; - break; - - case 'array': - $M = new Matrix($args[0]); - $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; - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * minusEquals - * A = A - B - * @param mixed $B Matrix/Array - * @return Matrix Sum - */ - function &minusEquals() { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'object': - $M = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - $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); - } - } - - return $this; - break; - - case 'array': - $M = new Matrix($args[0]); - $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); - } - } - - return $this; - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * arrayTimes - * Element-by-element multiplication - * Cij = Aij * Bij - * @param mixed $B Matrix/Array - * @return Matrix Matrix Cij - */ - function arrayTimes() { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'object': - $M = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - $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; - break; - - case 'array': - $M = new Matrix($args[0]); - $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; - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - - /** - * arrayTimesEquals - * Element-by-element multiplication - * Aij = Aij * Bij - * @param mixed $B Matrix/Array - * @return Matrix Matrix Aij - */ - function &arrayTimesEquals() { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'object': - $M = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - $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); - } - } - - return $this; - break; - - case 'array': - $M = new Matrix($args[0]); - $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); - } - } - - return $this; - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * arrayRightDivide - * Element-by-element right division - * A / B - * @param Matrix $B Matrix B - * @return Matrix Division result - */ - function arrayRightDivide() { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'object': - $M = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - $this->checkMatrixDimensions($M); - - for($i = 0; $i < $this->m; $i++) { - for($j = 0; $j < $this->n; $j++) { - $M->set($i, $j, $this->A[$i][$j] / $M->get($i, $j) ); - } - } - - return $M; - break; - - case 'array': - $M = new Matrix($args[0]); - $this->checkMatrixDimensions($M); - - for($i = 0; $i < $this->m; $i++) { - for($j = 0; $j < $this->n; $j++) { - $M->set($i, $j, $this->A[$i][$j] / $M->get($i, $j)); - } - } - - return $M; - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * arrayRightDivideEquals - * Element-by-element right division - * Aij = Aij / Bij - * @param mixed $B Matrix/Array - * @return Matrix Matrix Aij - */ - function &arrayRightDivideEquals() { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'object': - $M = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - $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; - break; - - case 'array': - $M = new Matrix($args[0]); - $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; - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * arrayLeftDivide - * Element-by-element Left division - * A / B - * @param Matrix $B Matrix B - * @return Matrix Division result - */ - function arrayLeftDivide() { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'object': - $M = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - $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; - break; - - case 'array': - $M = new Matrix($args[0]); - $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; - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * arrayLeftDivideEquals - * Element-by-element Left division - * Aij = Aij / Bij - * @param mixed $B Matrix/Array - * @return Matrix Matrix Aij - */ - function &arrayLeftDivideEquals() { - if( func_num_args() > 0 ) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - - switch( $match ) { - case 'object': - $M = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - $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; - break; - - case 'array': - $M = new Matrix($args[0]); - $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; - break; - - default: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else { - trigger_error(PolymorphicArgumentException, ERROR); - } - } - - /** - * times - * Matrix multiplication - * @param mixed $n Matrix/Array/Scalar - * @return Matrix Product - */ - function times() { - if(func_num_args() > 0) { - $args = func_get_args(); - $match = implode(",", array_map('gettype', $args)); - switch($match) { - case 'object': - $B = is_a($args[0], 'Matrix') ? $args[0] : trigger_error(ArgumentTypeException, ERROR); - if($this->n == $B->m) { - $C = new 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 - trigger_error(MatrixDimensionMismatch, FATAL); - break; - - case 'array': - $B = new Matrix($args[0]); - if($this->n == $B->m) { - $C = new 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 - trigger_error(MatrixDimensionMismatch, FATAL); - return $M; - break; - case 'integer': - $C = new 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 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 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: - trigger_error(PolymorphicArgumentException, ERROR); - break; - } - } else - trigger_error(PolymorphicArgumentException, ERROR); - } - - /** - * chol - * Cholesky decomposition - * @return Matrix Cholesky decomposition - */ - function chol() { - return new CholeskyDecomposition($this); - } - - /** - * lu - * LU decomposition - * @return Matrix LU decomposition - */ - function lu() { - return new LUDecomposition($this); - } - - /** - * qr - * QR decomposition - * @return Matrix QR decomposition - */ - function qr() { - return new QRDecomposition($this); - } - - - /** - * eig - * Eigenvalue decomposition - * @return Matrix Eigenvalue decomposition - */ - function eig() { - return new EigenvalueDecomposition($this); - } - - /** - * svd - * Singular value decomposition - * @return Singular value decomposition - */ - function svd() { - return new SingularValueDecomposition($this); - } - - /** - * Solve A*X = B. - * @param Matrix $B Right hand side - * @return Matrix ... Solution if A is square, least squares solution otherwise - */ - function solve($B) { - if ($this->m == $this->n) { - $LU = new LUDecomposition($this); - return $LU->solve($B); - } else { - $QR = new QRDecomposition($this); - return $QR->solve($B); - } - } - - /** - * Matrix inverse or pseudoinverse. - * @return Matrix ... Inverse(A) if A is square, pseudoinverse otherwise. - */ - function inverse() { - return $this->solve($this->identity($this->m, $this->m)); - } - - - /** - * det - * Calculate determinant - * @return float Determinant - */ - function det() { - $L = new LUDecomposition($this); - return $L->det(); - } - - /** - * Older debugging utility for backwards compatability. - * @return html version of matrix - */ - function mprint($A, $format="%01.2f", $width=2) { - $spacing = ""; - $m = count($A); - $n = count($A[0]); - for($i = 0; $i < $width; $i++) - $spacing .= " "; - for ($i = 0; $i < $m; $i++) { - for ($j = 0; $j < $n; $j++) { - $formatted = sprintf($format, $A[$i][$j]); - echo $formatted . $spacing; - } - echo "
"; - } - } - - /** - * Debugging utility. - * @return Output HTML representation of matrix - */ - function toHTML($width=2) { - print( ''); - for( $i = 0; $i < $this->m; $i++ ) { - print( '' ); - for( $j = 0; $j < $this->n; $j++ ) - print( '' ); - print( ''); - } - print( '
' . $this->A[$i][$j] . '
' ); - } +/** PHPExcel root directory */ +if (!defined('PHPEXCEL_ROOT')) { + /** + * @ignore + */ + define('PHPEXCEL_ROOT', dirname(__FILE__) . '/../../../'); } + +require_once PHPEXCEL_ROOT . 'PHPExcel/Shared/JAMA/utils/Error.php'; +require_once PHPEXCEL_ROOT . 'PHPExcel/Shared/JAMA/utils/Maths.php'; +require_once PHPEXCEL_ROOT . 'PHPExcel/Shared/JAMA/CholeskyDecomposition.php'; +require_once PHPEXCEL_ROOT . 'PHPExcel/Shared/JAMA/LUDecomposition.php'; +require_once PHPEXCEL_ROOT . 'PHPExcel/Shared/JAMA/QRDecomposition.php'; +require_once PHPEXCEL_ROOT . 'PHPExcel/Shared/JAMA/EigenvalueDecomposition.php'; +require_once PHPEXCEL_ROOT . 'PHPExcel/Shared/JAMA/SingularValueDecomposition.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 Matrix { + + /** + * 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) { + //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 constant-filled - m x n filled with c + case 'integer,integer,integer': + $this->m = $args[0]; + $this->n = $args[1]; + $this->A = array_fill(0, $this->m, array_fill(0, $this->n, $args[2])); + break; + //Rectangular matrix constant-filled - m x n filled with c + case 'integer,integer,double': + $this->m = $args[0]; + $this->n = $args[1]; + $this->A = array_fill(0, $this->m, array_fill(0, $this->n, $args[2])); + break; + //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; + //Rectangular matrix - m x n initialized from 2D array + case 'array,integer,integer': + $this->m = $args[1]; + $this->n = $args[2]; + $this->A = $args[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 Exception(JAMAError(ArrayLengthException)); + } + break; + default: + throw new Exception(JAMAError(PolymorphicArgumentException)); + break; + } + } else { + throw new Exception(JAMAError(PolymorphicArgumentException)); + } + } // function __construct() + + + /** + * getArray + * + * @return array Matrix array + */ + public function getArray() { + return $this->A; + } // function getArray() + + + /** + * getArrayCopy + * + * @return array Matrix array copy + */ + public function getArrayCopy() { + return $this->A; + } // function getArrayCopy() + + + /** + * constructWithCopy + * Construct a matrix from a copy of a 2-D array. + * + * @param double A[][] Two-dimensional array of doubles. + * @exception IllegalArgumentException All rows must have the same length + */ + public function constructWithCopy($A) { + $this->m = count($A); + $this->n = count($A[0]); + $newCopyMatrix = new Matrix($this->m, $this->n); + for ($i = 0; $i < $this->m; ++$i) { + if (count($A[$i]) != $this->n) { + throw new Exception(JAMAError(RowLengthException)); + } + for ($j = 0; $j < $this->n; ++$j) { + $newCopyMatrix->A[$i][$j] = $A[$i][$j]; + } + } + return $newCopyMatrix; + } // function constructWithCopy() + + + /** + * getColumnPackedCopy + * + * Get a column-packed array + * @return array Column-packed matrix array + */ + public function getColumnPackedCopy() { + $P = array(); + for($i = 0; $i < $this->m; ++$i) { + for($j = 0; $j < $this->n; ++$j) { + array_push($P, $this->A[$j][$i]); + } + } + return $P; + } // function getColumnPackedCopy() + + + /** + * getRowPackedCopy + * + * Get a row-packed array + * @return array Row-packed matrix array + */ + public function getRowPackedCopy() { + $P = array(); + for($i = 0; $i < $this->m; ++$i) { + for($j = 0; $j < $this->n; ++$j) { + array_push($P, $this->A[$i][$j]); + } + } + return $P; + } // function getRowPackedCopy() + + + /** + * 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 Exception(JAMAError(ArgumentBoundsException)); } + if ($j0 >= 0) { $n = $this->n - $j0; } else { throw new Exception(JAMAError(ArgumentBoundsException)); } + $R = new 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 Exception(JAMAError(ArgumentBoundsException)); } + if (($jF > $j0) && ($this->n >= $jF) && ($j0 >= 0)) { $n = $jF - $j0; } else { throw new Exception(JAMAError(ArgumentBoundsException)); } + $R = new 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 Exception(JAMAError(ArgumentBoundsException)); } + if (count($CL) > 0) { $n = count($CL); } else { throw new Exception(JAMAError(ArgumentBoundsException)); } + $R = new 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 Exception(JAMAError(ArgumentBoundsException)); } + if (count($CL) > 0) { $n = count($CL); } else { throw new Exception(JAMAError(ArgumentBoundsException)); } + $R = new 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 Exception(JAMAError(ArgumentBoundsException)); } + if (count($CL) > 0) { $n = count($CL); } else { throw new Exception(JAMAError(ArgumentBoundsException)); } + $R = new 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 Exception(JAMAError(ArgumentBoundsException)); } + if (($jF >= $j0) && ($this->n >= $jF) && ($j0 >= 0)) { $n = $jF - $j0; } else { throw new Exception(JAMAError(ArgumentBoundsException)); } + $R = new 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 Exception(JAMAError(PolymorphicArgumentException)); + break; + } + } else { + throw new Exception(JAMAError(PolymorphicArgumentException)); + } + } // function getMatrix() + + + /** + * setMatrix + * + * Set a submatrix + * @param int $i0 Initial row index + * @param int $j0 Initial column index + * @param mixed $S Matrix/Array submatrix + * ($i0, $j0, $S) $S = Matrix + * ($i0, $j0, $S) $S = Array + */ + public function setMatrix() { + if (func_num_args() > 0) { + $args = func_get_args(); + $match = implode(",", array_map('gettype', $args)); + + switch($match) { + case 'integer,integer,object': + if ($args[2] instanceof Matrix) { $M = $args[2]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + if (($args[0] + $M->m) <= $this->m) { $i0 = $args[0]; } else { throw new Exception(JAMAError(ArgumentBoundsException)); } + if (($args[1] + $M->n) <= $this->n) { $j0 = $args[1]; } else { throw new Exception(JAMAError(ArgumentBoundsException)); } + for($i = $i0; $i < $i0 + $M->m; ++$i) { + for($j = $j0; $j < $j0 + $M->n; ++$j) { + $this->A[$i][$j] = $M->get($i - $i0, $j - $j0); + } + } + break; + case 'integer,integer,array': + $M = new Matrix($args[2]); + if (($args[0] + $M->m) <= $this->m) { $i0 = $args[0]; } else { throw new Exception(JAMAError(ArgumentBoundsException)); } + if (($args[1] + $M->n) <= $this->n) { $j0 = $args[1]; } else { throw new Exception(JAMAError(ArgumentBoundsException)); } + for($i = $i0; $i < $i0 + $M->m; ++$i) { + for($j = $j0; $j < $j0 + $M->n; ++$j) { + $this->A[$i][$j] = $M->get($i - $i0, $j - $j0); + } + } + break; + default: + throw new Exception(JAMAError(PolymorphicArgumentException)); + break; + } + } else { + throw new Exception(JAMAError(PolymorphicArgumentException)); + } + } // function setMatrix() + + + /** + * checkMatrixDimensions + * + * Is matrix B the same size? + * @param Matrix $B Matrix B + * @return boolean + */ + public function checkMatrixDimensions($B = null) { + if ($B instanceof Matrix) { + if (($this->m == $B->getRowDimension()) && ($this->n == $B->getColumnDimension())) { + return true; + } else { + throw new Exception(JAMAError(MatrixDimensionException)); + } + } else { + throw new Exception(JAMAError(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; + /* + if (is_int($i) && is_int($j) && is_numeric($c)) { + if (($i < $this->m) && ($j < $this->n)) { + $this->A[$i][$j] = $c; + } else { + echo "A[$i][$j] = $c
"; + throw new Exception(JAMAError(ArgumentBoundsException)); + } + } else { + throw new Exception(JAMAError(ArgumentTypeException)); + } + */ + } // 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 Matrix($m, $n); + for($i = 0; $i < $m; ++$i) { + $R->set($i, $i, $c); + } + return $R; + } // function diagonal() + + + /** + * filled + * + * Generate a filled matrix + * @param int $m Row dimension + * @param int $n Column dimension + * @param int $c Fill constant + * @return Matrix Filled matrix + */ + public function filled($m = null, $n = null, $c = 0) { + if (is_int($m) && is_int($n) && is_numeric($c)) { + $R = new Matrix($m, $n, $c); + return $R; + } else { + throw new Exception(JAMAError(ArgumentTypeException)); + } + } // function filled() + + /** + * random + * + * Generate a random matrix + * @param int $m Row dimension + * @param int $n Column dimension + * @return Matrix Random matrix + */ + public function random($m = null, $n = null, $a = RAND_MIN, $b = RAND_MAX) { + if (is_int($m) && is_int($n) && is_numeric($a) && is_numeric($b)) { + $R = new Matrix($m, $n); + for($i = 0; $i < $m; ++$i) { + for($j = 0; $j < $n; ++$j) { + $R->set($i, $j, mt_rand($a, $b)); + } + } + return $R; + } else { + throw new Exception(JAMAError(ArgumentTypeException)); + } + } // function random() + + + /** + * packed + * + * Alias for getRowPacked + * @return array Packed array + */ + public function packed() { + return $this->getRowPacked(); + } // function packed() + + + /** + * 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 Exception(JAMAError(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 Exception(JAMAError(ArgumentTypeException)); + } + } // function getMatrixByCol() + + + /** + * transpose + * + * Tranpose matrix + * @return Matrix Transposed matrix + */ + public function transpose() { + $R = new 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() + + + /** + * norm1 + * + * One norm + * @return float Maximum column sum + */ + public function norm1() { + $r = 0; + for($j = 0; $j < $this->n; ++$j) { + $s = 0; + for($i = 0; $i < $this->m; ++$i) { + $s += abs($this->A[$i][$j]); + } + $r = ($r > $s) ? $r : $s; + } + return $r; + } // function norm1() + + + /** + * norm2 + * + * Maximum singular value + * @return float Maximum singular value + */ + public function norm2() { + } // function norm2() + + + /** + * normInf + * + * Infinite norm + * @return float Maximum row sum + */ + public function normInf() { + $r = 0; + for($i = 0; $i < $this->m; ++$i) { + $s = 0; + for($j = 0; $j < $this->n; ++$j) { + $s += abs($this->A[$i][$j]); + } + $r = ($r > $s) ? $r : $s; + } + return $r; + } // function normInf() + + + /** + * normF + * + * Frobenius norm + * @return float Square root of the sum of all elements squared + */ + public function normF() { + $f = 0; + for ($i = 0; $i < $this->m; ++$i) { + for ($j = 0; $j < $this->n; ++$j) { + $f = hypo($f,$this->A[$i][$j]); + } + } + return $f; + } // function normF() + + + /** + * Matrix rank + * + * @return effective numerical rank, obtained from SVD. + */ + public function rank () { + $svd = new SingularValueDecomposition($this); + return $svd->rank(); + } // function rank () + + + /** + * Matrix condition (2 norm) + * + * @return ratio of largest to smallest singular value. + */ + public function cond () { + $svd = new SingularValueDecomposition($this); + return $svd->cond(); + } // function cond () + + + /** + * 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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(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 Exception(JAMAError(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(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); + } + } + return $this; + } else { + throw new Exception(JAMAError(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(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 Exception(JAMAError(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(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); + } + } + return $this; + } else { + throw new Exception(JAMAError(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(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 Exception(JAMAError(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(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); + } + } + return $this; + } else { + throw new Exception(JAMAError(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(PolymorphicArgumentException)); + break; + } + $this->checkMatrixDimensions($M); + for($i = 0; $i < $this->m; ++$i) { + for($j = 0; $j < $this->n; ++$j) { + $M->set($i, $j, $this->A[$i][$j] / $M->get($i, $j)); + } + } + return $M; + } else { + throw new Exception(JAMAError(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(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 Exception(JAMAError(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(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 Exception(JAMAError(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(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 Exception(JAMAError(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 Matrix) { $B = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + if ($this->n == $B->m) { + $C = new 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 Exception(JAMAError(MatrixDimensionMismatch)); + } + break; + case 'array': + $B = new Matrix($args[0]); + if ($this->n == $B->m) { + $C = new 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 Exception(JAMAError(MatrixDimensionMismatch)); + } + return $M; + break; + case 'integer': + $C = new 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 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 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 Exception(JAMAError(PolymorphicArgumentException)); + break; + } + } else { + throw new Exception(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + break; + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(PolymorphicArgumentException)); + break; + } + $this->checkMatrixDimensions($M); + for($i = 0; $i < $this->m; ++$i) { + for($j = 0; $j < $this->n; ++$j) { + $this->A[$i][$j] = pow($this->A[$i][$j],$M->get($i, $j)); + } + } + return $this; + } else { + throw new Exception(JAMAError(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 Matrix) { $M = $args[0]; } else { throw new Exception(JAMAError(ArgumentTypeException)); } + case 'array': + $M = new Matrix($args[0]); + break; + default: + throw new Exception(JAMAError(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),'"').'"'; + $this->A[$i][$j] = trim($this->A[$i][$j],'"').trim($M->get($i, $j),'"'); + } + } + return $this; + } else { + throw new Exception(JAMAError(PolymorphicArgumentException)); + } + } // function concat() + + + /** + * chol + * + * Cholesky decomposition + * @return Matrix Cholesky decomposition + */ + public function chol() { + return new CholeskyDecomposition($this); + } // function chol() + + + /** + * lu + * + * LU decomposition + * @return Matrix LU decomposition + */ + public function lu() { + return new LUDecomposition($this); + } // function lu() + + + /** + * qr + * + * QR decomposition + * @return Matrix QR decomposition + */ + public function qr() { + return new QRDecomposition($this); + } // function qr() + + + /** + * eig + * + * Eigenvalue decomposition + * @return Matrix Eigenvalue decomposition + */ + public function eig() { + return new EigenvalueDecomposition($this); + } // function eig() + + + /** + * svd + * + * Singular value decomposition + * @return Singular value decomposition + */ + public function svd() { + return new SingularValueDecomposition($this); + } // function svd() + + + /** + * 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 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 LUDecomposition($this); + return $L->det(); + } // function det() + + + /** + * Older debugging utility for backwards compatability. + * + * @return html version of matrix + */ + public function mprint($A, $format="%01.2f", $width=2) { + $m = count($A); + $n = count($A[0]); + $spacing = str_repeat(' ',$width); + + for ($i = 0; $i < $m; ++$i) { + for ($j = 0; $j < $n; ++$j) { + $formatted = sprintf($format, $A[$i][$j]); + echo $formatted.$spacing; + } + echo "
"; + } + } // function mprint() + + + /** + * Debugging utility. + * + * @return Output HTML representation of matrix + */ + public function toHTML($width=2) { + print(''); + for($i = 0; $i < $this->m; ++$i) { + print(''); + for($j = 0; $j < $this->n; ++$j) { + print(''); + } + print(''); + } + print('
' . $this->A[$i][$j] . '
'); + } // function toHTML() + +} // class Matrix diff --git a/libraries/PHPExcel/PHPExcel/Shared/JAMA/QRDecomposition.php b/libraries/PHPExcel/PHPExcel/Shared/JAMA/QRDecomposition.php index f86b4e0f8..80680594d 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/JAMA/QRDecomposition.php +++ b/libraries/PHPExcel/PHPExcel/Shared/JAMA/QRDecomposition.php @@ -1,195 +1,232 @@ = 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 -*/ +/** + * @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 QRDecomposition { - /** - * Array for internal storage of decomposition. - * @var array - */ - var $QR = array(); - /** - * Row dimension. - * @var integer - */ - var $m; + /** + * Array for internal storage of decomposition. + * @var array + */ + private $QR = array(); - /** - * Column dimension. - * @var integer - */ - var $n; + /** + * Row dimension. + * @var integer + */ + private $m; - /** - * Array for internal storage of diagonal of R. - * @var array - */ - var $Rdiag = array(); + /** + * Column dimension. + * @var integer + */ + private $n; - /** - * QR Decomposition computed by Householder reflections. - * @param matrix $A Rectangular matrix - * @return Structure to access R and the Householder vectors and compute Q. - */ - function QRDecomposition($A) { - if( is_a($A, '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 - trigger_error(ArgumentTypeException, ERROR); - } + /** + * Array for internal storage of diagonal of R. + * @var array + */ + private $Rdiag = array(); - /** - * Is the matrix full rank? - * @return boolean true if R, and hence A, has full rank, else false. - */ - function isFullRank() { - for ($j = 0; $j < $this->n; $j++) - if ($this->Rdiag[$j] == 0) - return false; - return true; - } - /** - * Return the Householder vectors - * @return Matrix Lower trapezoidal matrix whose columns define the reflections - */ - 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 Matrix($H); - } + /** + * 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 Matrix) { + // Initialize. + $this->QR = $A->getArrayCopy(); + $this->m = $A->getRowDimension(); + $this->n = $A->getColumnDimension(); + // Main loop. + for ($k = 0; $k < $this->n; ++$k) { + // Compute 2-norm of k-th column without under/overflow. + $nrm = 0.0; + for ($i = $k; $i < $this->m; ++$i) { + $nrm = hypo($nrm, $this->QR[$i][$k]); + } + if ($nrm != 0.0) { + // Form k-th Householder vector. + if ($this->QR[$k][$k] < 0) { + $nrm = -$nrm; + } + for ($i = $k; $i < $this->m; ++$i) { + $this->QR[$i][$k] /= $nrm; + } + $this->QR[$k][$k] += 1.0; + // Apply transformation to remaining columns. + for ($j = $k+1; $j < $this->n; ++$j) { + $s = 0.0; + for ($i = $k; $i < $this->m; ++$i) { + $s += $this->QR[$i][$k] * $this->QR[$i][$j]; + } + $s = -$s/$this->QR[$k][$k]; + for ($i = $k; $i < $this->m; ++$i) { + $this->QR[$i][$j] += $s * $this->QR[$i][$k]; + } + } + } + $this->Rdiag[$k] = -$nrm; + } + } else { + throw new Exception(JAMAError(ArgumentTypeException)); + } + } // function __construct() - /** - * Return the upper triangular factor - * @return Matrix upper triangular factor - */ - 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]; - else if ($i == $j) - $R[$i][$j] = $this->Rdiag[$i]; - else - $R[$i][$j] = 0.0; - } - } - return new Matrix($R); - } - /** - * Generate and return the (economy-sized) orthogonal factor - * @return Matrix orthogonal factor - */ - 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 Matrix($Q); - } - - /** - * 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. - */ - 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 Matrix($X); - return ($X->getMatrix(0, $this->n-1, 0, $nx)); - } else - trigger_error(MatrixRankException, ERROR); - } else - trigger_error(MatrixDimensionException, ERROR); - } -} + /** + * 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 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 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 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 Matrix($X); + return ($X->getMatrix(0, $this->n-1, 0, $nx)); + } else { + throw new Exception(JAMAError(MatrixRankException)); + } + } else { + throw new Exception(JAMAError(MatrixDimensionException)); + } + } // function solve() + +} // class QRDecomposition diff --git a/libraries/PHPExcel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php b/libraries/PHPExcel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php index 4fec7116f..a4b096c59 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php +++ b/libraries/PHPExcel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php @@ -1,501 +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 -*/ + * @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 - */ - var $U = array(); + /** + * Internal storage of U. + * @var array + */ + private $U = array(); - /** - * Internal storage of V. - * @var array - */ - var $V = array(); + /** + * Internal storage of V. + * @var array + */ + private $V = array(); - /** - * Internal storage of singular values. - * @var array - */ - var $s = array(); + /** + * Internal storage of singular values. + * @var array + */ + private $s = array(); - /** - * Row dimension. - * @var int - */ - var $m; + /** + * Row dimension. + * @var int + */ + private $m; - /** - * Column dimension. - * @var int - */ - var $n; + /** + * 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. - */ - function SingularValueDecomposition ($Arg) { - // Initialize. + /** + * 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) { - $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)); + // 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. + // 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) { - 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]; + } - 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]; + } + } - 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 ($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]; + } + } + } + } - 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; - } - } - } + // 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

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; + } + } - 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++; + // Main iteration loop for the singular values. + $pp = $p - 1; + $iter = 0; + $eps = pow(2.0, -52.0); - // 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 + 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; - /* - echo "

Output A

"; - $A = new Matrix($A); - $A->toHTML(); - - echo "

Matrix U

"; - echo "
";
-    print_r($this->U);        
-    echo "
"; - - echo "

Matrix V

"; - echo "
";
-    print_r($this->V);        
-    echo "
"; - - echo "

Vector S

"; - echo "
";
-    print_r($this->s);        
-    echo "
"; - exit; - */ - - } // end constructor - - /** - * Return the left singular vectors - * @access public - * @return U - */ - function getU() { - return new Matrix($this->U, $this->m, min($this->m + 1, $this->n)); - } + // 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 - /** - * Return the right singular vectors - * @access public - * @return V - */ - function getV() { - return new Matrix($this->V); - } + } // end constructor - /** - * Return the one-dimensional array of singular values - * @access public - * @return diagonal of S. - */ - function getSingularValues() { - return $this->s; - } - /** - * Return the diagonal matrix of singular values - * @access public - * @return S - */ - 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); - } + /** + * 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)); + } - /** - * Two norm - * @access public - * @return max(S) - */ - function norm2() { - return $this->s[0]; - } - /** - * Two norm condition number - * @access public - * @return max(S)/min(S) - */ - function cond() { - return $this->s[0] / $this->s[min($this->m, $this->n) - 1]; - } + /** + * Return the right singular vectors + * + * @access public + * @return V + */ + public function getV() { + return new Matrix($this->V); + } - /** - * Effective numerical matrix rank - * @access public - * @return Number of nonnegligible singular values. - */ - 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; - } -} + + /** + * 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 diff --git a/libraries/PHPExcel/PHPExcel/Shared/JAMA/utils/Error.php b/libraries/PHPExcel/PHPExcel/Shared/JAMA/utils/Error.php index 29628a8e7..e73252b3d 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/JAMA/utils/Error.php +++ b/libraries/PHPExcel/PHPExcel/Shared/JAMA/utils/Error.php @@ -1,120 +1,82 @@ Error: ' . $error[$lang][$num] . '
' . $file . ' @ L' . $line . ''; - die(); - break; - - case WARNING: - echo '
Warning: ' . $error[$lang][$num] . '
' . $file . ' @ L' . $line . '
'; - break; - - case NOTICE: - //echo '
Notice: ' . $error[$lang][$num] . '
' . $file . ' @ L' . $line . '
'; - break; + * Custom error handler + * @param int $num Error number + */ +function JAMAError($errorNumber = null) { + global $error; - case E_NOTICE: - //echo '
Notice: ' . $error[$lang][$num] . '
' . $file . ' @ L' . $line . '
'; - break; - - case E_STRICT: - break; - - case E_WARNING: - break; - - default: - echo "
Unknown Error Type: $type - $file @ L{$line}
"; - die(); - break; - } - } else { - die( "Invalid arguments to JAMAError()" ); - } + if (isset($errorNumber)) { + if (isset($error[JAMALANG][$errorNumber])) { + return $error[JAMALANG][$errorNumber]; + } else { + return $error['EN'][$errorNumber]; + } + } else { + return ("Invalid argument to JAMAError()"); + } } - -// TODO MarkBaker -//set_error_handler('JAMAError'); -//error_reporting(ERROR | WARNING); - diff --git a/libraries/PHPExcel/PHPExcel/Shared/JAMA/utils/Maths.php b/libraries/PHPExcel/PHPExcel/Shared/JAMA/utils/Maths.php index dfe7733cb..f5e2a3721 100644 --- a/libraries/PHPExcel/PHPExcel/Shared/JAMA/utils/Maths.php +++ b/libraries/PHPExcel/PHPExcel/Shared/JAMA/utils/Maths.php @@ -1,40 +1,43 @@ abs($b)) { - $r = $b/$a; - $r = abs($a)* sqrt(1+$r*$r); - } else if ($b != 0) { - $r = $a/$b; - $r = abs($b)*sqrt(1+$r*$r); - } else - $r = 0.0; - return $r; -} + 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. -* + * 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 - trigger_error(ArgumentTypeException, ERROR); - } - return sqrt($s); + $s = 0; + foreach (func_get_args() as $d) { + if (is_numeric($d)) { + $s += pow($d, 2); + } else { + throw new Exception(JAMAError(ArgumentTypeException)); + } + } + return sqrt($s); } */