3 files changed, 478 insertions, 0 deletions

diff --git a/vendor/markbaker/matrix/classes/src/Decomposition/Decomposition.php b/vendor/markbaker/matrix/classes/src/Decomposition/Decomposition.php
new file mode 100644
index 0000000..d4e76f9
--- /dev/null
+++ b/vendor/markbaker/matrix/classes/src/Decomposition/Decomposition.php
@@ -0,0 +1,27 @@
+<?php

+

+namespace Matrix\Decomposition;

+

+use Matrix\Exception;

+use Matrix\Matrix;

+

+class Decomposition

+{

+    const LU = 'LU';

+    const QR = 'QR';

+

+    /**

+     * @throws Exception

+     */

+    public static function decomposition($type, Matrix $matrix)

+    {

+        switch (strtoupper($type)) {

+            case self::LU:

+                return new LU($matrix);

+            case self::QR:

+                return new QR($matrix);

+            default:

+                throw new Exception('Invalid Decomposition');

+        }

+    }

+}

diff --git a/vendor/markbaker/matrix/classes/src/Decomposition/LU.php b/vendor/markbaker/matrix/classes/src/Decomposition/LU.php
new file mode 100644
index 0000000..0a92ed0
--- /dev/null
+++ b/vendor/markbaker/matrix/classes/src/Decomposition/LU.php
@@ -0,0 +1,260 @@
+<?php

+

+namespace Matrix\Decomposition;

+

+use Matrix\Exception;

+use Matrix\Matrix;

+

+class LU

+{

+    private $luMatrix;

+    private $rows;

+    private $columns;

+

+    private $pivot = [];

+

+    public function __construct(Matrix $matrix)

+    {

+        $this->luMatrix = $matrix->toArray();

+        $this->rows = $matrix->rows;

+        $this->columns = $matrix->columns;

+

+        $this->buildPivot();

+    }

+

+    /**

+     * Get lower triangular factor.

+     *

+     * @return Matrix Lower triangular factor

+     */

+    public function getL(): Matrix

+    {

+        $lower = [];

+

+        $columns = min($this->rows, $this->columns);

+        for ($row = 0; $row < $this->rows; ++$row) {

+            for ($column = 0; $column < $columns; ++$column) {

+                if ($row > $column) {

+                    $lower[$row][$column] = $this->luMatrix[$row][$column];

+                } elseif ($row === $column) {

+                    $lower[$row][$column] = 1.0;

+                } else {

+                    $lower[$row][$column] = 0.0;

+                }

+            }

+        }

+

+        return new Matrix($lower);

+    }

+

+    /**

+     * Get upper triangular factor.

+     *

+     * @return Matrix Upper triangular factor

+     */

+    public function getU(): Matrix

+    {

+        $upper = [];

+

+        $rows = min($this->rows, $this->columns);

+        for ($row = 0; $row < $rows; ++$row) {

+            for ($column = 0; $column < $this->columns; ++$column) {

+                if ($row <= $column) {

+                    $upper[$row][$column] = $this->luMatrix[$row][$column];

+                } else {

+                    $upper[$row][$column] = 0.0;

+                }

+            }

+        }

+

+        return new Matrix($upper);

+    }

+

+    /**

+     * Return pivot permutation vector.

+     *

+     * @return Matrix Pivot matrix

+     */

+    public function getP(): Matrix

+    {

+        $pMatrix = [];

+

+        $pivots = $this->pivot;

+        $pivotCount = count($pivots);

+        foreach ($pivots as $row => $pivot) {

+            $pMatrix[$row] = array_fill(0, $pivotCount, 0);

+            $pMatrix[$row][$pivot] = 1;

+        }

+

+        return new Matrix($pMatrix);

+    }

+

+    /**

+     * Return pivot permutation vector.

+     *

+     * @return array Pivot vector

+     */

+    public function getPivot(): array

+    {

+        return $this->pivot;

+    }

+

+    /**

+     *    Is the matrix nonsingular?

+     *

+     * @return bool true if U, and hence A, is nonsingular

+     */

+    public function isNonsingular(): bool

+    {

+        for ($diagonal = 0; $diagonal < $this->columns; ++$diagonal) {

+            if ($this->luMatrix[$diagonal][$diagonal] === 0.0) {

+                return false;

+            }

+        }

+

+        return true;

+    }

+

+    private function buildPivot(): void

+    {

+        for ($row = 0; $row < $this->rows; ++$row) {

+            $this->pivot[$row] = $row;

+        }

+

+        for ($column = 0; $column < $this->columns; ++$column) {

+            $luColumn = $this->localisedReferenceColumn($column);

+

+            $this->applyTransformations($column, $luColumn);

+

+            $pivot = $this->findPivot($column, $luColumn);

+            if ($pivot !== $column) {

+                $this->pivotExchange($pivot, $column);

+            }

+

+            $this->computeMultipliers($column);

+

+            unset($luColumn);

+        }

+    }

+

+    private function localisedReferenceColumn($column): array

+    {

+        $luColumn = [];

+

+        for ($row = 0; $row < $this->rows; ++$row) {

+            $luColumn[$row] = &$this->luMatrix[$row][$column];

+        }

+

+        return $luColumn;

+    }

+

+    private function applyTransformations($column, array $luColumn): void

+    {

+        for ($row = 0; $row < $this->rows; ++$row) {

+            $luRow = $this->luMatrix[$row];

+            // Most of the time is spent in the following dot product.

+            $kmax = min($row, $column);

+            $sValue = 0.0;

+            for ($kValue = 0; $kValue < $kmax; ++$kValue) {

+                $sValue += $luRow[$kValue] * $luColumn[$kValue];

+            }

+            $luRow[$column] = $luColumn[$row] -= $sValue;

+        }

+    }

+

+    private function findPivot($column, array $luColumn): int

+    {

+        $pivot = $column;

+        for ($row = $column + 1; $row < $this->rows; ++$row) {

+            if (abs($luColumn[$row]) > abs($luColumn[$pivot])) {

+                $pivot = $row;

+            }

+        }

+

+        return $pivot;

+    }

+

+    private function pivotExchange($pivot, $column): void

+    {

+        for ($kValue = 0; $kValue < $this->columns; ++$kValue) {

+            $tValue = $this->luMatrix[$pivot][$kValue];

+            $this->luMatrix[$pivot][$kValue] = $this->luMatrix[$column][$kValue];

+            $this->luMatrix[$column][$kValue] = $tValue;

+        }

+

+        $lValue = $this->pivot[$pivot];

+        $this->pivot[$pivot] = $this->pivot[$column];

+        $this->pivot[$column] = $lValue;

+    }

+

+    private function computeMultipliers($diagonal): void

+    {

+        if (($diagonal < $this->rows) && ($this->luMatrix[$diagonal][$diagonal] != 0.0)) {

+            for ($row = $diagonal + 1; $row < $this->rows; ++$row) {

+                $this->luMatrix[$row][$diagonal] /= $this->luMatrix[$diagonal][$diagonal];

+            }

+        }

+    }

+

+    private function pivotB(Matrix $B): array

+    {

+        $X = [];

+        foreach ($this->pivot as $rowId) {

+            $row = $B->getRows($rowId + 1)->toArray();

+            $X[] = array_pop($row);

+        }

+

+        return $X;

+    }

+

+    /**

+     * Solve A*X = B.

+     *

+     * @param Matrix $B a Matrix with as many rows as A and any number of columns

+     *

+     * @throws Exception

+     *

+     * @return Matrix X so that L*U*X = B(piv,:)

+     */

+    public function solve(Matrix $B): Matrix

+    {

+        if ($B->rows !== $this->rows) {

+            throw new Exception('Matrix row dimensions are not equal');

+        }

+

+        if ($this->rows !== $this->columns) {

+            throw new Exception('LU solve() only works on square matrices');

+        }

+

+        if (!$this->isNonsingular()) {

+            throw new Exception('Can only perform operation on singular matrix');

+        }

+

+        // Copy right hand side with pivoting

+        $nx = $B->columns;

+        $X = $this->pivotB($B);

+

+        // Solve L*Y = B(piv,:)

+        for ($k = 0; $k < $this->columns; ++$k) {

+            for ($i = $k + 1; $i < $this->columns; ++$i) {

+                for ($j = 0; $j < $nx; ++$j) {

+                    $X[$i][$j] -= $X[$k][$j] * $this->luMatrix[$i][$k];

+                }

+            }

+        }

+

+        // Solve U*X = Y;

+        for ($k = $this->columns - 1; $k >= 0; --$k) {

+            for ($j = 0; $j < $nx; ++$j) {

+                $X[$k][$j] /= $this->luMatrix[$k][$k];

+            }

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

+                for ($j = 0; $j < $nx; ++$j) {

+                    $X[$i][$j] -= $X[$k][$j] * $this->luMatrix[$i][$k];

+                }

+            }

+        }

+

+        return new Matrix($X);

+    }

+}

diff --git a/vendor/markbaker/matrix/classes/src/Decomposition/QR.php b/vendor/markbaker/matrix/classes/src/Decomposition/QR.php
new file mode 100644
index 0000000..93bf414
--- /dev/null
+++ b/vendor/markbaker/matrix/classes/src/Decomposition/QR.php
@@ -0,0 +1,191 @@
+<?php

+

+namespace Matrix\Decomposition;

+

+use Matrix\Exception;

+use Matrix\Matrix;

+

+class QR

+{

+    private $qrMatrix;

+    private $rows;

+    private $columns;

+

+    private $rDiagonal = [];

+

+    public function __construct(Matrix $matrix)

+    {

+        $this->qrMatrix = $matrix->toArray();

+        $this->rows = $matrix->rows;

+        $this->columns = $matrix->columns;

+

+        $this->decompose();

+    }

+

+    public function getHouseholdVectors(): Matrix

+    {

+        $householdVectors = [];

+        for ($row = 0; $row < $this->rows; ++$row) {

+            for ($column = 0; $column < $this->columns; ++$column) {

+                if ($row >= $column) {

+                    $householdVectors[$row][$column] = $this->qrMatrix[$row][$column];

+                } else {

+                    $householdVectors[$row][$column] = 0.0;

+                }

+            }

+        }

+

+        return new Matrix($householdVectors);

+    }

+

+    public function getQ(): Matrix

+    {

+        $qGrid = [];

+

+        $rowCount = $this->rows;

+        for ($k = $this->columns - 1; $k >= 0; --$k) {

+            for ($i = 0; $i < $this->rows; ++$i) {

+                $qGrid[$i][$k] = 0.0;

+            }

+            $qGrid[$k][$k] = 1.0;

+            if ($this->columns > $this->rows) {

+                $qGrid = array_slice($qGrid, 0, $this->rows);

+            }

+

+            for ($j = $k; $j < $this->columns; ++$j) {

+                if (isset($this->qrMatrix[$k], $this->qrMatrix[$k][$k]) && $this->qrMatrix[$k][$k] != 0.0) {

+                    $s = 0.0;

+                    for ($i = $k; $i < $this->rows; ++$i) {

+                        $s += $this->qrMatrix[$i][$k] * $qGrid[$i][$j];

+                    }

+                    $s = -$s / $this->qrMatrix[$k][$k];

+                    for ($i = $k; $i < $this->rows; ++$i) {

+                        $qGrid[$i][$j] += $s * $this->qrMatrix[$i][$k];

+                    }

+                }

+            }

+        }

+

+        array_walk(

+            $qGrid,

+            function (&$row) use ($rowCount) {

+                $row = array_reverse($row);

+                $row = array_slice($row, 0, $rowCount);

+            }

+        );

+

+        return new Matrix($qGrid);

+    }

+

+    public function getR(): Matrix

+    {

+        $rGrid = [];

+

+        for ($row = 0; $row < $this->columns; ++$row) {

+            for ($column = 0; $column < $this->columns; ++$column) {

+                if ($row < $column) {

+                    $rGrid[$row][$column] = $this->qrMatrix[$row][$column] ?? 0.0;

+                } elseif ($row === $column) {

+                    $rGrid[$row][$column] = $this->rDiagonal[$row] ?? 0.0;

+                } else {

+                    $rGrid[$row][$column] = 0.0;

+                }

+            }

+        }

+

+        if ($this->columns > $this->rows) {

+            $rGrid = array_slice($rGrid, 0, $this->rows);

+        }

+

+        return new Matrix($rGrid);

+    }

+

+    private function hypo($a, $b): float

+    {

+        if (abs($a) > abs($b)) {

+            $r = $b / $a;

+            $r = abs($a) * sqrt(1 + $r * $r);

+        } elseif ($b != 0.0) {

+            $r = $a / $b;

+            $r = abs($b) * sqrt(1 + $r * $r);

+        } else {

+            $r = 0.0;

+        }

+

+        return $r;

+    }

+

+    /**

+     * QR Decomposition computed by Householder reflections.

+     */

+    private function decompose(): void

+    {

+        for ($k = 0; $k < $this->columns; ++$k) {

+            // Compute 2-norm of k-th column without under/overflow.

+            $norm = 0.0;

+            for ($i = $k; $i < $this->rows; ++$i) {

+                $norm = $this->hypo($norm, $this->qrMatrix[$i][$k]);

+            }

+            if ($norm != 0.0) {

+                // Form k-th Householder vector.

+                if ($this->qrMatrix[$k][$k] < 0.0) {

+                    $norm = -$norm;

+                }

+                for ($i = $k; $i < $this->rows; ++$i) {

+                    $this->qrMatrix[$i][$k] /= $norm;

+                }

+                $this->qrMatrix[$k][$k] += 1.0;

+                // Apply transformation to remaining columns.

+                for ($j = $k + 1; $j < $this->columns; ++$j) {

+                    $s = 0.0;

+                    for ($i = $k; $i < $this->rows; ++$i) {

+                        $s += $this->qrMatrix[$i][$k] * $this->qrMatrix[$i][$j];

+                    }

+                    $s = -$s / $this->qrMatrix[$k][$k];

+                    for ($i = $k; $i < $this->rows; ++$i) {

+                        $this->qrMatrix[$i][$j] += $s * $this->qrMatrix[$i][$k];

+                    }

+                }

+            }

+            $this->rDiagonal[$k] = -$norm;

+        }

+    }

+

+    public function isFullRank(): bool

+    {

+        for ($j = 0; $j < $this->columns; ++$j) {

+            if ($this->rDiagonal[$j] == 0.0) {

+                return false;

+            }

+        }

+

+        return true;

+    }

+

+    /**

+     * Least squares solution of A*X = B.

+     *

+     * @param Matrix $B a Matrix with as many rows as A and any number of columns

+     *

+     * @throws Exception

+     *

+     * @return Matrix matrix that minimizes the two norm of Q*R*X-B

+     */

+    public function solve(Matrix $B): Matrix

+    {

+        if ($B->rows !== $this->rows) {

+            throw new Exception('Matrix row dimensions are not equal');

+        }

+

+        if (!$this->isFullRank()) {

+            throw new Exception('Can only perform this operation on a full-rank matrix');

+        }

+

+        // Compute Y = transpose(Q)*B

+        $Y = $this->getQ()->transpose()

+            ->multiply($B);

+        // Solve R*X = Y;

+        return $this->getR()->inverse()

+            ->multiply($Y);

+    }

+}