Subversion Repositories eFlore/Applications.cel

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2388 jpm 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
16
 *	@license PHP v3.0
17
 *	@version 1.1
18
 */
19
class 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 PHPExcel_Calculation_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 PHPExcel_Calculation_Exception(self::MatrixRankException);
228
			}
229
		} else {
230
			throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
231
		}
232
	}	//	function solve()
233
 
234
}	//	PHPExcel_Shared_JAMA_class QRDecomposition