#### /PHPExcel/Shared/JAMA/QRDecomposition.php

https://bitbucket.org/nfredricks/wp-employee-time
PHP | 234 lines | 134 code | 18 blank | 82 comment | 43 complexity | 95f548ec797ed8f8c06d6d134d33ebbd MD5 | raw file
```  1<?php
2/**
3 *	@package JAMA
4 *
5 *	For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
6 *	orthogonal matrix Q and an n-by-n upper triangular matrix R so that
7 *	A = Q*R.
8 *
9 *	The QR decompostion always exists, even if the matrix does not have
10 *	full rank, so the constructor will never fail.  The primary use of the
11 *	QR decomposition is in the least squares solution of nonsquare systems
12 *	of simultaneous linear equations.  This will fail if isFullRank()
13 *	returns false.
14 *
15 *	@author  Paul Meagher
17 *	@version 1.1
18 */
19class PHPExcel_Shared_JAMA_QRDecomposition {
20
21	const MatrixRankException	= "Can only perform operation on full-rank matrix.";
22
23	/**
24	 *	Array for internal storage of decomposition.
25	 *	@var array
26	 */
27	private \$QR = array();
28
29	/**
30	 *	Row dimension.
31	 *	@var integer
32	 */
33	private \$m;
34
35	/**
36	*	Column dimension.
37	*	@var integer
38	*/
39	private \$n;
40
41	/**
42	 *	Array for internal storage of diagonal of R.
43	 *	@var  array
44	 */
45	private \$Rdiag = array();
46
47
48	/**
49	 *	QR Decomposition computed by Householder reflections.
50	 *
51	 *	@param matrix \$A Rectangular matrix
52	 *	@return Structure to access R and the Householder vectors and compute Q.
53	 */
54	public function __construct(\$A) {
55		if(\$A instanceof PHPExcel_Shared_JAMA_Matrix) {
56			// Initialize.
57			\$this->QR = \$A->getArrayCopy();
58			\$this->m  = \$A->getRowDimension();
59			\$this->n  = \$A->getColumnDimension();
60			// Main loop.
61			for (\$k = 0; \$k < \$this->n; ++\$k) {
62				// Compute 2-norm of k-th column without under/overflow.
63				\$nrm = 0.0;
64				for (\$i = \$k; \$i < \$this->m; ++\$i) {
65					\$nrm = hypo(\$nrm, \$this->QR[\$i][\$k]);
66				}
67				if (\$nrm != 0.0) {
68					// Form k-th Householder vector.
69					if (\$this->QR[\$k][\$k] < 0) {
70						\$nrm = -\$nrm;
71					}
72					for (\$i = \$k; \$i < \$this->m; ++\$i) {
73						\$this->QR[\$i][\$k] /= \$nrm;
74					}
75					\$this->QR[\$k][\$k] += 1.0;
76					// Apply transformation to remaining columns.
77					for (\$j = \$k+1; \$j < \$this->n; ++\$j) {
78						\$s = 0.0;
79						for (\$i = \$k; \$i < \$this->m; ++\$i) {
80							\$s += \$this->QR[\$i][\$k] * \$this->QR[\$i][\$j];
81						}
82						\$s = -\$s/\$this->QR[\$k][\$k];
83						for (\$i = \$k; \$i < \$this->m; ++\$i) {
84							\$this->QR[\$i][\$j] += \$s * \$this->QR[\$i][\$k];
85						}
86					}
87				}
88				\$this->Rdiag[\$k] = -\$nrm;
89			}
90		} else {
91			throw new Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
92		}
93	}	//	function __construct()
94
95
96	/**
97	 *	Is the matrix full rank?
98	 *
99	 *	@return boolean true if R, and hence A, has full rank, else false.
100	 */
101	public function isFullRank() {
102		for (\$j = 0; \$j < \$this->n; ++\$j) {
103			if (\$this->Rdiag[\$j] == 0) {
104				return false;
105			}
106		}
107		return true;
108	}	//	function isFullRank()
109
110
111	/**
112	 *	Return the Householder vectors
113	 *
114	 *	@return Matrix Lower trapezoidal matrix whose columns define the reflections
115	 */
116	public function getH() {
117		for (\$i = 0; \$i < \$this->m; ++\$i) {
118			for (\$j = 0; \$j < \$this->n; ++\$j) {
119				if (\$i >= \$j) {
120					\$H[\$i][\$j] = \$this->QR[\$i][\$j];
121				} else {
122					\$H[\$i][\$j] = 0.0;
123				}
124			}
125		}
126		return new PHPExcel_Shared_JAMA_Matrix(\$H);
127	}	//	function getH()
128
129
130	/**
131	 *	Return the upper triangular factor
132	 *
133	 *	@return Matrix upper triangular factor
134	 */
135	public function getR() {
136		for (\$i = 0; \$i < \$this->n; ++\$i) {
137			for (\$j = 0; \$j < \$this->n; ++\$j) {
138				if (\$i < \$j) {
139					\$R[\$i][\$j] = \$this->QR[\$i][\$j];
140				} elseif (\$i == \$j) {
141					\$R[\$i][\$j] = \$this->Rdiag[\$i];
142				} else {
143					\$R[\$i][\$j] = 0.0;
144				}
145			}
146		}
147		return new PHPExcel_Shared_JAMA_Matrix(\$R);
148	}	//	function getR()
149
150
151	/**
152	 *	Generate and return the (economy-sized) orthogonal factor
153	 *
154	 *	@return Matrix orthogonal factor
155	 */
156	public function getQ() {
157		for (\$k = \$this->n-1; \$k >= 0; --\$k) {
158			for (\$i = 0; \$i < \$this->m; ++\$i) {
159				\$Q[\$i][\$k] = 0.0;
160			}
161			\$Q[\$k][\$k] = 1.0;
162			for (\$j = \$k; \$j < \$this->n; ++\$j) {
163				if (\$this->QR[\$k][\$k] != 0) {
164					\$s = 0.0;
165					for (\$i = \$k; \$i < \$this->m; ++\$i) {
166						\$s += \$this->QR[\$i][\$k] * \$Q[\$i][\$j];
167					}
168					\$s = -\$s/\$this->QR[\$k][\$k];
169					for (\$i = \$k; \$i < \$this->m; ++\$i) {
170						\$Q[\$i][\$j] += \$s * \$this->QR[\$i][\$k];
171					}
172				}
173			}
174		}
175		/*
176		for(\$i = 0; \$i < count(\$Q); ++\$i) {
177			for(\$j = 0; \$j < count(\$Q); ++\$j) {
178				if(! isset(\$Q[\$i][\$j]) ) {
179					\$Q[\$i][\$j] = 0;
180				}
181			}
182		}
183		*/
184		return new PHPExcel_Shared_JAMA_Matrix(\$Q);
185	}	//	function getQ()
186
187
188	/**
189	 *	Least squares solution of A*X = B
190	 *
191	 *	@param Matrix \$B A Matrix with as many rows as A and any number of columns.
192	 *	@return Matrix Matrix that minimizes the two norm of Q*R*X-B.
193	 */
194	public function solve(\$B) {
195		if (\$B->getRowDimension() == \$this->m) {
196			if (\$this->isFullRank()) {
197				// Copy right hand side
198				\$nx = \$B->getColumnDimension();
199				\$X  = \$B->getArrayCopy();
200				// Compute Y = transpose(Q)*B
201				for (\$k = 0; \$k < \$this->n; ++\$k) {
202					for (\$j = 0; \$j < \$nx; ++\$j) {
203						\$s = 0.0;
204						for (\$i = \$k; \$i < \$this->m; ++\$i) {
205							\$s += \$this->QR[\$i][\$k] * \$X[\$i][\$j];
206						}
207						\$s = -\$s/\$this->QR[\$k][\$k];
208						for (\$i = \$k; \$i < \$this->m; ++\$i) {
209							\$X[\$i][\$j] += \$s * \$this->QR[\$i][\$k];
210						}
211					}
212				}
213				// Solve R*X = Y;
214				for (\$k = \$this->n-1; \$k >= 0; --\$k) {
215					for (\$j = 0; \$j < \$nx; ++\$j) {
216						\$X[\$k][\$j] /= \$this->Rdiag[\$k];
217					}
218					for (\$i = 0; \$i < \$k; ++\$i) {
219						for (\$j = 0; \$j < \$nx; ++\$j) {
220							\$X[\$i][\$j] -= \$X[\$k][\$j]* \$this->QR[\$i][\$k];
221						}
222					}
223				}
224				\$X = new PHPExcel_Shared_JAMA_Matrix(\$X);
225				return (\$X->getMatrix(0, \$this->n-1, 0, \$nx));
226			} else {
227				throw new Exception(self::MatrixRankException);
228			}
229		} else {
230			throw new Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
231		}
232	}	//	function solve()
233
234}	//	PHPExcel_Shared_JAMA_class QRDecomposition
```