LUDecomposition::solve() - Code Metrics - Inspection of "Add .gitattributes" - php-ai/php-ml - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Pull Request — master (#287)

by Marcin

created 2018-06-24 18:53 UTC

LUDecomposition::solve() B

↳ Parent: LUDecomposition

Complexity

Conditions	10
Paths	30

Size

Total Lines

Duplication

Lines	10
Ratio	27.03 %

Importance

Changes

Metric	Value
dl	10
loc	37
rs	7.6666
c	0
b	0
f	0
cc	10
nc	30
nop	1

How to fix Complexity

<?php

declare(strict_types=1);

/**
 * @package JAMA
 *
 * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
 * unit lower triangular matrix L, an n-by-n upper triangular matrix U,
 * and a permutation vector piv of length m so that A(piv,:) = L*U.
 * If m < n, then L is m-by-m and U is m-by-n.
 *
 * The LU decompostion with pivoting always exists, even if the matrix is
 * singular, so the constructor will never fail. The primary use of the
 * LU decomposition is in the solution of square systems of simultaneous
 * linear equations. This will fail if isNonsingular() returns false.
 *
 * @author Paul Meagher
 * @author Bartosz Matosiuk
 * @author Michael Bommarito
 *
 * @version 1.1
 *
 * @license PHP v3.0
 *
 *  Slightly changed to adapt the original code to PHP-ML library
 *  @date 2017/04/24
 *
 *  @author Mustafa Karabulut
 */

namespace Phpml\Math\LinearAlgebra;

use Phpml\Exception\MatrixException;
use Phpml\Math\Matrix;

class LUDecomposition
{
    /**
     * Decomposition storage
     *
     * @var array
     */
    private $LU = [];

    /**
     * Row dimension.
     *
     * @var int
     */
    private $m;

    /**
     * Column dimension.
     *
     * @var int
     */
    private $n;

    /**
     * Pivot sign.
     *
     * @var int
     */
    private $pivsign;

    /**
     * Internal storage of pivot vector.
     *
     * @var array
     */
    private $piv = [];

    /**
     * Constructs Structure to access L, U and piv.
     *
     * @param Matrix $A Rectangular matrix
     *
     * @throws MatrixException
     */
    public function __construct(Matrix $A)
    {
        if ($A->getRows() != $A->getColumns()) {
            throw new MatrixException('Matrix is not square matrix');
        }

        // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
        $this->LU = $A->toArray();
        $this->m = $A->getRows();
        $this->n = $A->getColumns();
        for ($i = 0; $i < $this->m; ++$i) {
            $this->piv[$i] = $i;
        }

        $this->pivsign = 1;
        $LUcolj = [];

        // Outer loop.
        for ($j = 0; $j < $this->n; ++$j) {
            // Make a copy of the j-th column to localize references.
            for ($i = 0; $i < $this->m; ++$i) {

                $LUcolj[$i] = &$this->LU[$i][$j];
            }

            // Apply previous transformations.
            for ($i = 0; $i < $this->m; ++$i) {
                $LUrowi = $this->LU[$i];
                // Most of the time is spent in the following dot product.
                $kmax = min($i, $j);
                $s = 0.0;
                for ($k = 0; $k < $kmax; ++$k) {
                    $s += $LUrowi[$k] * $LUcolj[$k];
                }

                $LUrowi[$j] = $LUcolj[$i] -= $s;
            }

            // Find pivot and exchange if necessary.
            $p = $j;
            for ($i = $j + 1; $i < $this->m; ++$i) {
                if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
                    $p = $i;
                }
            }

            if ($p != $j) {
                for ($k = 0; $k < $this->n; ++$k) {
                    $t = $this->LU[$p][$k];
                    $this->LU[$p][$k] = $this->LU[$j][$k];
                    $this->LU[$j][$k] = $t;
                }

                $k = $this->piv[$p];
                $this->piv[$p] = $this->piv[$j];
                $this->piv[$j] = $k;
                $this->pivsign = $this->pivsign * -1;
            }

            // Compute multipliers.
            if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
                for ($i = $j + 1; $i < $this->m; ++$i) {

                    $this->LU[$i][$j] /= $this->LU[$j][$j];
                }
            }
        }
    }

    /**
     * Get lower triangular factor.
     *
     * @return Matrix Lower triangular factor
     */
    public function getL(): Matrix
    {
        $L = [];
        for ($i = 0; $i < $this->m; ++$i) {
            for ($j = 0; $j < $this->n; ++$j) {
                if ($i > $j) {
                    $L[$i][$j] = $this->LU[$i][$j];
                } elseif ($i == $j) {
                    $L[$i][$j] = 1.0;
                } else {
                    $L[$i][$j] = 0.0;
                }
            }
        }

        return new Matrix($L);
    }

    /**
     * Get upper triangular factor.
     *
     * @return Matrix Upper triangular factor
     */
    public function getU(): Matrix
    {
        $U = [];
        for ($i = 0; $i < $this->n; ++$i) {
            for ($j = 0; $j < $this->n; ++$j) {
                if ($i <= $j) {
                    $U[$i][$j] = $this->LU[$i][$j];
                } else {
                    $U[$i][$j] = 0.0;
                }
            }
        }

        return new Matrix($U);
    }

    /**
     * Return pivot permutation vector.
     *
     * @return array Pivot vector
     */
    public function getPivot(): array
    {
        return $this->piv;
    }

    /**
     * Alias for getPivot
     *
     * @see getPivot
     */
    public function getDoublePivot()
    {
        return $this->getPivot();
    }

    /**
     * Is the matrix nonsingular?
     *
     * @return bool true if U, and hence A, is nonsingular.
     */
    public function isNonsingular(): bool
    {
        for ($j = 0; $j < $this->n; ++$j) {

            if ($this->LU[$j][$j] == 0) {
                return false;
            }
        }

        return true;
    }

    public function det(): float
    {
        $d = $this->pivsign;
        for ($j = 0; $j < $this->n; ++$j) {

            $d *= $this->LU[$j][$j];
        }

        return (float) $d;
    }

    /**
     * Solve A*X = B
     *
     * @param Matrix $B A Matrix with as many rows as A and any number of columns.
     *
     * @return array X so that L*U*X = B(piv,:)
     *
     * @throws MatrixException
     */
    public function solve(Matrix $B): array
    {
        if ($B->getRows() != $this->m) {
            throw new MatrixException('Matrix is not square matrix');
        }

        if (!$this->isNonsingular()) {
            throw new MatrixException('Matrix is singular');
        }

        // Copy right hand side with pivoting
        $nx = $B->getColumns();
        $X = $this->getSubMatrix($B->toArray(), $this->piv, 0, $nx - 1);
        // Solve L*Y = B(piv,:)
        for ($k = 0; $k < $this->n; ++$k) {
            for ($i = $k + 1; $i < $this->n; ++$i) {

                for ($j = 0; $j < $nx; ++$j) {
                    $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
                }
            }
        }

        // Solve U*X = Y;
        for ($k = $this->n - 1; $k >= 0; --$k) {
            for ($j = 0; $j < $nx; ++$j) {
                $X[$k][$j] /= $this->LU[$k][$k];
            }

            for ($i = 0; $i < $k; ++$i) {

                for ($j = 0; $j < $nx; ++$j) {
                    $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
                }
            }
        }

        return $X;
    }

    protected function getSubMatrix(array $matrix, array $RL, int $j0, int $jF): array
    {
        $m = count($RL);
        $n = $jF - $j0;
        $R = array_fill(0, $m, array_fill(0, $n + 1, 0.0));

        for ($i = 0; $i < $m; ++$i) {
            for ($j = $j0; $j <= $jF; ++$j) {
                $R[$i][$j - $j0] = $matrix[$RL[$i]][$j];
            }
        }

        return $R;
    }
}


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 new MatrixException('Matrix is not square matrix');
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 Report Bug Copy Issue Report 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		for ($k = 0; $k < $this->n; ++$k) {
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 Report Bug Copy Issue Report 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 Report Bug Copy Issue Report 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		public function det(): float
229		{
230		$d = $this->pivsign;
231	View Code Duplication	for ($j = 0; $j < $this->n; ++$j) {
		0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by Report Bug Copy Issue Report 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...
232		$d *= $this->LU[$j][$j];
233		}
234
235		return (float) $d;
236		}
237
238		/**
239		* Solve A*X = B
240		*
241		* @param Matrix $B A Matrix with as many rows as A and any number of columns.
242		*
243		* @return array X so that LUX = B(piv,:)
244		*
245		* @throws MatrixException
246		*/
247		public function solve(Matrix $B): array
248		{
249		if ($B->getRows() != $this->m) {
250		throw new MatrixException('Matrix is not square matrix');
251		}
252
253		if (!$this->isNonsingular()) {
254		throw new MatrixException('Matrix is singular');
255		}
256
257		// Copy right hand side with pivoting
258		$nx = $B->getColumns();
259		$X = $this->getSubMatrix($B->toArray(), $this->piv, 0, $nx - 1);
260		// Solve L*Y = B(piv,:)
261		for ($k = 0; $k < $this->n; ++$k) {
262	View Code Duplication	for ($i = $k + 1; $i < $this->n; ++$i) {
		0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by Report Bug Copy Issue Report 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...
263		for ($j = 0; $j < $nx; ++$j) {
264		$X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
265		}
266		}
267		}
268
269		// Solve U*X = Y;
270		for ($k = $this->n - 1; $k >= 0; --$k) {
271		for ($j = 0; $j < $nx; ++$j) {
272		$X[$k][$j] /= $this->LU[$k][$k];
273		}
274
275	View Code Duplication	for ($i = 0; $i < $k; ++$i) {
		0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by Report Bug Copy Issue Report 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...
276		for ($j = 0; $j < $nx; ++$j) {
277		$X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
278		}
279		}
280		}
281
282		return $X;
283		}
284
285		protected function getSubMatrix(array $matrix, array $RL, int $j0, int $jF): array
286		{
287		$m = count($RL);
288		$n = $jF - $j0;
289		$R = array_fill(0, $m, array_fill(0, $n + 1, 0.0));
290
291		for ($i = 0; $i < $m; ++$i) {
292		for ($j = $j0; $j <= $jF; ++$j) {
293		$R[$i][$j - $j0] = $matrix[$RL[$i]][$j];
294		}
295		}
296
297		return $R;
298		}
299		}
300

php-ai / php-ml

Pull Request — master (#287)

LUDecomposition::solve() B

Complexity

Size

Duplication

Importance

How to fix Complexity

Long Method

Duplication Side-by-Side

Filter issues like