Passed
Push — master ( e83f7b...d953ef )
 1 = n, the LU decomposition is an m-by-n 9 * unit lower triangular matrix L, an n-by-n upper triangular matrix U, 10 * and a permutation vector piv of length m so that A(piv,:) = L*U. 11 * If m < n, then L is m-by-m and U is m-by-n. 12 * 13 * The LU decompostion with pivoting always exists, even if the matrix is 14 * singular, so the constructor will never fail. The primary use of the 15 * LU decomposition is in the solution of square systems of simultaneous 16 * linear equations. This will fail if isNonsingular() returns false. 17 * 18 * @author Paul Meagher 19 * @author Bartosz Matosiuk 20 * @author Michael Bommarito 21 * 22 * @version 1.1 23 * 24 * @license PHP v3.0 25 * 26 * Slightly changed to adapt the original code to PHP-ML library 27 * @date 2017/04/24 28 * 29 * @author Mustafa Karabulut 30 */ 31 32 namespace Phpml\Math\LinearAlgebra; 33 34 use Phpml\Exception\MatrixException; 35 use Phpml\Math\Matrix; 36 37 class LUDecomposition 38 { 39 /** 40 * Decomposition storage 41 * 42 * @var array 43 */ 44 private \$LU = []; 45 46 /** 47 * Row dimension. 48 * 49 * @var int 50 */ 51 private \$m; 52 53 /** 54 * Column dimension. 55 * 56 * @var int 57 */ 58 private \$n; 59 60 /** 61 * Pivot sign. 62 * 63 * @var int 64 */ 65 private \$pivsign; 66 67 /** 68 * Internal storage of pivot vector. 69 * 70 * @var array 71 */ 72 private \$piv = []; 73 74 /** 75 * Constructs Structure to access L, U and piv. 76 * 77 * @param Matrix \$A Rectangular matrix 78 * 79 * @throws MatrixException 80 */ 81 public function __construct(Matrix \$A) 82 { 83 if (\$A->getRows() != \$A->getColumns()) { 84 throw MatrixException::notSquareMatrix(); 85 } 86 87 // Use a "left-looking", dot-product, Crout/Doolittle algorithm. 88 \$this->LU = \$A->toArray(); 89 \$this->m = \$A->getRows(); 90 \$this->n = \$A->getColumns(); 91 for (\$i = 0; \$i < \$this->m; ++\$i) { 92 \$this->piv[\$i] = \$i; 93 } 94 95 \$this->pivsign = 1; 96 \$LUcolj = []; 97 98 // Outer loop. 99 for (\$j = 0; \$j < \$this->n; ++\$j) { 100 // Make a copy of the j-th column to localize references. 101 View Code Duplication for (\$i = 0; \$i < \$this->m; ++\$i) { 0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 102 \$LUcolj[\$i] = &\$this->LU[\$i][\$j]; 103 } 104 105 // Apply previous transformations. 106 for (\$i = 0; \$i < \$this->m; ++\$i) { 107 \$LUrowi = \$this->LU[\$i]; 108 // Most of the time is spent in the following dot product. 109 \$kmax = min(\$i, \$j); 110 \$s = 0.0; 111 for (\$k = 0; \$k < \$kmax; ++\$k) { 112 \$s += \$LUrowi[\$k] * \$LUcolj[\$k]; 113 } 114 115 \$LUrowi[\$j] = \$LUcolj[\$i] -= \$s; 116 } 117 118 // Find pivot and exchange if necessary. 119 \$p = \$j; 120 for (\$i = \$j + 1; \$i < \$this->m; ++\$i) { 121 if (abs(\$LUcolj[\$i]) > abs(\$LUcolj[\$p])) { 122 \$p = \$i; 123 } 124 } 125 126 if (\$p != \$j) { 127 View Code Duplication for (\$k = 0; \$k < \$this->n; ++\$k) { 0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 128 \$t = \$this->LU[\$p][\$k]; 129 \$this->LU[\$p][\$k] = \$this->LU[\$j][\$k]; 130 \$this->LU[\$j][\$k] = \$t; 131 } 132 133 \$k = \$this->piv[\$p]; 134 \$this->piv[\$p] = \$this->piv[\$j]; 135 \$this->piv[\$j] = \$k; 136 \$this->pivsign = \$this->pivsign * -1; 137 } 138 139 // Compute multipliers. 140 if ((\$j < \$this->m) && (\$this->LU[\$j][\$j] != 0.0)) { 141 View Code Duplication for (\$i = \$j + 1; \$i < \$this->m; ++\$i) { 0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 142 \$this->LU[\$i][\$j] /= \$this->LU[\$j][\$j]; 143 } 144 } 145 } 146 } 147 148 /** 149 * Get lower triangular factor. 150 * 151 * @return Matrix Lower triangular factor 152 */ 153 public function getL(): Matrix 154 { 155 \$L = []; 156 for (\$i = 0; \$i < \$this->m; ++\$i) { 157 for (\$j = 0; \$j < \$this->n; ++\$j) { 158 if (\$i > \$j) { 159 \$L[\$i][\$j] = \$this->LU[\$i][\$j]; 160 } elseif (\$i == \$j) { 161 \$L[\$i][\$j] = 1.0; 162 } else { 163 \$L[\$i][\$j] = 0.0; 164 } 165 } 166 } 167 168 return new Matrix(\$L); 169 } 170 171 /** 172 * Get upper triangular factor. 173 * 174 * @return Matrix Upper triangular factor 175 */ 176 public function getU(): Matrix 177 { 178 \$U = []; 179 for (\$i = 0; \$i < \$this->n; ++\$i) { 180 for (\$j = 0; \$j < \$this->n; ++\$j) { 181 if (\$i <= \$j) { 182 \$U[\$i][\$j] = \$this->LU[\$i][\$j]; 183 } else { 184 \$U[\$i][\$j] = 0.0; 185 } 186 } 187 } 188 189 return new Matrix(\$U); 190 } 191 192 /** 193 * Return pivot permutation vector. 194 * 195 * @return array Pivot vector 196 */ 197 public function getPivot(): array 198 { 199 return \$this->piv; 200 } 201 202 /** 203 * Alias for getPivot 204 * 205 * @see getPivot 206 */ 207 public function getDoublePivot() 208 { 209 return \$this->getPivot(); 210 } 211 212 /** 213 * Is the matrix nonsingular? 214 * 215 * @return bool true if U, and hence A, is nonsingular. 216 */ 217 public function isNonsingular(): bool 218 { 219 View Code Duplication for (\$j = 0; \$j < \$this->n; ++\$j) { 0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 220 if (\$this->LU[\$j][\$j] == 0) { 221 return false; 222 } 223 } 224 225 return true; 226 } 227 228 /** 229 * Count determinants 230 * 231 * @return float|int d matrix determinant 232 * 233 * @throws MatrixException 234 */ 235 public function det() 236 { 237 if (\$this->m !== \$this->n) { 238 throw MatrixException::notSquareMatrix(); 239 } 240 241 \$d = \$this->pivsign; 242 View Code Duplication for (\$j = 0; \$j < \$this->n; ++\$j) { 0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 243 \$d *= \$this->LU[\$j][\$j]; 244 } 245 246 return \$d; 247 } 248 249 /** 250 * Solve A*X = B 251 * 252 * @param Matrix \$B A Matrix with as many rows as A and any number of columns. 253 * 254 * @return array X so that L*U*X = B(piv,:) 255 * 256 * @throws MatrixException 257 */ 258 public function solve(Matrix \$B): array 259 { 260 if (\$B->getRows() != \$this->m) { 261 throw MatrixException::notSquareMatrix(); 262 } 263 264 if (!\$this->isNonsingular()) { 265 throw MatrixException::singularMatrix(); 266 } 267 268 // Copy right hand side with pivoting 269 \$nx = \$B->getColumns(); 270 \$X = \$this->getSubMatrix(\$B->toArray(), \$this->piv, 0, \$nx - 1); 271 // Solve L*Y = B(piv,:) 272 for (\$k = 0; \$k < \$this->n; ++\$k) { 273 View Code Duplication for (\$i = \$k + 1; \$i < \$this->n; ++\$i) { 0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 274 for (\$j = 0; \$j < \$nx; ++\$j) { 275 \$X[\$i][\$j] -= \$X[\$k][\$j] * \$this->LU[\$i][\$k]; 276 } 277 } 278 } 279 280 // Solve U*X = Y; 281 for (\$k = \$this->n - 1; \$k >= 0; --\$k) { 282 for (\$j = 0; \$j < \$nx; ++\$j) { 283 \$X[\$k][\$j] /= \$this->LU[\$k][\$k]; 284 } 285 286 View Code Duplication for (\$i = 0; \$i < \$k; ++\$i) { 0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 287 for (\$j = 0; \$j < \$nx; ++\$j) { 288 \$X[\$i][\$j] -= \$X[\$k][\$j] * \$this->LU[\$i][\$k]; 289 } 290 } 291 } 292 293 return \$X; 294 } 295 296 protected function getSubMatrix(array \$matrix, array \$RL, int \$j0, int \$jF): array 297 { 298 \$m = count(\$RL); 299 \$n = \$jF - \$j0; 300 \$R = array_fill(0, \$m, array_fill(0, \$n + 1, 0.0)); 301 302 for (\$i = 0; \$i < \$m; ++\$i) { 303 for (\$j = \$j0; \$j <= \$jF; ++\$j) { 304 \$R[\$i][\$j - \$j0] = \$matrix[\$RL[\$i]][\$j]; 305 } 306 } 307 308 return \$R; 309 } 310 } 311