upgrade to PHPExcel 1.7.0
This commit is contained in:
Dieter Adriaenssens
@@ -22,7 +22,7 @@
|
|||||||
* @package PHPExcel_Shared_Escher
|
* @package PHPExcel_Shared_Escher
|
||||||
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
||||||
* @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL
|
* @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
|
||||||
*/
|
*/
|
||||||
|
|
||||||
/**
|
/**
|
||||||
|
|||||||
@@ -22,7 +22,7 @@
|
|||||||
* @package PHPExcel_Shared_Escher
|
* @package PHPExcel_Shared_Escher
|
||||||
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
||||||
* @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL
|
* @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
|
||||||
*/
|
*/
|
||||||
|
|
||||||
/**
|
/**
|
||||||
|
|||||||
@@ -22,7 +22,7 @@
|
|||||||
* @package PHPExcel_Shared_Escher
|
* @package PHPExcel_Shared_Escher
|
||||||
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
||||||
* @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL
|
* @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
|
||||||
*/
|
*/
|
||||||
|
|
||||||
/**
|
/**
|
||||||
|
|||||||
@@ -22,7 +22,7 @@
|
|||||||
* @package PHPExcel_Shared_Escher
|
* @package PHPExcel_Shared_Escher
|
||||||
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
||||||
* @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL
|
* @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
|
||||||
*/
|
*/
|
||||||
|
|
||||||
/**
|
/**
|
||||||
|
|||||||
@@ -22,7 +22,7 @@
|
|||||||
* @package PHPExcel_Shared_Escher
|
* @package PHPExcel_Shared_Escher
|
||||||
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
||||||
* @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL
|
* @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
|
||||||
*/
|
*/
|
||||||
|
|
||||||
/**
|
/**
|
||||||
|
|||||||
@@ -22,7 +22,7 @@
|
|||||||
* @package PHPExcel_Shared_Escher
|
* @package PHPExcel_Shared_Escher
|
||||||
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
||||||
* @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL
|
* @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
|
||||||
*/
|
*/
|
||||||
|
|
||||||
/**
|
/**
|
||||||
|
|||||||
@@ -22,7 +22,7 @@
|
|||||||
* @package PHPExcel_Shared_Escher
|
* @package PHPExcel_Shared_Escher
|
||||||
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
|
||||||
* @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL
|
* @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
|
||||||
*/
|
*/
|
||||||
|
|
||||||
/**
|
/**
|
||||||
|
|||||||
@@ -1,133 +1,149 @@
|
|||||||
<?php
|
<?php
|
||||||
/**
|
/**
|
||||||
* @package JAMA
|
* @package JAMA
|
||||||
*
|
*
|
||||||
* Cholesky decomposition class
|
* Cholesky decomposition class
|
||||||
*
|
*
|
||||||
* For a symmetric, positive definite matrix A, the Cholesky decomposition
|
* For a symmetric, positive definite matrix A, the Cholesky decomposition
|
||||||
* is an lower triangular matrix L so that A = L*L'.
|
* is an lower triangular matrix L so that A = L*L'.
|
||||||
*
|
*
|
||||||
* If the matrix is not symmetric or positive definite, the constructor
|
* If the matrix is not symmetric or positive definite, the constructor
|
||||||
* returns a partial decomposition and sets an internal flag that may
|
* returns a partial decomposition and sets an internal flag that may
|
||||||
* be queried by the isSPD() method.
|
* be queried by the isSPD() method.
|
||||||
*
|
*
|
||||||
* @author Paul Meagher
|
* @author Paul Meagher
|
||||||
* @author Michael Bommarito
|
* @author Michael Bommarito
|
||||||
* @version 1.2
|
* @version 1.2
|
||||||
*/
|
*/
|
||||||
class CholeskyDecomposition {
|
class CholeskyDecomposition {
|
||||||
/**
|
|
||||||
* Decomposition storage
|
|
||||||
* @var array
|
|
||||||
* @access private
|
|
||||||
*/
|
|
||||||
var $L = array();
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Matrix row and column dimension
|
* Decomposition storage
|
||||||
* @var int
|
* @var array
|
||||||
* @access private
|
* @access private
|
||||||
*/
|
*/
|
||||||
var $m;
|
private $L = array();
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Symmetric positive definite flag
|
* Matrix row and column dimension
|
||||||
* @var boolean
|
* @var int
|
||||||
* @access private
|
* @access private
|
||||||
*/
|
*/
|
||||||
var $isspd = true;
|
private $m;
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* CholeskyDecomposition
|
* Symmetric positive definite flag
|
||||||
* Class constructor - decomposes symmetric positive definite matrix
|
* @var boolean
|
||||||
* @param mixed Matrix square symmetric positive definite matrix
|
* @access private
|
||||||
*/
|
*/
|
||||||
function CholeskyDecomposition( $A = null ) {
|
private $isspd = true;
|
||||||
if( is_a($A, '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 ) {
|
* CholeskyDecomposition
|
||||||
$this->L[$i][$i] = sqrt( $sum );
|
*
|
||||||
} else {
|
* Class constructor - decomposes symmetric positive definite matrix
|
||||||
$this->isspd = false;
|
* @param mixed Matrix square symmetric positive definite matrix
|
||||||
}
|
*/
|
||||||
} else {
|
public function __construct($A = null) {
|
||||||
if( $this->L[$i][$i] != 0 )
|
if ($A instanceof Matrix) {
|
||||||
$this->L[$j][$i] = $sum / $this->L[$i][$i];
|
$this->L = $A->getArray();
|
||||||
}
|
$this->m = $A->getRowDimension();
|
||||||
}
|
|
||||||
|
|
||||||
for ($k = $i+1; $k < $this->m; $k++)
|
for($i = 0; $i < $this->m; ++$i) {
|
||||||
$this->L[$i][$k] = 0.0;
|
for($j = $i; $j < $this->m; ++$j) {
|
||||||
}
|
for($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
|
||||||
} else {
|
$sum -= $this->L[$i][$k] * $this->L[$j][$k];
|
||||||
trigger_error(ArgumentTypeException, ERROR);
|
}
|
||||||
}
|
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) {
|
||||||
* Is the matrix symmetric and positive definite?
|
$this->L[$i][$k] = 0.0;
|
||||||
* @return boolean
|
}
|
||||||
*/
|
}
|
||||||
function isSPD () {
|
} else {
|
||||||
return $this->isspd;
|
throw new Exception(JAMAError(ArgumentTypeException));
|
||||||
}
|
}
|
||||||
|
} // function __construct()
|
||||||
|
|
||||||
/**
|
|
||||||
* getL
|
|
||||||
* Return triangular factor.
|
|
||||||
* @return Matrix Lower triangular matrix
|
|
||||||
*/
|
|
||||||
function getL () {
|
|
||||||
return new Matrix($this->L);
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Solve A*X = B
|
* Is the matrix symmetric and positive definite?
|
||||||
* @param $B Row-equal matrix
|
*
|
||||||
* @return Matrix L * L' * X = B
|
* @return boolean
|
||||||
*/
|
*/
|
||||||
function solve ( $B = null ) {
|
public function isSPD() {
|
||||||
if( is_a($B, 'Matrix') ) {
|
return $this->isspd;
|
||||||
if ($B->getRowDimension() == $this->m) {
|
} // function isSPD()
|
||||||
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];
|
* getL
|
||||||
}
|
*
|
||||||
|
* Return triangular factor.
|
||||||
|
* @return Matrix Lower triangular matrix
|
||||||
|
*/
|
||||||
|
public function getL() {
|
||||||
|
return new Matrix($this->L);
|
||||||
|
} // function getL()
|
||||||
|
|
||||||
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++)
|
* Solve A*X = B
|
||||||
$X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
|
*
|
||||||
}
|
* @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();
|
||||||
|
|
||||||
return new Matrix($X, $this->m, $nx);
|
for ($k = 0; $k < $this->m; ++$k) {
|
||||||
} else {
|
for ($i = $k + 1; $i < $this->m; ++$i) {
|
||||||
trigger_error(MatrixSPDException, ERROR);
|
for ($j = 0; $j < $nx; ++$j) {
|
||||||
}
|
$X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
|
||||||
} else {
|
}
|
||||||
trigger_error(MatrixDimensionException, ERROR);
|
}
|
||||||
}
|
for ($j = 0; $j < $nx; ++$j) {
|
||||||
} else {
|
$X[$k][$j] /= $this->L[$k][$k];
|
||||||
trigger_error(ArgumentTypeException, ERROR);
|
}
|
||||||
}
|
}
|
||||||
}
|
|
||||||
}
|
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
|
||||||
|
|||||||
File diff suppressed because it is too large
Load Diff
@@ -1,222 +1,255 @@
|
|||||||
<?php
|
<?php
|
||||||
/**
|
/**
|
||||||
* @package JAMA
|
* @package JAMA
|
||||||
*
|
*
|
||||||
* For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
|
* 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,
|
* 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.
|
* 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.
|
* 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
|
* The LU decompostion with pivoting always exists, even if the matrix is
|
||||||
* singular, so the constructor will never fail. The primary use of the
|
* singular, so the constructor will never fail. The primary use of the
|
||||||
* LU decomposition is in the solution of square systems of simultaneous
|
* LU decomposition is in the solution of square systems of simultaneous
|
||||||
* linear equations. This will fail if isNonsingular() returns false.
|
* linear equations. This will fail if isNonsingular() returns false.
|
||||||
*
|
*
|
||||||
* @author Paul Meagher
|
* @author Paul Meagher
|
||||||
* @author Bartosz Matosiuk
|
* @author Bartosz Matosiuk
|
||||||
* @author Michael Bommarito
|
* @author Michael Bommarito
|
||||||
* @version 1.1
|
* @version 1.1
|
||||||
* @license PHP v3.0
|
* @license PHP v3.0
|
||||||
*/
|
*/
|
||||||
class LUDecomposition {
|
class LUDecomposition {
|
||||||
/**
|
|
||||||
* Decomposition storage
|
|
||||||
* @var array
|
|
||||||
*/
|
|
||||||
var $LU = array();
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Row dimension.
|
* Decomposition storage
|
||||||
* @var int
|
* @var array
|
||||||
*/
|
*/
|
||||||
var $m;
|
private $LU = array();
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Column dimension.
|
* Row dimension.
|
||||||
* @var int
|
* @var int
|
||||||
*/
|
*/
|
||||||
var $n;
|
private $m;
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Pivot sign.
|
* Column dimension.
|
||||||
* @var int
|
* @var int
|
||||||
*/
|
*/
|
||||||
var $pivsign;
|
private $n;
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Internal storage of pivot vector.
|
* Pivot sign.
|
||||||
* @var array
|
* @var int
|
||||||
*/
|
*/
|
||||||
var $piv = array();
|
private $pivsign;
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* LU Decomposition constructor.
|
* Internal storage of pivot vector.
|
||||||
* @param $A Rectangular matrix
|
* @var array
|
||||||
* @return Structure to access L, U and piv.
|
*/
|
||||||
*/
|
private $piv = array();
|
||||||
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);
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Get upper triangular factor.
|
* LU Decomposition constructor.
|
||||||
* @return array Upper triangular factor
|
*
|
||||||
*/
|
* @param $A Rectangular matrix
|
||||||
function getU () {
|
* @return Structure to access L, U and piv.
|
||||||
for ($i = 0; $i < $this->n; $i++) {
|
*/
|
||||||
for ($j = 0; $j < $this->n; $j++) {
|
public function __construct($A) {
|
||||||
if ($i <= $j)
|
if ($A instanceof Matrix) {
|
||||||
$U[$i][$j] = $this->LU[$i][$j];
|
// Use a "left-looking", dot-product, Crout/Doolittle algorithm.
|
||||||
else
|
$this->LU = $A->getArrayCopy();
|
||||||
$U[$i][$j] = 0.0;
|
$this->m = $A->getRowDimension();
|
||||||
}
|
$this->n = $A->getColumnDimension();
|
||||||
}
|
for ($i = 0; $i < $this->m; ++$i) {
|
||||||
return new Matrix($U);
|
$this->piv[$i] = $i;
|
||||||
}
|
}
|
||||||
|
$this->pivsign = 1;
|
||||||
|
$LUrowi = $LUcolj = array();
|
||||||
|
|
||||||
/**
|
// Outer loop.
|
||||||
* Return pivot permutation vector.
|
for ($j = 0; $j < $this->n; ++$j) {
|
||||||
* @return array Pivot vector
|
// Make a copy of the j-th column to localize references.
|
||||||
*/
|
for ($i = 0; $i < $this->m; ++$i) {
|
||||||
function getPivot () {
|
$LUcolj[$i] = &$this->LU[$i][$j];
|
||||||
return $this->piv;
|
}
|
||||||
}
|
// 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()
|
||||||
|
|
||||||
/**
|
|
||||||
* Alias for getPivot
|
|
||||||
* @see getPivot
|
|
||||||
*/
|
|
||||||
function getDoublePivot () {
|
|
||||||
return $this->getPivot();
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Is the matrix nonsingular?
|
* Get lower triangular factor.
|
||||||
* @return true if U, and hence A, is nonsingular.
|
*
|
||||||
*/
|
* @return array Lower triangular factor
|
||||||
function isNonsingular () {
|
*/
|
||||||
for ($j = 0; $j < $this->n; $j++) {
|
public function getL() {
|
||||||
if ($this->LU[$j][$j] == 0)
|
for ($i = 0; $i < $this->m; ++$i) {
|
||||||
return false;
|
for ($j = 0; $j < $this->n; ++$j) {
|
||||||
}
|
if ($i > $j) {
|
||||||
return true;
|
$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()
|
||||||
|
|
||||||
/**
|
|
||||||
* 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
|
* Get upper triangular factor.
|
||||||
* @param $B A Matrix with as many rows as A and any number of columns.
|
*
|
||||||
* @return X so that L*U*X = B(piv,:)
|
* @return array Upper triangular factor
|
||||||
* @exception IllegalArgumentException Matrix row dimensions must agree.
|
*/
|
||||||
* @exception RuntimeException Matrix is singular.
|
public function getU() {
|
||||||
*/
|
for ($i = 0; $i < $this->n; ++$i) {
|
||||||
function solve($B) {
|
for ($j = 0; $j < $this->n; ++$j) {
|
||||||
if ($B->getRowDimension() == $this->m) {
|
if ($i <= $j) {
|
||||||
if ($this->isNonsingular()) {
|
$U[$i][$j] = $this->LU[$i][$j];
|
||||||
// Copy right hand side with pivoting
|
} else {
|
||||||
$nx = $B->getColumnDimension();
|
$U[$i][$j] = 0.0;
|
||||||
$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++)
|
return new Matrix($U);
|
||||||
for ($j = 0; $j < $nx; $j++)
|
} // function getU()
|
||||||
$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++)
|
* Return pivot permutation vector.
|
||||||
$X->A[$k][$j] /= $this->LU[$k][$k];
|
*
|
||||||
for ($i = 0; $i < $k; $i++)
|
* @return array Pivot vector
|
||||||
for ($j = 0; $j < $nx; $j++)
|
*/
|
||||||
$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
|
public function getPivot() {
|
||||||
}
|
return $this->piv;
|
||||||
return $X;
|
} // function getPivot()
|
||||||
} else {
|
|
||||||
trigger_error(MatrixSingularException, ERROR);
|
|
||||||
}
|
/**
|
||||||
} else {
|
* Alias for getPivot
|
||||||
trigger_error(MatrixSquareException, ERROR);
|
*
|
||||||
}
|
* @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
|
||||||
|
|||||||
File diff suppressed because it is too large
Load Diff
@@ -1,195 +1,232 @@
|
|||||||
<?php
|
<?php
|
||||||
/**
|
/**
|
||||||
* @package JAMA
|
* @package JAMA
|
||||||
*
|
*
|
||||||
* For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
|
* 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
|
* orthogonal matrix Q and an n-by-n upper triangular matrix R so that
|
||||||
* A = Q*R.
|
* A = Q*R.
|
||||||
*
|
*
|
||||||
* The QR decompostion always exists, even if the matrix does not have
|
* 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
|
* full rank, so the constructor will never fail. The primary use of the
|
||||||
* QR decomposition is in the least squares solution of nonsquare systems
|
* QR decomposition is in the least squares solution of nonsquare systems
|
||||||
* of simultaneous linear equations. This will fail if isFullRank()
|
* of simultaneous linear equations. This will fail if isFullRank()
|
||||||
* returns false.
|
* returns false.
|
||||||
*
|
*
|
||||||
* @author Paul Meagher
|
* @author Paul Meagher
|
||||||
* @license PHP v3.0
|
* @license PHP v3.0
|
||||||
* @version 1.1
|
* @version 1.1
|
||||||
*/
|
*/
|
||||||
class QRDecomposition {
|
class QRDecomposition {
|
||||||
/**
|
|
||||||
* Array for internal storage of decomposition.
|
|
||||||
* @var array
|
|
||||||
*/
|
|
||||||
var $QR = array();
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Row dimension.
|
* Array for internal storage of decomposition.
|
||||||
* @var integer
|
* @var array
|
||||||
*/
|
*/
|
||||||
var $m;
|
private $QR = array();
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Column dimension.
|
* Row dimension.
|
||||||
* @var integer
|
* @var integer
|
||||||
*/
|
*/
|
||||||
var $n;
|
private $m;
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Array for internal storage of diagonal of R.
|
* Column dimension.
|
||||||
* @var array
|
* @var integer
|
||||||
*/
|
*/
|
||||||
var $Rdiag = array();
|
private $n;
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* QR Decomposition computed by Householder reflections.
|
* Array for internal storage of diagonal of R.
|
||||||
* @param matrix $A Rectangular matrix
|
* @var array
|
||||||
* @return Structure to access R and the Householder vectors and compute Q.
|
*/
|
||||||
*/
|
private $Rdiag = array();
|
||||||
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);
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* 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
|
* QR Decomposition computed by Householder reflections.
|
||||||
* @return Matrix Lower trapezoidal matrix whose columns define the reflections
|
*
|
||||||
*/
|
* @param matrix $A Rectangular matrix
|
||||||
function getH() {
|
* @return Structure to access R and the Householder vectors and compute Q.
|
||||||
for ($i = 0; $i < $this->m; $i++) {
|
*/
|
||||||
for ($j = 0; $j < $this->n; $j++) {
|
public function __construct($A) {
|
||||||
if ($i >= $j)
|
if($A instanceof Matrix) {
|
||||||
$H[$i][$j] = $this->QR[$i][$j];
|
// Initialize.
|
||||||
else
|
$this->QR = $A->getArrayCopy();
|
||||||
$H[$i][$j] = 0.0;
|
$this->m = $A->getRowDimension();
|
||||||
}
|
$this->n = $A->getColumnDimension();
|
||||||
}
|
// Main loop.
|
||||||
return new Matrix($H);
|
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
|
* Is the matrix full rank?
|
||||||
* @return Matrix orthogonal factor
|
*
|
||||||
*/
|
* @return boolean true if R, and hence A, has full rank, else false.
|
||||||
function getQ() {
|
*/
|
||||||
for ($k = $this->n-1; $k >= 0; $k--) {
|
public function isFullRank() {
|
||||||
for ($i = 0; $i < $this->m; $i++)
|
for ($j = 0; $j < $this->n; ++$j) {
|
||||||
$Q[$i][$k] = 0.0;
|
if ($this->Rdiag[$j] == 0) {
|
||||||
$Q[$k][$k] = 1.0;
|
return false;
|
||||||
for ($j = $k; $j < $this->n; $j++) {
|
}
|
||||||
if ($this->QR[$k][$k] != 0) {
|
}
|
||||||
$s = 0.0;
|
return true;
|
||||||
for ($i = $k; $i < $this->m; $i++)
|
} // function isFullRank()
|
||||||
$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 the Householder vectors
|
||||||
* @return Matrix Matrix that minimizes the two norm of Q*R*X-B.
|
*
|
||||||
*/
|
* @return Matrix Lower trapezoidal matrix whose columns define the reflections
|
||||||
function solve($B) {
|
*/
|
||||||
if ($B->getRowDimension() == $this->m) {
|
public function getH() {
|
||||||
if ($this->isFullRank()) {
|
for ($i = 0; $i < $this->m; ++$i) {
|
||||||
// Copy right hand side
|
for ($j = 0; $j < $this->n; ++$j) {
|
||||||
$nx = $B->getColumnDimension();
|
if ($i >= $j) {
|
||||||
$X = $B->getArrayCopy();
|
$H[$i][$j] = $this->QR[$i][$j];
|
||||||
// Compute Y = transpose(Q)*B
|
} else {
|
||||||
for ($k = 0; $k < $this->n; $k++) {
|
$H[$i][$j] = 0.0;
|
||||||
for ($j = 0; $j < $nx; $j++) {
|
}
|
||||||
$s = 0.0;
|
}
|
||||||
for ($i = $k; $i < $this->m; $i++)
|
}
|
||||||
$s += $this->QR[$i][$k] * $X[$i][$j];
|
return new Matrix($H);
|
||||||
$s = -$s/$this->QR[$k][$k];
|
} // function getH()
|
||||||
for ($i = $k; $i < $this->m; $i++)
|
|
||||||
$X[$i][$j] += $s * $this->QR[$i][$k];
|
|
||||||
}
|
/**
|
||||||
}
|
* Return the upper triangular factor
|
||||||
// Solve R*X = Y;
|
*
|
||||||
for ($k = $this->n-1; $k >= 0; $k--) {
|
* @return Matrix upper triangular factor
|
||||||
for ($j = 0; $j < $nx; $j++)
|
*/
|
||||||
$X[$k][$j] /= $this->Rdiag[$k];
|
public function getR() {
|
||||||
for ($i = 0; $i < $k; $i++)
|
for ($i = 0; $i < $this->n; ++$i) {
|
||||||
for ($j = 0; $j < $nx; $j++)
|
for ($j = 0; $j < $this->n; ++$j) {
|
||||||
$X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];
|
if ($i < $j) {
|
||||||
}
|
$R[$i][$j] = $this->QR[$i][$j];
|
||||||
$X = new Matrix($X);
|
} elseif ($i == $j) {
|
||||||
return ($X->getMatrix(0, $this->n-1, 0, $nx));
|
$R[$i][$j] = $this->Rdiag[$i];
|
||||||
} else
|
} else {
|
||||||
trigger_error(MatrixRankException, ERROR);
|
$R[$i][$j] = 0.0;
|
||||||
} else
|
}
|
||||||
trigger_error(MatrixDimensionException, ERROR);
|
}
|
||||||
}
|
}
|
||||||
}
|
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
|
||||||
|
|||||||
@@ -1,501 +1,526 @@
|
|||||||
<?php
|
<?php
|
||||||
/**
|
/**
|
||||||
* @package JAMA
|
* @package JAMA
|
||||||
*
|
*
|
||||||
* For an m-by-n matrix A with m >= n, the singular value decomposition is
|
* 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 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'.
|
* 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
|
* The singular values, sigma[$k] = S[$k][$k], are ordered so that
|
||||||
* sigma[0] >= sigma[1] >= ... >= sigma[n-1].
|
* sigma[0] >= sigma[1] >= ... >= sigma[n-1].
|
||||||
*
|
*
|
||||||
* The singular value decompostion always exists, so the constructor will
|
* The singular value decompostion always exists, so the constructor will
|
||||||
* never fail. The matrix condition number and the effective numerical
|
* never fail. The matrix condition number and the effective numerical
|
||||||
* rank can be computed from this decomposition.
|
* rank can be computed from this decomposition.
|
||||||
*
|
*
|
||||||
* @author Paul Meagher
|
* @author Paul Meagher
|
||||||
* @license PHP v3.0
|
* @license PHP v3.0
|
||||||
* @version 1.1
|
* @version 1.1
|
||||||
*/
|
*/
|
||||||
class SingularValueDecomposition {
|
class SingularValueDecomposition {
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Internal storage of U.
|
* Internal storage of U.
|
||||||
* @var array
|
* @var array
|
||||||
*/
|
*/
|
||||||
var $U = array();
|
private $U = array();
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Internal storage of V.
|
* Internal storage of V.
|
||||||
* @var array
|
* @var array
|
||||||
*/
|
*/
|
||||||
var $V = array();
|
private $V = array();
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Internal storage of singular values.
|
* Internal storage of singular values.
|
||||||
* @var array
|
* @var array
|
||||||
*/
|
*/
|
||||||
var $s = array();
|
private $s = array();
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Row dimension.
|
* Row dimension.
|
||||||
* @var int
|
* @var int
|
||||||
*/
|
*/
|
||||||
var $m;
|
private $m;
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Column dimension.
|
* Column dimension.
|
||||||
* @var int
|
* @var int
|
||||||
*/
|
*/
|
||||||
var $n;
|
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();
|
// Initialize.
|
||||||
$this->m = $Arg->getRowDimension();
|
$A = $Arg->getArrayCopy();
|
||||||
$this->n = $Arg->getColumnDimension();
|
$this->m = $Arg->getRowDimension();
|
||||||
$nu = min($this->m, $this->n);
|
$this->n = $Arg->getColumnDimension();
|
||||||
$e = array();
|
$nu = min($this->m, $this->n);
|
||||||
$work = array();
|
$e = array();
|
||||||
$wantu = true;
|
$work = array();
|
||||||
$wantv = true;
|
$wantu = true;
|
||||||
$nct = min($this->m - 1, $this->n);
|
$wantv = true;
|
||||||
$nrt = max(0, min($this->n - 2, $this->m));
|
$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
|
// Reduce A to bidiagonal form, storing the diagonal elements
|
||||||
// in s and the super-diagonal elements in e.
|
// 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) {
|
for ($j = $k + 1; $j < $this->n; ++$j) {
|
||||||
// Compute the transformation for the k-th column and
|
if (($k < $nct) & ($this->s[$k] != 0.0)) {
|
||||||
// place the k-th diagonal in s[$k].
|
// Apply the transformation.
|
||||||
// Compute 2-norm of k-th column without under/overflow.
|
$t = 0;
|
||||||
$this->s[$k] = 0;
|
for ($i = $k; $i < $this->m; ++$i) {
|
||||||
for ($i = $k; $i < $this->m; $i++)
|
$t += $A[$i][$k] * $A[$i][$j];
|
||||||
$this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
|
}
|
||||||
if ($this->s[$k] != 0.0) {
|
$t = -$t / $A[$k][$k];
|
||||||
if ($A[$k][$k] < 0.0)
|
for ($i = $k; $i < $this->m; ++$i) {
|
||||||
$this->s[$k] = -$this->s[$k];
|
$A[$i][$j] += $t * $A[$i][$k];
|
||||||
for ($i = $k; $i < $this->m; $i++)
|
}
|
||||||
$A[$i][$k] /= $this->s[$k];
|
// Place the k-th row of A into e for the
|
||||||
$A[$k][$k] += 1.0;
|
// subsequent calculation of the row transformation.
|
||||||
}
|
$e[$j] = $A[$k][$j];
|
||||||
$this->s[$k] = -$this->s[$k];
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
for ($j = $k + 1; $j < $this->n; $j++) {
|
if ($wantu AND ($k < $nct)) {
|
||||||
if (($k < $nct) & ($this->s[$k] != 0.0)) {
|
// Place the transformation in U for subsequent back
|
||||||
// Apply the transformation.
|
// multiplication.
|
||||||
$t = 0;
|
for ($i = $k; $i < $this->m; ++$i) {
|
||||||
for ($i = $k; $i < $this->m; $i++)
|
$this->U[$i][$k] = $A[$i][$k];
|
||||||
$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)) {
|
if ($k < $nrt) {
|
||||||
// Place the transformation in U for subsequent back
|
// Compute the k-th row transformation and place the
|
||||||
// multiplication.
|
// k-th super-diagonal in e[$k].
|
||||||
for ($i = $k; $i < $this->m; $i++)
|
// Compute 2-norm without under/overflow.
|
||||||
$this->U[$i][$k] = $A[$i][$k];
|
$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) {
|
// Set up the final bidiagonal matrix or order p.
|
||||||
// Compute the k-th row transformation and place the
|
$p = min($this->n, $this->m + 1);
|
||||||
// k-th super-diagonal in e[$k].
|
if ($nct < $this->n) {
|
||||||
// Compute 2-norm without under/overflow.
|
$this->s[$nct] = $A[$nct][$nct];
|
||||||
$e[$k] = 0;
|
}
|
||||||
for ($i = $k + 1; $i < $this->n; $i++)
|
if ($this->m < $p) {
|
||||||
$e[$k] = hypo($e[$k], $e[$i]);
|
$this->s[$p-1] = 0.0;
|
||||||
if ($e[$k] != 0.0) {
|
}
|
||||||
if ($e[$k+1] < 0.0)
|
if ($nrt + 1 < $p) {
|
||||||
$e[$k] = -$e[$k];
|
$e[$nrt] = $A[$nrt][$p-1];
|
||||||
for ($i = $k + 1; $i < $this->n; $i++)
|
}
|
||||||
$e[$i] /= $e[$k];
|
$e[$p-1] = 0.0;
|
||||||
$e[$k+1] += 1.0;
|
// If required, generate U.
|
||||||
}
|
if ($wantu) {
|
||||||
$e[$k] = -$e[$k];
|
for ($j = $nct; $j < $nu; ++$j) {
|
||||||
if (($k+1 < $this->m) AND ($e[$k] != 0.0)) {
|
for ($i = 0; $i < $this->m; ++$i) {
|
||||||
// Apply the transformation.
|
$this->U[$i][$j] = 0.0;
|
||||||
for ($i = $k+1; $i < $this->m; $i++)
|
}
|
||||||
$work[$i] = 0.0;
|
$this->U[$j][$j] = 1.0;
|
||||||
for ($j = $k+1; $j < $this->n; $j++)
|
}
|
||||||
for ($i = $k+1; $i < $this->m; $i++)
|
for ($k = $nct - 1; $k >= 0; --$k) {
|
||||||
$work[$i] += $e[$j] * $A[$i][$j];
|
if ($this->s[$k] != 0.0) {
|
||||||
for ($j = $k + 1; $j < $this->n; $j++) {
|
for ($j = $k + 1; $j < $nu; ++$j) {
|
||||||
$t = -$e[$j] / $e[$k+1];
|
$t = 0;
|
||||||
for ($i = $k + 1; $i < $this->m; $i++)
|
for ($i = $k; $i < $this->m; ++$i) {
|
||||||
$A[$i][$j] += $t * $work[$i];
|
$t += $this->U[$i][$k] * $this->U[$i][$j];
|
||||||
}
|
}
|
||||||
}
|
$t = -$t / $this->U[$k][$k];
|
||||||
if ($wantv) {
|
for ($i = $k; $i < $this->m; ++$i) {
|
||||||
// Place the transformation in V for subsequent
|
$this->U[$i][$j] += $t * $this->U[$i][$k];
|
||||||
// back multiplication.
|
}
|
||||||
for ($i = $k + 1; $i < $this->n; $i++)
|
}
|
||||||
$this->V[$i][$k] = $e[$i];
|
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.
|
// If required, generate V.
|
||||||
$p = min($this->n, $this->m + 1);
|
if ($wantv) {
|
||||||
if ($nct < $this->n)
|
for ($k = $this->n - 1; $k >= 0; --$k) {
|
||||||
$this->s[$nct] = $A[$nct][$nct];
|
if (($k < $nrt) AND ($e[$k] != 0.0)) {
|
||||||
if ($this->m < $p)
|
for ($j = $k + 1; $j < $nu; ++$j) {
|
||||||
$this->s[$p-1] = 0.0;
|
$t = 0;
|
||||||
if ($nrt + 1 < $p)
|
for ($i = $k + 1; $i < $this->n; ++$i) {
|
||||||
$e[$nrt] = $A[$nrt][$p-1];
|
$t += $this->V[$i][$k]* $this->V[$i][$j];
|
||||||
$e[$p-1] = 0.0;
|
}
|
||||||
// If required, generate U.
|
$t = -$t / $this->V[$k+1][$k];
|
||||||
if ($wantu) {
|
for ($i = $k + 1; $i < $this->n; ++$i) {
|
||||||
for ($j = $nct; $j < $nu; $j++) {
|
$this->V[$i][$j] += $t * $this->V[$i][$k];
|
||||||
for ($i = 0; $i < $this->m; $i++)
|
}
|
||||||
$this->U[$i][$j] = 0.0;
|
}
|
||||||
$this->U[$j][$j] = 1.0;
|
}
|
||||||
}
|
for ($i = 0; $i < $this->n; ++$i) {
|
||||||
for ($k = $nct - 1; $k >= 0; $k--) {
|
$this->V[$i][$k] = 0.0;
|
||||||
if ($this->s[$k] != 0.0) {
|
}
|
||||||
for ($j = $k + 1; $j < $nu; $j++) {
|
$this->V[$k][$k] = 1.0;
|
||||||
$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.
|
// Main iteration loop for the singular values.
|
||||||
if ($wantv) {
|
$pp = $p - 1;
|
||||||
for ($k = $this->n - 1; $k >= 0; $k--) {
|
$iter = 0;
|
||||||
if (($k < $nrt) AND ($e[$k] != 0.0)) {
|
$eps = pow(2.0, -52.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.
|
while ($p > 0) {
|
||||||
$pp = $p - 1;
|
// Here is where a test for too many iterations would go.
|
||||||
$iter = 0;
|
// This section of the program inspects for negligible
|
||||||
$eps = pow(2.0, -52.0);
|
// elements in the s and e arrays. On completion the
|
||||||
while ($p > 0) {
|
// variables kase and k are set as follows:
|
||||||
|
// kase = 1 if s(p) and e[k-1] are negligible and k<p
|
||||||
|
// kase = 2 if s(k) is negligible and k<p
|
||||||
|
// kase = 3 if e[k-1] is negligible, k<p, and
|
||||||
|
// s(k), ..., s(p) are not negligible (qr step).
|
||||||
|
// kase = 4 if e(p-1) is negligible (convergence).
|
||||||
|
for ($k = $p - 2; $k >= -1; --$k) {
|
||||||
|
if ($k == -1) {
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) {
|
||||||
|
$e[$k] = 0.0;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if ($k == $p - 2) {
|
||||||
|
$kase = 4;
|
||||||
|
} else {
|
||||||
|
for ($ks = $p - 1; $ks >= $k; --$ks) {
|
||||||
|
if ($ks == $k) {
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
$t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.);
|
||||||
|
if (abs($this->s[$ks]) <= $eps * $t) {
|
||||||
|
$this->s[$ks] = 0.0;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if ($ks == $k) {
|
||||||
|
$kase = 3;
|
||||||
|
} else if ($ks == $p-1) {
|
||||||
|
$kase = 1;
|
||||||
|
} else {
|
||||||
|
$kase = 2;
|
||||||
|
$k = $ks;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
++$k;
|
||||||
|
|
||||||
// Here is where a test for too many iterations would go.
|
// Perform the task indicated by kase.
|
||||||
// This section of the program inspects for negligible
|
switch ($kase) {
|
||||||
// elements in the s and e arrays. On completion the
|
// Deflate negligible s(p).
|
||||||
// variables kase and k are set as follows:
|
case 1:
|
||||||
// kase = 1 if s(p) and e[k-1] are negligible and k<p
|
$f = $e[$p-2];
|
||||||
// kase = 2 if s(k) is negligible and k<p
|
$e[$p-2] = 0.0;
|
||||||
// kase = 3 if e[k-1] is negligible, k<p, and
|
for ($j = $p - 2; $j >= $k; --$j) {
|
||||||
// s(k), ..., s(p) are not negligible (qr step).
|
$t = hypo($this->s[$j],$f);
|
||||||
// kase = 4 if e(p-1) is negligible (convergence).
|
$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
|
||||||
|
|
||||||
for ($k = $p - 2; $k >= -1; $k--) {
|
} // end constructor
|
||||||
if ($k == -1)
|
|
||||||
break;
|
|
||||||
if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) {
|
|
||||||
$e[$k] = 0.0;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if ($k == $p - 2)
|
|
||||||
$kase = 4;
|
|
||||||
else {
|
|
||||||
for ($ks = $p - 1; $ks >= $k; $ks--) {
|
|
||||||
if ($ks == $k)
|
|
||||||
break;
|
|
||||||
$t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.);
|
|
||||||
if (abs($this->s[$ks]) <= $eps * $t) {
|
|
||||||
$this->s[$ks] = 0.0;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if ($ks == $k)
|
|
||||||
$kase = 3;
|
|
||||||
else if ($ks == $p-1)
|
|
||||||
$kase = 1;
|
|
||||||
else {
|
|
||||||
$kase = 2;
|
|
||||||
$k = $ks;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
$k++;
|
|
||||||
|
|
||||||
// Perform the task indicated by kase.
|
|
||||||
switch ($kase) {
|
|
||||||
// Deflate negligible s(p).
|
|
||||||
case 1:
|
|
||||||
$f = $e[$p-2];
|
|
||||||
$e[$p-2] = 0.0;
|
|
||||||
for ($j = $p - 2; $j >= $k; $j--) {
|
|
||||||
$t = hypo($this->s[$j],$f);
|
|
||||||
$cs = $this->s[$j] / $t;
|
|
||||||
$sn = $f / $t;
|
|
||||||
$this->s[$j] = $t;
|
|
||||||
if ($j != $k) {
|
|
||||||
$f = -$sn * $e[$j-1];
|
|
||||||
$e[$j-1] = $cs * $e[$j-1];
|
|
||||||
}
|
|
||||||
if ($wantv) {
|
|
||||||
for ($i = 0; $i < $this->n; $i++) {
|
|
||||||
$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1];
|
|
||||||
$this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1];
|
|
||||||
$this->V[$i][$j] = $t;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
break;
|
|
||||||
// Split at negligible s(k).
|
|
||||||
case 2:
|
|
||||||
$f = $e[$k-1];
|
|
||||||
$e[$k-1] = 0.0;
|
|
||||||
for ($j = $k; $j < $p; $j++) {
|
|
||||||
$t = hypo($this->s[$j], $f);
|
|
||||||
$cs = $this->s[$j] / $t;
|
|
||||||
$sn = $f / $t;
|
|
||||||
$this->s[$j] = $t;
|
|
||||||
$f = -$sn * $e[$j];
|
|
||||||
$e[$j] = $cs * $e[$j];
|
|
||||||
if ($wantu) {
|
|
||||||
for ($i = 0; $i < $this->m; $i++) {
|
|
||||||
$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1];
|
|
||||||
$this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1];
|
|
||||||
$this->U[$i][$j] = $t;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
break;
|
|
||||||
// Perform one qr step.
|
|
||||||
case 3:
|
|
||||||
// Calculate the shift.
|
|
||||||
$scale = max(max(max(max(
|
|
||||||
abs($this->s[$p-1]),abs($this->s[$p-2])),abs($e[$p-2])),
|
|
||||||
abs($this->s[$k])), abs($e[$k]));
|
|
||||||
$sp = $this->s[$p-1] / $scale;
|
|
||||||
$spm1 = $this->s[$p-2] / $scale;
|
|
||||||
$epm1 = $e[$p-2] / $scale;
|
|
||||||
$sk = $this->s[$k] / $scale;
|
|
||||||
$ek = $e[$k] / $scale;
|
|
||||||
$b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
|
|
||||||
$c = ($sp * $epm1) * ($sp * $epm1);
|
|
||||||
$shift = 0.0;
|
|
||||||
if (($b != 0.0) || ($c != 0.0)) {
|
|
||||||
$shift = sqrt($b * $b + $c);
|
|
||||||
if ($b < 0.0)
|
|
||||||
$shift = -$shift;
|
|
||||||
$shift = $c / ($b + $shift);
|
|
||||||
}
|
|
||||||
$f = ($sk + $sp) * ($sk - $sp) + $shift;
|
|
||||||
$g = $sk * $ek;
|
|
||||||
// Chase zeros.
|
|
||||||
for ($j = $k; $j < $p-1; $j++) {
|
|
||||||
$t = hypo($f,$g);
|
|
||||||
$cs = $f/$t;
|
|
||||||
$sn = $g/$t;
|
|
||||||
if ($j != $k)
|
|
||||||
$e[$j-1] = $t;
|
|
||||||
$f = $cs * $this->s[$j] + $sn * $e[$j];
|
|
||||||
$e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
|
|
||||||
$g = $sn * $this->s[$j+1];
|
|
||||||
$this->s[$j+1] = $cs * $this->s[$j+1];
|
|
||||||
if ($wantv) {
|
|
||||||
for ($i = 0; $i < $this->n; $i++) {
|
|
||||||
$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1];
|
|
||||||
$this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1];
|
|
||||||
$this->V[$i][$j] = $t;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
$t = hypo($f,$g);
|
|
||||||
$cs = $f/$t;
|
|
||||||
$sn = $g/$t;
|
|
||||||
$this->s[$j] = $t;
|
|
||||||
$f = $cs * $e[$j] + $sn * $this->s[$j+1];
|
|
||||||
$this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1];
|
|
||||||
$g = $sn * $e[$j+1];
|
|
||||||
$e[$j+1] = $cs * $e[$j+1];
|
|
||||||
if ($wantu && ($j < $this->m - 1)) {
|
|
||||||
for ($i = 0; $i < $this->m; $i++) {
|
|
||||||
$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1];
|
|
||||||
$this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1];
|
|
||||||
$this->U[$i][$j] = $t;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
$e[$p-2] = $f;
|
|
||||||
$iter = $iter + 1;
|
|
||||||
break;
|
|
||||||
// Convergence.
|
|
||||||
case 4:
|
|
||||||
// Make the singular values positive.
|
|
||||||
if ($this->s[$k] <= 0.0) {
|
|
||||||
$this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);
|
|
||||||
if ($wantv) {
|
|
||||||
for ($i = 0; $i <= $pp; $i++)
|
|
||||||
$this->V[$i][$k] = -$this->V[$i][$k];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// Order the singular values.
|
|
||||||
while ($k < $pp) {
|
|
||||||
if ($this->s[$k] >= $this->s[$k+1])
|
|
||||||
break;
|
|
||||||
$t = $this->s[$k];
|
|
||||||
$this->s[$k] = $this->s[$k+1];
|
|
||||||
$this->s[$k+1] = $t;
|
|
||||||
if ($wantv AND ($k < $this->n - 1)) {
|
|
||||||
for ($i = 0; $i < $this->n; $i++) {
|
|
||||||
$t = $this->V[$i][$k+1];
|
|
||||||
$this->V[$i][$k+1] = $this->V[$i][$k];
|
|
||||||
$this->V[$i][$k] = $t;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if ($wantu AND ($k < $this->m-1)) {
|
|
||||||
for ($i = 0; $i < $this->m; $i++) {
|
|
||||||
$t = $this->U[$i][$k+1];
|
|
||||||
$this->U[$i][$k+1] = $this->U[$i][$k];
|
|
||||||
$this->U[$i][$k] = $t;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
$k++;
|
|
||||||
}
|
|
||||||
$iter = 0;
|
|
||||||
$p--;
|
|
||||||
break;
|
|
||||||
} // end switch
|
|
||||||
} // end while
|
|
||||||
|
|
||||||
/*
|
/**
|
||||||
echo "<p>Output A</p>";
|
* Return the left singular vectors
|
||||||
$A = new Matrix($A);
|
*
|
||||||
$A->toHTML();
|
* @access public
|
||||||
|
* @return U
|
||||||
|
*/
|
||||||
|
public function getU() {
|
||||||
|
return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
|
||||||
|
}
|
||||||
|
|
||||||
echo "<p>Matrix U</p>";
|
|
||||||
echo "<pre>";
|
|
||||||
print_r($this->U);
|
|
||||||
echo "</pre>";
|
|
||||||
|
|
||||||
echo "<p>Matrix V</p>";
|
/**
|
||||||
echo "<pre>";
|
* Return the right singular vectors
|
||||||
print_r($this->V);
|
*
|
||||||
echo "</pre>";
|
* @access public
|
||||||
|
* @return V
|
||||||
|
*/
|
||||||
|
public function getV() {
|
||||||
|
return new Matrix($this->V);
|
||||||
|
}
|
||||||
|
|
||||||
echo "<p>Vector S</p>";
|
|
||||||
echo "<pre>";
|
|
||||||
print_r($this->s);
|
|
||||||
echo "</pre>";
|
|
||||||
exit;
|
|
||||||
*/
|
|
||||||
|
|
||||||
} // end constructor
|
/**
|
||||||
|
* Return the one-dimensional array of singular values
|
||||||
|
*
|
||||||
|
* @access public
|
||||||
|
* @return diagonal of S.
|
||||||
|
*/
|
||||||
|
public function getSingularValues() {
|
||||||
|
return $this->s;
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
|
||||||
* Return the left singular vectors
|
|
||||||
* @access public
|
|
||||||
* @return U
|
|
||||||
*/
|
|
||||||
function getU() {
|
|
||||||
return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Return the right singular vectors
|
* Return the diagonal matrix of singular values
|
||||||
* @access public
|
*
|
||||||
* @return V
|
* @access public
|
||||||
*/
|
* @return S
|
||||||
function getV() {
|
*/
|
||||||
return new Matrix($this->V);
|
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);
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
|
||||||
* 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
|
* Two norm
|
||||||
* @access public
|
*
|
||||||
* @return S
|
* @access public
|
||||||
*/
|
* @return max(S)
|
||||||
function getS() {
|
*/
|
||||||
for ($i = 0; $i < $this->n; $i++) {
|
public function norm2() {
|
||||||
for ($j = 0; $j < $this->n; $j++)
|
return $this->s[0];
|
||||||
$S[$i][$j] = 0.0;
|
}
|
||||||
$S[$i][$i] = $this->s[$i];
|
|
||||||
}
|
|
||||||
return new Matrix($S);
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Two norm
|
|
||||||
* @access public
|
|
||||||
* @return max(S)
|
|
||||||
*/
|
|
||||||
function norm2() {
|
|
||||||
return $this->s[0];
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Two norm condition number
|
* Two norm condition number
|
||||||
* @access public
|
*
|
||||||
* @return max(S)/min(S)
|
* @access public
|
||||||
*/
|
* @return max(S)/min(S)
|
||||||
function cond() {
|
*/
|
||||||
return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
|
public function cond() {
|
||||||
}
|
return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
|
||||||
* Effective numerical matrix rank
|
/**
|
||||||
* @access public
|
* Effective numerical matrix rank
|
||||||
* @return Number of nonnegligible singular values.
|
*
|
||||||
*/
|
* @access public
|
||||||
function rank() {
|
* @return Number of nonnegligible singular values.
|
||||||
$eps = pow(2.0, -52.0);
|
*/
|
||||||
$tol = max($this->m, $this->n) * $this->s[0] * $eps;
|
public function rank() {
|
||||||
$r = 0;
|
$eps = pow(2.0, -52.0);
|
||||||
for ($i = 0; $i < count($this->s); $i++) {
|
$tol = max($this->m, $this->n) * $this->s[0] * $eps;
|
||||||
if ($this->s[$i] > $tol)
|
$r = 0;
|
||||||
$r++;
|
for ($i = 0; $i < count($this->s); ++$i) {
|
||||||
}
|
if ($this->s[$i] > $tol) {
|
||||||
return $r;
|
++$r;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
return $r;
|
||||||
|
}
|
||||||
|
|
||||||
|
} // class SingularValueDecomposition
|
||||||
|
|||||||
@@ -1,120 +1,82 @@
|
|||||||
<?php
|
<?php
|
||||||
/**
|
/**
|
||||||
* @package JAMA
|
* @package JAMA
|
||||||
*
|
*
|
||||||
* Error handling
|
* Error handling
|
||||||
* @author Michael Bommarito
|
* @author Michael Bommarito
|
||||||
* @version 01292005
|
* @version 01292005
|
||||||
*/
|
*/
|
||||||
|
|
||||||
//Language constant
|
//Language constant
|
||||||
define('LANG', 'EN');
|
define('JAMALANG', 'EN');
|
||||||
|
|
||||||
|
|
||||||
//Error type constants
|
|
||||||
define('ERROR', E_USER_ERROR);
|
|
||||||
define('WARNING', E_USER_WARNING);
|
|
||||||
define('NOTICE', E_USER_NOTICE);
|
|
||||||
|
|
||||||
//All errors may be defined by the following format:
|
//All errors may be defined by the following format:
|
||||||
//define('ExceptionName', N);
|
//define('ExceptionName', N);
|
||||||
//$error['lang'][N] = 'Error message';
|
//$error['lang'][ExceptionName] = 'Error message';
|
||||||
$error = array();
|
$error = array();
|
||||||
|
|
||||||
/*
|
/*
|
||||||
I've used Babelfish and a little poor knowledge of Romance/Germanic languages for the translations
|
I've used Babelfish and a little poor knowledge of Romance/Germanic languages for the translations here.
|
||||||
here. Feel free to correct anything that looks amiss to you.
|
Feel free to correct anything that looks amiss to you.
|
||||||
*/
|
*/
|
||||||
|
|
||||||
define('PolymorphicArgumentException', -1);
|
define('PolymorphicArgumentException', -1);
|
||||||
$error['EN'][-1] = "Invalid argument pattern for polymorphic function.";
|
$error['EN'][PolymorphicArgumentException] = "Invalid argument pattern for polymorphic function.";
|
||||||
$error['FR'][-1] = "Modèle inadmissible d'argument pour la fonction polymorphe.".
|
$error['FR'][PolymorphicArgumentException] = "Modèle inadmissible d'argument pour la fonction polymorphe.".
|
||||||
$error['DE'][-1] = "Unzulässiges Argumentmuster für polymorphe Funktion.";
|
$error['DE'][PolymorphicArgumentException] = "Unzulässiges Argumentmuster für polymorphe Funktion.";
|
||||||
|
|
||||||
define('ArgumentTypeException', -2);
|
define('ArgumentTypeException', -2);
|
||||||
$error['EN'][-2] = "Invalid argument type.";
|
$error['EN'][ArgumentTypeException] = "Invalid argument type.";
|
||||||
$error['FR'][-2] = "Type inadmissible d'argument.";
|
$error['FR'][ArgumentTypeException] = "Type inadmissible d'argument.";
|
||||||
$error['DE'][-2] = "Unzulässige Argumentart.";
|
$error['DE'][ArgumentTypeException] = "Unzulässige Argumentart.";
|
||||||
|
|
||||||
define('ArgumentBoundsException', -3);
|
define('ArgumentBoundsException', -3);
|
||||||
$error['EN'][-3] = "Invalid argument range.";
|
$error['EN'][ArgumentBoundsException] = "Invalid argument range.";
|
||||||
$error['FR'][-3] = "Gamme inadmissible d'argument.";
|
$error['FR'][ArgumentBoundsException] = "Gamme inadmissible d'argument.";
|
||||||
$error['DE'][-3] = "Unzulässige Argumentstrecke.";
|
$error['DE'][ArgumentBoundsException] = "Unzulässige Argumentstrecke.";
|
||||||
|
|
||||||
define('MatrixDimensionException', -4);
|
define('MatrixDimensionException', -4);
|
||||||
$error['EN'][-4] = "Matrix dimensions are not equal.";
|
$error['EN'][MatrixDimensionException] = "Matrix dimensions are not equal.";
|
||||||
$error['FR'][-4] = "Les dimensions de Matrix ne sont pas égales.";
|
$error['FR'][MatrixDimensionException] = "Les dimensions de Matrix ne sont pas égales.";
|
||||||
$error['DE'][-4] = "Matrixmaße sind nicht gleich.";
|
$error['DE'][MatrixDimensionException] = "Matrixmaße sind nicht gleich.";
|
||||||
|
|
||||||
define('PrecisionLossException', -5);
|
define('PrecisionLossException', -5);
|
||||||
$error['EN'][-5] = "Significant precision loss detected.";
|
$error['EN'][PrecisionLossException] = "Significant precision loss detected.";
|
||||||
$error['FR'][-5] = "Perte significative de précision détectée.";
|
$error['FR'][PrecisionLossException] = "Perte significative de précision détectée.";
|
||||||
$error['DE'][-5] = "Bedeutender Präzision Verlust ermittelte.";
|
$error['DE'][PrecisionLossException] = "Bedeutender Präzision Verlust ermittelte.";
|
||||||
|
|
||||||
define('MatrixSPDException', -6);
|
define('MatrixSPDException', -6);
|
||||||
$error['EN'][-6] = "Can only perform operation on symmetric positive definite matrix.";
|
$error['EN'][MatrixSPDException] = "Can only perform operation on symmetric positive definite matrix.";
|
||||||
$error['FR'][-6] = "Perte significative de précision détectée.";
|
$error['FR'][MatrixSPDException] = "Perte significative de précision détectée.";
|
||||||
$error['DE'][-6] = "Bedeutender Präzision Verlust ermittelte.";
|
$error['DE'][MatrixSPDException] = "Bedeutender Präzision Verlust ermittelte.";
|
||||||
|
|
||||||
define('MatrixSingularException', -7);
|
define('MatrixSingularException', -7);
|
||||||
$error['EN'][-7] = "Can only perform operation on singular matrix.";
|
$error['EN'][MatrixSingularException] = "Can only perform operation on singular matrix.";
|
||||||
|
|
||||||
define('MatrixRankException', -8);
|
define('MatrixRankException', -8);
|
||||||
$error['EN'][-8] = "Can only perform operation on full-rank matrix.";
|
$error['EN'][MatrixRankException] = "Can only perform operation on full-rank matrix.";
|
||||||
|
|
||||||
define('ArrayLengthException', -9);
|
define('ArrayLengthException', -9);
|
||||||
$error['EN'][-9] = "Array length must be a multiple of m.";
|
$error['EN'][ArrayLengthException] = "Array length must be a multiple of m.";
|
||||||
|
|
||||||
define('RowLengthException', -10);
|
define('RowLengthException', -10);
|
||||||
$error['EN'][-10] = "All rows must have the same length.";
|
$error['EN'][RowLengthException] = "All rows must have the same length.";
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Custom error handler
|
* Custom error handler
|
||||||
* @param int $type Error type: {ERROR, WARNING, NOTICE}
|
* @param int $num Error number
|
||||||
* @param int $num Error number
|
*/
|
||||||
* @param string $file File in which the error occured
|
function JAMAError($errorNumber = null) {
|
||||||
* @param int $line Line on which the error occured
|
global $error;
|
||||||
*/
|
|
||||||
function JAMAError( $type = null, $num = null, $file = null, $line = null, $context = null ) {
|
|
||||||
global $error;
|
|
||||||
|
|
||||||
$lang = LANG;
|
if (isset($errorNumber)) {
|
||||||
if( isset($type) && isset($num) && isset($file) && isset($line) ) {
|
if (isset($error[JAMALANG][$errorNumber])) {
|
||||||
switch( $type ) {
|
return $error[JAMALANG][$errorNumber];
|
||||||
case ERROR:
|
} else {
|
||||||
echo '<div class="errror"><b>Error:</b> ' . $error[$lang][$num] . '<br />' . $file . ' @ L' . $line . '</div>';
|
return $error['EN'][$errorNumber];
|
||||||
die();
|
}
|
||||||
break;
|
} else {
|
||||||
|
return ("Invalid argument to JAMAError()");
|
||||||
case WARNING:
|
}
|
||||||
echo '<div class="warning"><b>Warning:</b> ' . $error[$lang][$num] . '<br />' . $file . ' @ L' . $line . '</div>';
|
|
||||||
break;
|
|
||||||
|
|
||||||
case NOTICE:
|
|
||||||
//echo '<div class="notice"><b>Notice:</b> ' . $error[$lang][$num] . '<br />' . $file . ' @ L' . $line . '</div>';
|
|
||||||
break;
|
|
||||||
|
|
||||||
case E_NOTICE:
|
|
||||||
//echo '<div class="errror"><b>Notice:</b> ' . $error[$lang][$num] . '<br />' . $file . ' @ L' . $line . '</div>';
|
|
||||||
break;
|
|
||||||
|
|
||||||
case E_STRICT:
|
|
||||||
break;
|
|
||||||
|
|
||||||
case E_WARNING:
|
|
||||||
break;
|
|
||||||
|
|
||||||
default:
|
|
||||||
echo "<div class=\"error\"><b>Unknown Error Type:</b> $type - $file @ L{$line}</div>";
|
|
||||||
die();
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
} else {
|
|
||||||
die( "Invalid arguments to JAMAError()" );
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// TODO MarkBaker
|
|
||||||
//set_error_handler('JAMAError');
|
|
||||||
//error_reporting(ERROR | WARNING);
|
|
||||||
|
|
||||||
|
|||||||
@@ -1,40 +1,43 @@
|
|||||||
<?php
|
<?php
|
||||||
/**
|
/**
|
||||||
* @package JAMA
|
* @package JAMA
|
||||||
*
|
*
|
||||||
* Pythagorean Theorem:
|
* Pythagorean Theorem:
|
||||||
*
|
*
|
||||||
* a = 3
|
* a = 3
|
||||||
* b = 4
|
* b = 4
|
||||||
* r = sqrt(square(a) + square(b))
|
* r = sqrt(square(a) + square(b))
|
||||||
* r = 5
|
* r = 5
|
||||||
*
|
*
|
||||||
* r = sqrt(a^2 + b^2) without under/overflow.
|
* r = sqrt(a^2 + b^2) without under/overflow.
|
||||||
*/
|
*/
|
||||||
function hypo($a, $b) {
|
function hypo($a, $b) {
|
||||||
if (abs($a) > abs($b)) {
|
if (abs($a) > abs($b)) {
|
||||||
$r = $b/$a;
|
$r = $b / $a;
|
||||||
$r = abs($a)* sqrt(1+$r*$r);
|
$r = abs($a) * sqrt(1 + $r * $r);
|
||||||
} else if ($b != 0) {
|
} elseif ($b != 0) {
|
||||||
$r = $a/$b;
|
$r = $a / $b;
|
||||||
$r = abs($b)*sqrt(1+$r*$r);
|
$r = abs($b) * sqrt(1 + $r * $r);
|
||||||
} else
|
} else {
|
||||||
$r = 0.0;
|
$r = 0.0;
|
||||||
return $r;
|
}
|
||||||
}
|
return $r;
|
||||||
|
} // function hypo()
|
||||||
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Mike Bommarito's version.
|
* Mike Bommarito's version.
|
||||||
* Compute n-dimensional hyotheneuse.
|
* Compute n-dimensional hyotheneuse.
|
||||||
*
|
*
|
||||||
function hypot() {
|
function hypot() {
|
||||||
$s = 0;
|
$s = 0;
|
||||||
foreach (func_get_args() as $d) {
|
foreach (func_get_args() as $d) {
|
||||||
if (is_numeric($d))
|
if (is_numeric($d)) {
|
||||||
$s += pow($d, 2);
|
$s += pow($d, 2);
|
||||||
else
|
} else {
|
||||||
trigger_error(ArgumentTypeException, ERROR);
|
throw new Exception(JAMAError(ArgumentTypeException));
|
||||||
}
|
}
|
||||||
return sqrt($s);
|
}
|
||||||
|
return sqrt($s);
|
||||||
}
|
}
|
||||||
*/
|
*/
|
||||||
|
|||||||
Reference in New Issue
Block a user