diff options
Diffstat (limited to 'admin/survey/excel/PHPExcel/Shared/JAMA/QRDecomposition.php')
| -rw-r--r-- | admin/survey/excel/PHPExcel/Shared/JAMA/QRDecomposition.php | 234 |
1 files changed, 234 insertions, 0 deletions
diff --git a/admin/survey/excel/PHPExcel/Shared/JAMA/QRDecomposition.php b/admin/survey/excel/PHPExcel/Shared/JAMA/QRDecomposition.php new file mode 100644 index 0000000..4492fa5 --- /dev/null +++ b/admin/survey/excel/PHPExcel/Shared/JAMA/QRDecomposition.php @@ -0,0 +1,234 @@ +<?php
+/**
+ * @package JAMA
+ *
+ * For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
+ * orthogonal matrix Q and an n-by-n upper triangular matrix R so that
+ * A = Q*R.
+ *
+ * The QR decompostion always exists, even if the matrix does not have
+ * full rank, so the constructor will never fail. The primary use of the
+ * QR decomposition is in the least squares solution of nonsquare systems
+ * of simultaneous linear equations. This will fail if isFullRank()
+ * returns false.
+ *
+ * @author Paul Meagher
+ * @license PHP v3.0
+ * @version 1.1
+ */
+class PHPExcel_Shared_JAMA_QRDecomposition {
+
+ const MatrixRankException = "Can only perform operation on full-rank matrix.";
+
+ /**
+ * Array for internal storage of decomposition.
+ * @var array
+ */
+ private $QR = array();
+
+ /**
+ * Row dimension.
+ * @var integer
+ */
+ private $m;
+
+ /**
+ * Column dimension.
+ * @var integer
+ */
+ private $n;
+
+ /**
+ * Array for internal storage of diagonal of R.
+ * @var array
+ */
+ private $Rdiag = array();
+
+
+ /**
+ * QR Decomposition computed by Householder reflections.
+ *
+ * @param matrix $A Rectangular matrix
+ * @return Structure to access R and the Householder vectors and compute Q.
+ */
+ public function __construct($A) {
+ if($A instanceof PHPExcel_Shared_JAMA_Matrix) {
+ // Initialize.
+ $this->QR = $A->getArrayCopy();
+ $this->m = $A->getRowDimension();
+ $this->n = $A->getColumnDimension();
+ // Main loop.
+ for ($k = 0; $k < $this->n; ++$k) {
+ // Compute 2-norm of k-th column without under/overflow.
+ $nrm = 0.0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $nrm = hypo($nrm, $this->QR[$i][$k]);
+ }
+ if ($nrm != 0.0) {
+ // Form k-th Householder vector.
+ if ($this->QR[$k][$k] < 0) {
+ $nrm = -$nrm;
+ }
+ for ($i = $k; $i < $this->m; ++$i) {
+ $this->QR[$i][$k] /= $nrm;
+ }
+ $this->QR[$k][$k] += 1.0;
+ // Apply transformation to remaining columns.
+ for ($j = $k+1; $j < $this->n; ++$j) {
+ $s = 0.0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $s += $this->QR[$i][$k] * $this->QR[$i][$j];
+ }
+ $s = -$s/$this->QR[$k][$k];
+ for ($i = $k; $i < $this->m; ++$i) {
+ $this->QR[$i][$j] += $s * $this->QR[$i][$k];
+ }
+ }
+ }
+ $this->Rdiag[$k] = -$nrm;
+ }
+ } else {
+ throw new Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
+ }
+ } // function __construct()
+
+
+ /**
+ * Is the matrix full rank?
+ *
+ * @return boolean true if R, and hence A, has full rank, else false.
+ */
+ public function isFullRank() {
+ for ($j = 0; $j < $this->n; ++$j) {
+ if ($this->Rdiag[$j] == 0) {
+ return false;
+ }
+ }
+ return true;
+ } // function isFullRank()
+
+
+ /**
+ * Return the Householder vectors
+ *
+ * @return Matrix Lower trapezoidal matrix whose columns define the reflections
+ */
+ public function getH() {
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ if ($i >= $j) {
+ $H[$i][$j] = $this->QR[$i][$j];
+ } else {
+ $H[$i][$j] = 0.0;
+ }
+ }
+ }
+ return new PHPExcel_Shared_JAMA_Matrix($H);
+ } // function getH()
+
+
+ /**
+ * Return the upper triangular factor
+ *
+ * @return Matrix upper triangular factor
+ */
+ public function getR() {
+ for ($i = 0; $i < $this->n; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ if ($i < $j) {
+ $R[$i][$j] = $this->QR[$i][$j];
+ } elseif ($i == $j) {
+ $R[$i][$j] = $this->Rdiag[$i];
+ } else {
+ $R[$i][$j] = 0.0;
+ }
+ }
+ }
+ return new PHPExcel_Shared_JAMA_Matrix($R);
+ } // function getR()
+
+
+ /**
+ * Generate and return the (economy-sized) orthogonal factor
+ *
+ * @return Matrix orthogonal factor
+ */
+ public function getQ() {
+ for ($k = $this->n-1; $k >= 0; --$k) {
+ for ($i = 0; $i < $this->m; ++$i) {
+ $Q[$i][$k] = 0.0;
+ }
+ $Q[$k][$k] = 1.0;
+ for ($j = $k; $j < $this->n; ++$j) {
+ if ($this->QR[$k][$k] != 0) {
+ $s = 0.0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $s += $this->QR[$i][$k] * $Q[$i][$j];
+ }
+ $s = -$s/$this->QR[$k][$k];
+ for ($i = $k; $i < $this->m; ++$i) {
+ $Q[$i][$j] += $s * $this->QR[$i][$k];
+ }
+ }
+ }
+ }
+ /*
+ for($i = 0; $i < count($Q); ++$i) {
+ for($j = 0; $j < count($Q); ++$j) {
+ if(! isset($Q[$i][$j]) ) {
+ $Q[$i][$j] = 0;
+ }
+ }
+ }
+ */
+ return new PHPExcel_Shared_JAMA_Matrix($Q);
+ } // function getQ()
+
+
+ /**
+ * Least squares solution of A*X = B
+ *
+ * @param Matrix $B A Matrix with as many rows as A and any number of columns.
+ * @return Matrix Matrix that minimizes the two norm of Q*R*X-B.
+ */
+ public function solve($B) {
+ if ($B->getRowDimension() == $this->m) {
+ if ($this->isFullRank()) {
+ // Copy right hand side
+ $nx = $B->getColumnDimension();
+ $X = $B->getArrayCopy();
+ // Compute Y = transpose(Q)*B
+ for ($k = 0; $k < $this->n; ++$k) {
+ for ($j = 0; $j < $nx; ++$j) {
+ $s = 0.0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $s += $this->QR[$i][$k] * $X[$i][$j];
+ }
+ $s = -$s/$this->QR[$k][$k];
+ for ($i = $k; $i < $this->m; ++$i) {
+ $X[$i][$j] += $s * $this->QR[$i][$k];
+ }
+ }
+ }
+ // Solve R*X = Y;
+ for ($k = $this->n-1; $k >= 0; --$k) {
+ for ($j = 0; $j < $nx; ++$j) {
+ $X[$k][$j] /= $this->Rdiag[$k];
+ }
+ for ($i = 0; $i < $k; ++$i) {
+ for ($j = 0; $j < $nx; ++$j) {
+ $X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];
+ }
+ }
+ }
+ $X = new PHPExcel_Shared_JAMA_Matrix($X);
+ return ($X->getMatrix(0, $this->n-1, 0, $nx));
+ } else {
+ throw new Exception(self::MatrixRankException);
+ }
+ } else {
+ throw new Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
+ }
+ } // function solve()
+
+} // PHPExcel_Shared_JAMA_class QRDecomposition
|