Subversion Repositories eFlore/Applications.cel

<?php

/**

 *	@package JAMA

 *

 *	For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n

 *	orthogonal matrix Q and an n-by-n upper triangular matrix R so that

 *	A = Q*R.

 *

 *	The QR decompostion always exists, even if the matrix does not have

 *	full rank, so the constructor will never fail.  The primary use of the

 *	QR decomposition is in the least squares solution of nonsquare systems

 *	of simultaneous linear equations.  This will fail if isFullRank()

 *	returns false.

 *

 *	@author  Paul Meagher

 *	@license PHP v3.0

 *	@version 1.1

 */

class PHPExcel_Shared_JAMA_QRDecomposition {

	const MatrixRankException	= "Can only perform operation on full-rank matrix.";

	/**

	 *	Array for internal storage of decomposition.

	 *	@var array

	 */

	private $QR = array();

	/**

	 *	Row dimension.

	 *	@var integer

	 */

	private $m;

	/**

	*	Column dimension.

	*	@var integer

	*/

	private $n;

	/**

	 *	Array for internal storage of diagonal of R.

	 *	@var  array

	 */

	private $Rdiag = array();

	/**

	 *	QR Decomposition computed by Householder reflections.

	 *

	 *	@param matrix $A Rectangular matrix

	 *	@return Structure to access R and the Householder vectors and compute Q.

	 */

	public function __construct($A) {

		if($A instanceof PHPExcel_Shared_JAMA_Matrix) {

			// Initialize.

			$this->QR = $A->getArrayCopy();

			$this->m  = $A->getRowDimension();

			$this->n  = $A->getColumnDimension();

			// Main loop.

			for ($k = 0; $k < $this->n; ++$k) {

				// Compute 2-norm of k-th column without under/overflow.

				$nrm = 0.0;

				for ($i = $k; $i < $this->m; ++$i) {

					$nrm = hypo($nrm, $this->QR[$i][$k]);

				}

				if ($nrm != 0.0) {

					// Form k-th Householder vector.

					if ($this->QR[$k][$k] < 0) {

						$nrm = -$nrm;

					}

					for ($i = $k; $i < $this->m; ++$i) {

						$this->QR[$i][$k] /= $nrm;

					}

					$this->QR[$k][$k] += 1.0;

					// Apply transformation to remaining columns.

					for ($j = $k+1; $j < $this->n; ++$j) {

						$s = 0.0;

						for ($i = $k; $i < $this->m; ++$i) {

							$s += $this->QR[$i][$k] * $this->QR[$i][$j];

						}

						$s = -$s/$this->QR[$k][$k];

						for ($i = $k; $i < $this->m; ++$i) {

							$this->QR[$i][$j] += $s * $this->QR[$i][$k];

						}

					}

				}

				$this->Rdiag[$k] = -$nrm;

			}

		} else {

			throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);

		}

	}	//	function __construct()

	/**

	 *	Is the matrix full rank?

	 *

	 *	@return boolean true if R, and hence A, has full rank, else false.

	 */

	public function isFullRank() {

		for ($j = 0; $j < $this->n; ++$j) {

			if ($this->Rdiag[$j] == 0) {

				return false;

			}

		}

		return true;

	}	//	function isFullRank()

	/**

	 *	Return the Householder vectors

	 *

	 *	@return Matrix Lower trapezoidal matrix whose columns define the reflections

	 */

	public function getH() {

		for ($i = 0; $i < $this->m; ++$i) {

			for ($j = 0; $j < $this->n; ++$j) {

				if ($i >= $j) {

					$H[$i][$j] = $this->QR[$i][$j];

				} else {

					$H[$i][$j] = 0.0;

				}

			}

		}

		return new PHPExcel_Shared_JAMA_Matrix($H);

	}	//	function getH()

	/**

	 *	Return the upper triangular factor

	 *

	 *	@return Matrix upper triangular factor

	 */

	public function getR() {

		for ($i = 0; $i < $this->n; ++$i) {

			for ($j = 0; $j < $this->n; ++$j) {

				if ($i < $j) {

					$R[$i][$j] = $this->QR[$i][$j];

				} elseif ($i == $j) {

					$R[$i][$j] = $this->Rdiag[$i];

				} else {

					$R[$i][$j] = 0.0;

				}

			}

		}

		return new PHPExcel_Shared_JAMA_Matrix($R);

	}	//	function getR()

	/**

	 *	Generate and return the (economy-sized) orthogonal factor

	 *

	 *	@return Matrix orthogonal factor

	 */

	public function getQ() {

		for ($k = $this->n-1; $k >= 0; --$k) {

			for ($i = 0; $i < $this->m; ++$i) {

				$Q[$i][$k] = 0.0;

			}

			$Q[$k][$k] = 1.0;

			for ($j = $k; $j < $this->n; ++$j) {

				if ($this->QR[$k][$k] != 0) {

					$s = 0.0;

					for ($i = $k; $i < $this->m; ++$i) {

						$s += $this->QR[$i][$k] * $Q[$i][$j];

					}

					$s = -$s/$this->QR[$k][$k];

					for ($i = $k; $i < $this->m; ++$i) {

						$Q[$i][$j] += $s * $this->QR[$i][$k];

					}

				}

			}

		}

		/*

		for($i = 0; $i < count($Q); ++$i) {

			for($j = 0; $j < count($Q); ++$j) {

				if(! isset($Q[$i][$j]) ) {

					$Q[$i][$j] = 0;

				}

			}

		}

		*/

		return new PHPExcel_Shared_JAMA_Matrix($Q);

	}	//	function getQ()

	/**

	 *	Least squares solution of A*X = B

	 *

	 *	@param Matrix $B A Matrix with as many rows as A and any number of columns.

	 *	@return Matrix Matrix that minimizes the two norm of Q*R*X-B.

	 */

	public function solve($B) {

		if ($B->getRowDimension() == $this->m) {

			if ($this->isFullRank()) {

				// Copy right hand side

				$nx = $B->getColumnDimension();

				$X  = $B->getArrayCopy();

				// Compute Y = transpose(Q)*B

				for ($k = 0; $k < $this->n; ++$k) {

					for ($j = 0; $j < $nx; ++$j) {

						$s = 0.0;

						for ($i = $k; $i < $this->m; ++$i) {

							$s += $this->QR[$i][$k] * $X[$i][$j];

						}

						$s = -$s/$this->QR[$k][$k];

						for ($i = $k; $i < $this->m; ++$i) {

							$X[$i][$j] += $s * $this->QR[$i][$k];

						}

					}

				}

				// Solve R*X = Y;

				for ($k = $this->n-1; $k >= 0; --$k) {

					for ($j = 0; $j < $nx; ++$j) {

						$X[$k][$j] /= $this->Rdiag[$k];

					}

					for ($i = 0; $i < $k; ++$i) {

						for ($j = 0; $j < $nx; ++$j) {

							$X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];

						}

					}

				}

				$X = new PHPExcel_Shared_JAMA_Matrix($X);

				return ($X->getMatrix(0, $this->n-1, 0, $nx));

			} else {

				throw new PHPExcel_Calculation_Exception(self::MatrixRankException);

			}

		} else {

			throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);

		}

	}	//	function solve()

}	//	PHPExcel_Shared_JAMA_class QRDecomposition