Passed
Push — master ( 12b8b1...5b373f )
by Arkadiusz
04:53
created

LUDecomposition::__construct()   C

Complexity

Conditions 14
Paths 147

Size

Total Lines 59
Code Lines 36

Duplication

Lines 11
Ratio 18.64 %

Importance

Changes 0
Metric Value
dl 11
loc 59
rs 5.8008
c 0
b 0
f 0
cc 14
eloc 36
nc 147
nop 1

How to fix   Long Method    Complexity   

Long Method

Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.

For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.

Commonly applied refactorings include:

1
<?php declare(strict_types=1);
2
/**
3
 *	@package JAMA
4
 *
5
 *	For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
6
 *	unit lower triangular matrix L, an n-by-n upper triangular matrix U,
7
 *	and a permutation vector piv of length m so that A(piv,:) = L*U.
8
 *	If m < n, then L is m-by-m and U is m-by-n.
9
 *
10
 *	The LU decompostion with pivoting always exists, even if the matrix is
11
 *	singular, so the constructor will never fail. The primary use of the
12
 *	LU decomposition is in the solution of square systems of simultaneous
13
 *	linear equations. This will fail if isNonsingular() returns false.
14
 *
15
 *	@author Paul Meagher
16
 *	@author Bartosz Matosiuk
17
 *	@author Michael Bommarito
18
 *	@version 1.1
19
 *	@license PHP v3.0
20
 *
21
 *  Slightly changed to adapt the original code to PHP-ML library
22
 *  @date 2017/04/24
23
 *  @author Mustafa Karabulut
24
 */
25
26
namespace Phpml\Math\LinearAlgebra;
27
28
use Phpml\Math\Matrix;
29
use Phpml\Exception\MatrixException;
30
31
class LUDecomposition
32
{
33
    /**
34
     *	Decomposition storage
35
     *	@var array
36
     */
37
    private $LU = [];
38
39
    /**
40
     *	Row dimension.
41
     *	@var int
42
     */
43
    private $m;
44
45
    /**
46
     *	Column dimension.
47
     *	@var int
48
     */
49
    private $n;
50
51
    /**
52
     *	Pivot sign.
53
     *	@var int
54
     */
55
    private $pivsign;
56
57
    /**
58
     *	Internal storage of pivot vector.
59
     *	@var array
60
     */
61
    private $piv = [];
62
63
64
    /**
65
     *	LU Decomposition constructor.
66
     *
67
     *	@param $A Rectangular matrix
68
     *	@return Structure to access L, U and piv.
0 ignored issues
show
Comprehensibility Best Practice introduced by
Adding a @return annotation to constructors is generally not recommended as a constructor does not have a meaningful return value.

Adding a @return annotation to a constructor is not recommended, since a constructor does not have a meaningful return value.

Please refer to the PHP core documentation on constructors.

Loading history...
69
     */
70
    public function __construct(Matrix $A)
71
    {
72
        if ($A->getRows() != $A->getColumns()) {
73
            throw MatrixException::notSquareMatrix();
74
        }
75
76
        // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
77
        $this->LU = $A->toArray();
78
        $this->m  = $A->getRows();
79
        $this->n  = $A->getColumns();
80
        for ($i = 0; $i < $this->m; ++$i) {
81
            $this->piv[$i] = $i;
82
        }
83
        $this->pivsign = 1;
84
        $LUrowi = $LUcolj = [];
0 ignored issues
show
Unused Code introduced by
$LUrowi is not used, you could remove the assignment.

This check looks for variable assignements that are either overwritten by other assignments or where the variable is not used subsequently.

$myVar = 'Value';
$higher = false;

if (rand(1, 6) > 3) {
    $higher = true;
} else {
    $higher = false;
}

Both the $myVar assignment in line 1 and the $higher assignment in line 2 are dead. The first because $myVar is never used and the second because $higher is always overwritten for every possible time line.

Loading history...
85
86
        // Outer loop.
87
        for ($j = 0; $j < $this->n; ++$j) {
88
            // Make a copy of the j-th column to localize references.
89 View Code Duplication
            for ($i = 0; $i < $this->m; ++$i) {
0 ignored issues
show
Duplication introduced 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...
90
                $LUcolj[$i] = &$this->LU[$i][$j];
91
            }
92
            // Apply previous transformations.
93
            for ($i = 0; $i < $this->m; ++$i) {
94
                $LUrowi = $this->LU[$i];
95
                // Most of the time is spent in the following dot product.
96
                $kmax = min($i, $j);
97
                $s = 0.0;
98
                for ($k = 0; $k < $kmax; ++$k) {
99
                    $s += $LUrowi[$k] * $LUcolj[$k];
100
                }
101
                $LUrowi[$j] = $LUcolj[$i] -= $s;
102
            }
103
            // Find pivot and exchange if necessary.
104
            $p = $j;
105
            for ($i = $j+1; $i < $this->m; ++$i) {
106
                if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
107
                    $p = $i;
108
                }
109
            }
110
            if ($p != $j) {
111 View Code Duplication
                for ($k = 0; $k < $this->n; ++$k) {
0 ignored issues
show
Duplication introduced 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...
112
                    $t = $this->LU[$p][$k];
113
                    $this->LU[$p][$k] = $this->LU[$j][$k];
114
                    $this->LU[$j][$k] = $t;
115
                }
116
                $k = $this->piv[$p];
117
                $this->piv[$p] = $this->piv[$j];
118
                $this->piv[$j] = $k;
119
                $this->pivsign = $this->pivsign * -1;
120
            }
121
            // Compute multipliers.
122
            if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
123 View Code Duplication
                for ($i = $j+1; $i < $this->m; ++$i) {
0 ignored issues
show
Duplication introduced 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...
124
                    $this->LU[$i][$j] /= $this->LU[$j][$j];
125
                }
126
            }
127
        }
128
    }    //	function __construct()
0 ignored issues
show
Unused Code Comprehensibility introduced by
50% of this comment could be valid code. Did you maybe forget this after debugging?

Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it.

The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production.

This check looks for comments that seem to be mostly valid code and reports them.

Loading history...
129
130
131
    /**
132
     *	Get lower triangular factor.
133
     *
134
     *	@return array Lower triangular factor
135
     */
136
    public function getL()
137
    {
138
        $L = [];
139
        for ($i = 0; $i < $this->m; ++$i) {
140
            for ($j = 0; $j < $this->n; ++$j) {
141
                if ($i > $j) {
142
                    $L[$i][$j] = $this->LU[$i][$j];
143
                } elseif ($i == $j) {
144
                    $L[$i][$j] = 1.0;
145
                } else {
146
                    $L[$i][$j] = 0.0;
147
                }
148
            }
149
        }
150
        return new Matrix($L);
151
    }    //	function getL()
0 ignored issues
show
Unused Code Comprehensibility introduced by
50% of this comment could be valid code. Did you maybe forget this after debugging?

Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it.

The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production.

This check looks for comments that seem to be mostly valid code and reports them.

Loading history...
152
153
154
    /**
155
     *	Get upper triangular factor.
156
     *
157
     *	@return array Upper triangular factor
158
     */
159
    public function getU()
160
    {
161
        $U = [];
162
        for ($i = 0; $i < $this->n; ++$i) {
163
            for ($j = 0; $j < $this->n; ++$j) {
164
                if ($i <= $j) {
165
                    $U[$i][$j] = $this->LU[$i][$j];
166
                } else {
167
                    $U[$i][$j] = 0.0;
168
                }
169
            }
170
        }
171
        return new Matrix($U);
172
    }    //	function getU()
0 ignored issues
show
Unused Code Comprehensibility introduced by
50% of this comment could be valid code. Did you maybe forget this after debugging?

Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it.

The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production.

This check looks for comments that seem to be mostly valid code and reports them.

Loading history...
173
174
175
    /**
176
     *	Return pivot permutation vector.
177
     *
178
     *	@return array Pivot vector
179
     */
180
    public function getPivot()
181
    {
182
        return $this->piv;
183
    }    //	function getPivot()
0 ignored issues
show
Unused Code Comprehensibility introduced by
50% of this comment could be valid code. Did you maybe forget this after debugging?

Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it.

The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production.

This check looks for comments that seem to be mostly valid code and reports them.

Loading history...
184
185
186
    /**
187
     *	Alias for getPivot
188
     *
189
     *	@see getPivot
190
     */
191
    public function getDoublePivot()
192
    {
193
        return $this->getPivot();
194
    }    //	function getDoublePivot()
0 ignored issues
show
Unused Code Comprehensibility introduced by
50% of this comment could be valid code. Did you maybe forget this after debugging?

Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it.

The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production.

This check looks for comments that seem to be mostly valid code and reports them.

Loading history...
195
196
197
    /**
198
     *	Is the matrix nonsingular?
199
     *
200
     *	@return true if U, and hence A, is nonsingular.
201
     */
202
    public function isNonsingular()
203
    {
204 View Code Duplication
        for ($j = 0; $j < $this->n; ++$j) {
0 ignored issues
show
Duplication introduced 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...
205
            if ($this->LU[$j][$j] == 0) {
206
                return false;
207
            }
208
        }
209
        return true;
210
    }    //	function isNonsingular()
0 ignored issues
show
Unused Code Comprehensibility introduced by
50% of this comment could be valid code. Did you maybe forget this after debugging?

Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it.

The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production.

This check looks for comments that seem to be mostly valid code and reports them.

Loading history...
211
212
213
    /**
214
     *	Count determinants
215
     *
216
     *	@return array d matrix deterninat
217
     */
218
    public function det()
219
    {
220
        if ($this->m == $this->n) {
221
            $d = $this->pivsign;
222 View Code Duplication
            for ($j = 0; $j < $this->n; ++$j) {
0 ignored issues
show
Duplication introduced 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...
223
                $d *= $this->LU[$j][$j];
224
            }
225
            return $d;
226
        } else {
227
            throw MatrixException::notSquareMatrix();
228
        }
229
    }    //	function det()
0 ignored issues
show
Unused Code Comprehensibility introduced by
50% of this comment could be valid code. Did you maybe forget this after debugging?

Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it.

The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production.

This check looks for comments that seem to be mostly valid code and reports them.

Loading history...
230
231
232
    /**
233
     *	Solve A*X = B
234
     *
235
     * @param Matrix $B A Matrix with as many rows as A and any number of columns.
236
     *
237
     * @return array X so that L*U*X = B(piv,:)
238
     *
239
     * @throws MatrixException
240
     */
241
    public function solve(Matrix $B)
242
    {
243
        if ($B->getRows() != $this->m) {
244
            throw MatrixException::notSquareMatrix();
245
        }
246
247
        if (! $this->isNonsingular()) {
248
            throw MatrixException::singularMatrix();
249
        }
250
251
        // Copy right hand side with pivoting
252
        $nx = $B->getColumns();
253
        $X  = $this->getSubMatrix($B->toArray(), $this->piv, 0, $nx-1);
254
        // Solve L*Y = B(piv,:)
255
        for ($k = 0; $k < $this->n; ++$k) {
256 View Code Duplication
            for ($i = $k+1; $i < $this->n; ++$i) {
0 ignored issues
show
Duplication introduced 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...
257
                for ($j = 0; $j < $nx; ++$j) {
258
                    $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
259
                }
260
            }
261
        }
262
        // Solve U*X = Y;
263
        for ($k = $this->n-1; $k >= 0; --$k) {
264
            for ($j = 0; $j < $nx; ++$j) {
265
                $X[$k][$j] /= $this->LU[$k][$k];
266
            }
267 View Code Duplication
            for ($i = 0; $i < $k; ++$i) {
0 ignored issues
show
Duplication introduced 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...
268
                for ($j = 0; $j < $nx; ++$j) {
269
                    $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
270
                }
271
            }
272
        }
273
        return $X;
274
    }    //	function solve()
0 ignored issues
show
Unused Code Comprehensibility introduced by
50% of this comment could be valid code. Did you maybe forget this after debugging?

Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it.

The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production.

This check looks for comments that seem to be mostly valid code and reports them.

Loading history...
275
276
    /**
277
     * @param Matrix $matrix
278
     * @param int $j0
279
     * @param int $jF
280
     *
281
     * @return array
282
     */
283
    protected function getSubMatrix(array $matrix, array $RL, int $j0, int $jF)
284
    {
285
        $m = count($RL);
286
        $n = $jF - $j0;
287
        $R = array_fill(0, $m, array_fill(0, $n+1, 0.0));
288
289
        for ($i = 0; $i < $m; ++$i) {
290
            for ($j = $j0; $j <= $jF; ++$j) {
291
                $R[$i][$j - $j0]= $matrix[ $RL[$i] ][$j];
292
            }
293
        }
294
295
        return $R;
296
    }
297
}    //	class LUDecomposition
298