LUDecomposition::getPivot() - Code Metrics - php-ai/php-ml - Measure and Improve Code Quality continuously with Scrutinizer

LUDecomposition::getPivot() A
last analyzed 2020-05-15 05:48 UTC

↳ Parent: LUDecomposition

Complexity

Conditions	1
Paths	1

Size

Total Lines	3
Code Lines	1

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	1
dl	0
loc	3
rs	10
c	0
b	0
f	0
cc	1
nc	1
nop	0

<?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] ?? 0) > abs($LUcolj[$p] ?? 0)) {
                    $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 *= -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(): array
    {
        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			for ($i = 0; $i < $this->m; ++$i) {
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] ?? 0) > abs($LUcolj[$p] ?? 0)) {
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 *= -1;
137			}
138
139			// Compute multipliers.
140			if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
141			for ($i = $j + 1; $i < $this->m; ++$i) {
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(): array
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			for ($j = 0; $j < $this->n; ++$j) {
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			for ($j = 0; $j < $this->n; ++$j) {
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			for ($i = $k + 1; $i < $this->n; ++$i) {
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			for ($i = 0; $i < $k; ++$i) {
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

LUDecomposition::getPivot() A last analyzed 2020-05-15 05:48 UTC

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like

LUDecomposition::getPivot() A
last analyzed 2020-05-15 05:48 UTC