Passed
Push — master ( e83f7b...d953ef )
by Arkadiusz
03:28
created

src/Phpml/Math/LinearAlgebra/LUDecomposition.php (7 issues)

Upgrade to new PHP Analysis Engine

These results are based on our legacy PHP analysis, consider migrating to our new PHP analysis engine instead. Learn more

1
<?php
2
3
declare(strict_types=1);
4
5
/**
6
 * @package JAMA
7
 *
8
 * For an m-by-n matrix A with m >= 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
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
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
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
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
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
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
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