diff options
Diffstat (limited to '')
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);
+ }
+}
|