diff options
Diffstat (limited to 'vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA')
8 files changed, 3318 insertions, 0 deletions
diff --git a/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/CHANGELOG.TXT b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/CHANGELOG.TXT new file mode 100644 index 0000000..2fc9cd4 --- /dev/null +++ b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/CHANGELOG.TXT @@ -0,0 +1,16 @@ +Mar 1, 2005 11:15 AST by PM
+
++ For consistency, renamed Math.php to Maths.java, utils to util,
+ tests to test, docs to doc -
+
++ Removed conditional logic from top of Matrix class.
+
++ Switched to using hypo function in Maths.php for all php-hypot calls.
+ NOTE TO SELF: Need to make sure that all decompositions have been
+ switched over to using the bundled hypo.
+
+Feb 25, 2005 at 10:00 AST by PM
+
++ Recommend using simpler Error.php instead of JAMA_Error.php but
+ can be persuaded otherwise.
+
diff --git a/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/CholeskyDecomposition.php b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/CholeskyDecomposition.php new file mode 100644 index 0000000..4f7c666 --- /dev/null +++ b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/CholeskyDecomposition.php @@ -0,0 +1,147 @@ +<?php
+
+namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
+
+use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
+
+/**
+ * Cholesky decomposition class.
+ *
+ * For a symmetric, positive definite matrix A, the Cholesky decomposition
+ * is an lower triangular matrix L so that A = L*L'.
+ *
+ * If the matrix is not symmetric or positive definite, the constructor
+ * returns a partial decomposition and sets an internal flag that may
+ * be queried by the isSPD() method.
+ *
+ * @author Paul Meagher
+ * @author Michael Bommarito
+ *
+ * @version 1.2
+ */
+class CholeskyDecomposition
+{
+ /**
+ * Decomposition storage.
+ *
+ * @var array
+ */
+ private $L = [];
+
+ /**
+ * Matrix row and column dimension.
+ *
+ * @var int
+ */
+ private $m;
+
+ /**
+ * Symmetric positive definite flag.
+ *
+ * @var bool
+ */
+ private $isspd = true;
+
+ /**
+ * CholeskyDecomposition.
+ *
+ * Class constructor - decomposes symmetric positive definite matrix
+ *
+ * @param Matrix $A Matrix square symmetric positive definite matrix
+ */
+ public function __construct(Matrix $A)
+ {
+ $this->L = $A->getArray();
+ $this->m = $A->getRowDimension();
+
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = $i; $j < $this->m; ++$j) {
+ for ($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
+ $sum -= $this->L[$i][$k] * $this->L[$j][$k];
+ }
+ if ($i == $j) {
+ if ($sum >= 0) {
+ $this->L[$i][$i] = sqrt($sum);
+ } else {
+ $this->isspd = false;
+ }
+ } else {
+ if ($this->L[$i][$i] != 0) {
+ $this->L[$j][$i] = $sum / $this->L[$i][$i];
+ }
+ }
+ }
+
+ for ($k = $i + 1; $k < $this->m; ++$k) {
+ $this->L[$i][$k] = 0.0;
+ }
+ }
+ }
+
+ /**
+ * Is the matrix symmetric and positive definite?
+ *
+ * @return bool
+ */
+ public function isSPD()
+ {
+ return $this->isspd;
+ }
+
+ /**
+ * getL.
+ *
+ * Return triangular factor.
+ *
+ * @return Matrix Lower triangular matrix
+ */
+ public function getL()
+ {
+ return new Matrix($this->L);
+ }
+
+ /**
+ * Solve A*X = B.
+ *
+ * @param $B Row-equal matrix
+ *
+ * @return Matrix L * L' * X = B
+ */
+ public function solve(Matrix $B)
+ {
+ if ($B->getRowDimension() == $this->m) {
+ if ($this->isspd) {
+ $X = $B->getArrayCopy();
+ $nx = $B->getColumnDimension();
+
+ for ($k = 0; $k < $this->m; ++$k) {
+ for ($i = $k + 1; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $nx; ++$j) {
+ $X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
+ }
+ }
+ for ($j = 0; $j < $nx; ++$j) {
+ $X[$k][$j] /= $this->L[$k][$k];
+ }
+ }
+
+ for ($k = $this->m - 1; $k >= 0; --$k) {
+ for ($j = 0; $j < $nx; ++$j) {
+ $X[$k][$j] /= $this->L[$k][$k];
+ }
+ for ($i = 0; $i < $k; ++$i) {
+ for ($j = 0; $j < $nx; ++$j) {
+ $X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
+ }
+ }
+ }
+
+ return new Matrix($X, $this->m, $nx);
+ }
+
+ throw new CalculationException(Matrix::MATRIX_SPD_EXCEPTION);
+ }
+
+ throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
+ }
+}
diff --git a/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/EigenvalueDecomposition.php b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/EigenvalueDecomposition.php new file mode 100644 index 0000000..a8fc240 --- /dev/null +++ b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/EigenvalueDecomposition.php @@ -0,0 +1,863 @@ +<?php
+
+namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
+
+/**
+ * Class to obtain eigenvalues and eigenvectors of a real matrix.
+ *
+ * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D
+ * is diagonal and the eigenvector matrix V is orthogonal (i.e.
+ * A = V.times(D.times(V.transpose())) and V.times(V.transpose())
+ * equals the identity matrix).
+ *
+ * If A is not symmetric, then the eigenvalue matrix D is block diagonal
+ * with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
+ * lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
+ * columns of V represent the eigenvectors in the sense that A*V = V*D,
+ * i.e. A.times(V) equals V.times(D). The matrix V may be badly
+ * conditioned, or even singular, so the validity of the equation
+ * A = V*D*inverse(V) depends upon V.cond().
+ *
+ * @author Paul Meagher
+ *
+ * @version 1.1
+ */
+class EigenvalueDecomposition
+{
+ /**
+ * Row and column dimension (square matrix).
+ *
+ * @var int
+ */
+ private $n;
+
+ /**
+ * Arrays for internal storage of eigenvalues.
+ *
+ * @var array
+ */
+ private $d = [];
+
+ private $e = [];
+
+ /**
+ * Array for internal storage of eigenvectors.
+ *
+ * @var array
+ */
+ private $V = [];
+
+ /**
+ * Array for internal storage of nonsymmetric Hessenberg form.
+ *
+ * @var array
+ */
+ private $H = [];
+
+ /**
+ * Working storage for nonsymmetric algorithm.
+ *
+ * @var array
+ */
+ private $ort;
+
+ /**
+ * Used for complex scalar division.
+ *
+ * @var float
+ */
+ private $cdivr;
+
+ private $cdivi;
+
+ /**
+ * Symmetric Householder reduction to tridiagonal form.
+ */
+ private function tred2(): void
+ {
+ // This is derived from the Algol procedures tred2 by
+ // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+ $this->d = $this->V[$this->n - 1];
+ // Householder reduction to tridiagonal form.
+ for ($i = $this->n - 1; $i > 0; --$i) {
+ $i_ = $i - 1;
+ // Scale to avoid under/overflow.
+ $h = $scale = 0.0;
+ $scale += array_sum(array_map('abs', $this->d));
+ if ($scale == 0.0) {
+ $this->e[$i] = $this->d[$i_];
+ $this->d = array_slice($this->V[$i_], 0, $i_);
+ for ($j = 0; $j < $i; ++$j) {
+ $this->V[$j][$i] = $this->V[$i][$j] = 0.0;
+ }
+ } else {
+ // Generate Householder vector.
+ for ($k = 0; $k < $i; ++$k) {
+ $this->d[$k] /= $scale;
+ $h += $this->d[$k] ** 2;
+ }
+ $f = $this->d[$i_];
+ $g = sqrt($h);
+ if ($f > 0) {
+ $g = -$g;
+ }
+ $this->e[$i] = $scale * $g;
+ $h = $h - $f * $g;
+ $this->d[$i_] = $f - $g;
+ for ($j = 0; $j < $i; ++$j) {
+ $this->e[$j] = 0.0;
+ }
+ // Apply similarity transformation to remaining columns.
+ for ($j = 0; $j < $i; ++$j) {
+ $f = $this->d[$j];
+ $this->V[$j][$i] = $f;
+ $g = $this->e[$j] + $this->V[$j][$j] * $f;
+ for ($k = $j + 1; $k <= $i_; ++$k) {
+ $g += $this->V[$k][$j] * $this->d[$k];
+ $this->e[$k] += $this->V[$k][$j] * $f;
+ }
+ $this->e[$j] = $g;
+ }
+ $f = 0.0;
+ for ($j = 0; $j < $i; ++$j) {
+ $this->e[$j] /= $h;
+ $f += $this->e[$j] * $this->d[$j];
+ }
+ $hh = $f / (2 * $h);
+ for ($j = 0; $j < $i; ++$j) {
+ $this->e[$j] -= $hh * $this->d[$j];
+ }
+ for ($j = 0; $j < $i; ++$j) {
+ $f = $this->d[$j];
+ $g = $this->e[$j];
+ for ($k = $j; $k <= $i_; ++$k) {
+ $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
+ }
+ $this->d[$j] = $this->V[$i - 1][$j];
+ $this->V[$i][$j] = 0.0;
+ }
+ }
+ $this->d[$i] = $h;
+ }
+
+ // Accumulate transformations.
+ for ($i = 0; $i < $this->n - 1; ++$i) {
+ $this->V[$this->n - 1][$i] = $this->V[$i][$i];
+ $this->V[$i][$i] = 1.0;
+ $h = $this->d[$i + 1];
+ if ($h != 0.0) {
+ for ($k = 0; $k <= $i; ++$k) {
+ $this->d[$k] = $this->V[$k][$i + 1] / $h;
+ }
+ for ($j = 0; $j <= $i; ++$j) {
+ $g = 0.0;
+ for ($k = 0; $k <= $i; ++$k) {
+ $g += $this->V[$k][$i + 1] * $this->V[$k][$j];
+ }
+ for ($k = 0; $k <= $i; ++$k) {
+ $this->V[$k][$j] -= $g * $this->d[$k];
+ }
+ }
+ }
+ for ($k = 0; $k <= $i; ++$k) {
+ $this->V[$k][$i + 1] = 0.0;
+ }
+ }
+
+ $this->d = $this->V[$this->n - 1];
+ $this->V[$this->n - 1] = array_fill(0, $j, 0.0);
+ $this->V[$this->n - 1][$this->n - 1] = 1.0;
+ $this->e[0] = 0.0;
+ }
+
+ /**
+ * Symmetric tridiagonal QL algorithm.
+ *
+ * This is derived from the Algol procedures tql2, by
+ * Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ * Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ * Fortran subroutine in EISPACK.
+ */
+ private function tql2(): void
+ {
+ for ($i = 1; $i < $this->n; ++$i) {
+ $this->e[$i - 1] = $this->e[$i];
+ }
+ $this->e[$this->n - 1] = 0.0;
+ $f = 0.0;
+ $tst1 = 0.0;
+ $eps = 2.0 ** (-52.0);
+
+ for ($l = 0; $l < $this->n; ++$l) {
+ // Find small subdiagonal element
+ $tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
+ $m = $l;
+ while ($m < $this->n) {
+ if (abs($this->e[$m]) <= $eps * $tst1) {
+ break;
+ }
+ ++$m;
+ }
+ // If m == l, $this->d[l] is an eigenvalue,
+ // otherwise, iterate.
+ if ($m > $l) {
+ $iter = 0;
+ do {
+ // Could check iteration count here.
+ ++$iter;
+ // Compute implicit shift
+ $g = $this->d[$l];
+ $p = ($this->d[$l + 1] - $g) / (2.0 * $this->e[$l]);
+ $r = hypo($p, 1.0);
+ if ($p < 0) {
+ $r *= -1;
+ }
+ $this->d[$l] = $this->e[$l] / ($p + $r);
+ $this->d[$l + 1] = $this->e[$l] * ($p + $r);
+ $dl1 = $this->d[$l + 1];
+ $h = $g - $this->d[$l];
+ for ($i = $l + 2; $i < $this->n; ++$i) {
+ $this->d[$i] -= $h;
+ }
+ $f += $h;
+ // Implicit QL transformation.
+ $p = $this->d[$m];
+ $c = 1.0;
+ $c2 = $c3 = $c;
+ $el1 = $this->e[$l + 1];
+ $s = $s2 = 0.0;
+ for ($i = $m - 1; $i >= $l; --$i) {
+ $c3 = $c2;
+ $c2 = $c;
+ $s2 = $s;
+ $g = $c * $this->e[$i];
+ $h = $c * $p;
+ $r = hypo($p, $this->e[$i]);
+ $this->e[$i + 1] = $s * $r;
+ $s = $this->e[$i] / $r;
+ $c = $p / $r;
+ $p = $c * $this->d[$i] - $s * $g;
+ $this->d[$i + 1] = $h + $s * ($c * $g + $s * $this->d[$i]);
+ // Accumulate transformation.
+ for ($k = 0; $k < $this->n; ++$k) {
+ $h = $this->V[$k][$i + 1];
+ $this->V[$k][$i + 1] = $s * $this->V[$k][$i] + $c * $h;
+ $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
+ }
+ }
+ $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
+ $this->e[$l] = $s * $p;
+ $this->d[$l] = $c * $p;
+ // Check for convergence.
+ } while (abs($this->e[$l]) > $eps * $tst1);
+ }
+ $this->d[$l] = $this->d[$l] + $f;
+ $this->e[$l] = 0.0;
+ }
+
+ // Sort eigenvalues and corresponding vectors.
+ for ($i = 0; $i < $this->n - 1; ++$i) {
+ $k = $i;
+ $p = $this->d[$i];
+ for ($j = $i + 1; $j < $this->n; ++$j) {
+ if ($this->d[$j] < $p) {
+ $k = $j;
+ $p = $this->d[$j];
+ }
+ }
+ if ($k != $i) {
+ $this->d[$k] = $this->d[$i];
+ $this->d[$i] = $p;
+ for ($j = 0; $j < $this->n; ++$j) {
+ $p = $this->V[$j][$i];
+ $this->V[$j][$i] = $this->V[$j][$k];
+ $this->V[$j][$k] = $p;
+ }
+ }
+ }
+ }
+
+ /**
+ * Nonsymmetric reduction to Hessenberg form.
+ *
+ * This is derived from the Algol procedures orthes and ortran,
+ * by Martin and Wilkinson, Handbook for Auto. Comp.,
+ * Vol.ii-Linear Algebra, and the corresponding
+ * Fortran subroutines in EISPACK.
+ */
+ private function orthes(): void
+ {
+ $low = 0;
+ $high = $this->n - 1;
+
+ for ($m = $low + 1; $m <= $high - 1; ++$m) {
+ // Scale column.
+ $scale = 0.0;
+ for ($i = $m; $i <= $high; ++$i) {
+ $scale = $scale + abs($this->H[$i][$m - 1]);
+ }
+ if ($scale != 0.0) {
+ // Compute Householder transformation.
+ $h = 0.0;
+ for ($i = $high; $i >= $m; --$i) {
+ $this->ort[$i] = $this->H[$i][$m - 1] / $scale;
+ $h += $this->ort[$i] * $this->ort[$i];
+ }
+ $g = sqrt($h);
+ if ($this->ort[$m] > 0) {
+ $g *= -1;
+ }
+ $h -= $this->ort[$m] * $g;
+ $this->ort[$m] -= $g;
+ // Apply Householder similarity transformation
+ // H = (I -u * u' / h) * H * (I -u * u') / h)
+ for ($j = $m; $j < $this->n; ++$j) {
+ $f = 0.0;
+ for ($i = $high; $i >= $m; --$i) {
+ $f += $this->ort[$i] * $this->H[$i][$j];
+ }
+ $f /= $h;
+ for ($i = $m; $i <= $high; ++$i) {
+ $this->H[$i][$j] -= $f * $this->ort[$i];
+ }
+ }
+ for ($i = 0; $i <= $high; ++$i) {
+ $f = 0.0;
+ for ($j = $high; $j >= $m; --$j) {
+ $f += $this->ort[$j] * $this->H[$i][$j];
+ }
+ $f = $f / $h;
+ for ($j = $m; $j <= $high; ++$j) {
+ $this->H[$i][$j] -= $f * $this->ort[$j];
+ }
+ }
+ $this->ort[$m] = $scale * $this->ort[$m];
+ $this->H[$m][$m - 1] = $scale * $g;
+ }
+ }
+
+ // Accumulate transformations (Algol's ortran).
+ for ($i = 0; $i < $this->n; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
+ }
+ }
+ for ($m = $high - 1; $m >= $low + 1; --$m) {
+ if ($this->H[$m][$m - 1] != 0.0) {
+ for ($i = $m + 1; $i <= $high; ++$i) {
+ $this->ort[$i] = $this->H[$i][$m - 1];
+ }
+ for ($j = $m; $j <= $high; ++$j) {
+ $g = 0.0;
+ for ($i = $m; $i <= $high; ++$i) {
+ $g += $this->ort[$i] * $this->V[$i][$j];
+ }
+ // Double division avoids possible underflow
+ $g = ($g / $this->ort[$m]) / $this->H[$m][$m - 1];
+ for ($i = $m; $i <= $high; ++$i) {
+ $this->V[$i][$j] += $g * $this->ort[$i];
+ }
+ }
+ }
+ }
+ }
+
+ /**
+ * Performs complex division.
+ *
+ * @param mixed $xr
+ * @param mixed $xi
+ * @param mixed $yr
+ * @param mixed $yi
+ */
+ private function cdiv($xr, $xi, $yr, $yi): void
+ {
+ if (abs($yr) > abs($yi)) {
+ $r = $yi / $yr;
+ $d = $yr + $r * $yi;
+ $this->cdivr = ($xr + $r * $xi) / $d;
+ $this->cdivi = ($xi - $r * $xr) / $d;
+ } else {
+ $r = $yr / $yi;
+ $d = $yi + $r * $yr;
+ $this->cdivr = ($r * $xr + $xi) / $d;
+ $this->cdivi = ($r * $xi - $xr) / $d;
+ }
+ }
+
+ /**
+ * Nonsymmetric reduction from Hessenberg to real Schur form.
+ *
+ * Code is derived from the Algol procedure hqr2,
+ * by Martin and Wilkinson, Handbook for Auto. Comp.,
+ * Vol.ii-Linear Algebra, and the corresponding
+ * Fortran subroutine in EISPACK.
+ */
+ private function hqr2(): void
+ {
+ // Initialize
+ $nn = $this->n;
+ $n = $nn - 1;
+ $low = 0;
+ $high = $nn - 1;
+ $eps = 2.0 ** (-52.0);
+ $exshift = 0.0;
+ $p = $q = $r = $s = $z = 0;
+ // Store roots isolated by balanc and compute matrix norm
+ $norm = 0.0;
+
+ for ($i = 0; $i < $nn; ++$i) {
+ if (($i < $low) || ($i > $high)) {
+ $this->d[$i] = $this->H[$i][$i];
+ $this->e[$i] = 0.0;
+ }
+ for ($j = max($i - 1, 0); $j < $nn; ++$j) {
+ $norm = $norm + abs($this->H[$i][$j]);
+ }
+ }
+
+ // Outer loop over eigenvalue index
+ $iter = 0;
+ while ($n >= $low) {
+ // Look for single small sub-diagonal element
+ $l = $n;
+ while ($l > $low) {
+ $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]);
+ if ($s == 0.0) {
+ $s = $norm;
+ }
+ if (abs($this->H[$l][$l - 1]) < $eps * $s) {
+ break;
+ }
+ --$l;
+ }
+ // Check for convergence
+ // One root found
+ if ($l == $n) {
+ $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
+ $this->d[$n] = $this->H[$n][$n];
+ $this->e[$n] = 0.0;
+ --$n;
+ $iter = 0;
+ // Two roots found
+ } elseif ($l == $n - 1) {
+ $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
+ $p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0;
+ $q = $p * $p + $w;
+ $z = sqrt(abs($q));
+ $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
+ $this->H[$n - 1][$n - 1] = $this->H[$n - 1][$n - 1] + $exshift;
+ $x = $this->H[$n][$n];
+ // Real pair
+ if ($q >= 0) {
+ if ($p >= 0) {
+ $z = $p + $z;
+ } else {
+ $z = $p - $z;
+ }
+ $this->d[$n - 1] = $x + $z;
+ $this->d[$n] = $this->d[$n - 1];
+ if ($z != 0.0) {
+ $this->d[$n] = $x - $w / $z;
+ }
+ $this->e[$n - 1] = 0.0;
+ $this->e[$n] = 0.0;
+ $x = $this->H[$n][$n - 1];
+ $s = abs($x) + abs($z);
+ $p = $x / $s;
+ $q = $z / $s;
+ $r = sqrt($p * $p + $q * $q);
+ $p = $p / $r;
+ $q = $q / $r;
+ // Row modification
+ for ($j = $n - 1; $j < $nn; ++$j) {
+ $z = $this->H[$n - 1][$j];
+ $this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j];
+ $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
+ }
+ // Column modification
+ for ($i = 0; $i <= $n; ++$i) {
+ $z = $this->H[$i][$n - 1];
+ $this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n];
+ $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
+ }
+ // Accumulate transformations
+ for ($i = $low; $i <= $high; ++$i) {
+ $z = $this->V[$i][$n - 1];
+ $this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n];
+ $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
+ }
+ // Complex pair
+ } else {
+ $this->d[$n - 1] = $x + $p;
+ $this->d[$n] = $x + $p;
+ $this->e[$n - 1] = $z;
+ $this->e[$n] = -$z;
+ }
+ $n = $n - 2;
+ $iter = 0;
+ // No convergence yet
+ } else {
+ // Form shift
+ $x = $this->H[$n][$n];
+ $y = 0.0;
+ $w = 0.0;
+ if ($l < $n) {
+ $y = $this->H[$n - 1][$n - 1];
+ $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
+ }
+ // Wilkinson's original ad hoc shift
+ if ($iter == 10) {
+ $exshift += $x;
+ for ($i = $low; $i <= $n; ++$i) {
+ $this->H[$i][$i] -= $x;
+ }
+ $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]);
+ $x = $y = 0.75 * $s;
+ $w = -0.4375 * $s * $s;
+ }
+ // MATLAB's new ad hoc shift
+ if ($iter == 30) {
+ $s = ($y - $x) / 2.0;
+ $s = $s * $s + $w;
+ if ($s > 0) {
+ $s = sqrt($s);
+ if ($y < $x) {
+ $s = -$s;
+ }
+ $s = $x - $w / (($y - $x) / 2.0 + $s);
+ for ($i = $low; $i <= $n; ++$i) {
+ $this->H[$i][$i] -= $s;
+ }
+ $exshift += $s;
+ $x = $y = $w = 0.964;
+ }
+ }
+ // Could check iteration count here.
+ $iter = $iter + 1;
+ // Look for two consecutive small sub-diagonal elements
+ $m = $n - 2;
+ while ($m >= $l) {
+ $z = $this->H[$m][$m];
+ $r = $x - $z;
+ $s = $y - $z;
+ $p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1];
+ $q = $this->H[$m + 1][$m + 1] - $z - $r - $s;
+ $r = $this->H[$m + 2][$m + 1];
+ $s = abs($p) + abs($q) + abs($r);
+ $p = $p / $s;
+ $q = $q / $s;
+ $r = $r / $s;
+ if ($m == $l) {
+ break;
+ }
+ if (
+ abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) <
+ $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1])))
+ ) {
+ break;
+ }
+ --$m;
+ }
+ for ($i = $m + 2; $i <= $n; ++$i) {
+ $this->H[$i][$i - 2] = 0.0;
+ if ($i > $m + 2) {
+ $this->H[$i][$i - 3] = 0.0;
+ }
+ }
+ // Double QR step involving rows l:n and columns m:n
+ for ($k = $m; $k <= $n - 1; ++$k) {
+ $notlast = ($k != $n - 1);
+ if ($k != $m) {
+ $p = $this->H[$k][$k - 1];
+ $q = $this->H[$k + 1][$k - 1];
+ $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0);
+ $x = abs($p) + abs($q) + abs($r);
+ if ($x != 0.0) {
+ $p = $p / $x;
+ $q = $q / $x;
+ $r = $r / $x;
+ }
+ }
+ if ($x == 0.0) {
+ break;
+ }
+ $s = sqrt($p * $p + $q * $q + $r * $r);
+ if ($p < 0) {
+ $s = -$s;
+ }
+ if ($s != 0) {
+ if ($k != $m) {
+ $this->H[$k][$k - 1] = -$s * $x;
+ } elseif ($l != $m) {
+ $this->H[$k][$k - 1] = -$this->H[$k][$k - 1];
+ }
+ $p = $p + $s;
+ $x = $p / $s;
+ $y = $q / $s;
+ $z = $r / $s;
+ $q = $q / $p;
+ $r = $r / $p;
+ // Row modification
+ for ($j = $k; $j < $nn; ++$j) {
+ $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j];
+ if ($notlast) {
+ $p = $p + $r * $this->H[$k + 2][$j];
+ $this->H[$k + 2][$j] = $this->H[$k + 2][$j] - $p * $z;
+ }
+ $this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
+ $this->H[$k + 1][$j] = $this->H[$k + 1][$j] - $p * $y;
+ }
+ // Column modification
+ $iMax = min($n, $k + 3);
+ for ($i = 0; $i <= $iMax; ++$i) {
+ $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k + 1];
+ if ($notlast) {
+ $p = $p + $z * $this->H[$i][$k + 2];
+ $this->H[$i][$k + 2] = $this->H[$i][$k + 2] - $p * $r;
+ }
+ $this->H[$i][$k] = $this->H[$i][$k] - $p;
+ $this->H[$i][$k + 1] = $this->H[$i][$k + 1] - $p * $q;
+ }
+ // Accumulate transformations
+ for ($i = $low; $i <= $high; ++$i) {
+ $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1];
+ if ($notlast) {
+ $p = $p + $z * $this->V[$i][$k + 2];
+ $this->V[$i][$k + 2] = $this->V[$i][$k + 2] - $p * $r;
+ }
+ $this->V[$i][$k] = $this->V[$i][$k] - $p;
+ $this->V[$i][$k + 1] = $this->V[$i][$k + 1] - $p * $q;
+ }
+ } // ($s != 0)
+ } // k loop
+ } // check convergence
+ } // while ($n >= $low)
+
+ // Backsubstitute to find vectors of upper triangular form
+ if ($norm == 0.0) {
+ return;
+ }
+
+ for ($n = $nn - 1; $n >= 0; --$n) {
+ $p = $this->d[$n];
+ $q = $this->e[$n];
+ // Real vector
+ if ($q == 0) {
+ $l = $n;
+ $this->H[$n][$n] = 1.0;
+ for ($i = $n - 1; $i >= 0; --$i) {
+ $w = $this->H[$i][$i] - $p;
+ $r = 0.0;
+ for ($j = $l; $j <= $n; ++$j) {
+ $r = $r + $this->H[$i][$j] * $this->H[$j][$n];
+ }
+ if ($this->e[$i] < 0.0) {
+ $z = $w;
+ $s = $r;
+ } else {
+ $l = $i;
+ if ($this->e[$i] == 0.0) {
+ if ($w != 0.0) {
+ $this->H[$i][$n] = -$r / $w;
+ } else {
+ $this->H[$i][$n] = -$r / ($eps * $norm);
+ }
+ // Solve real equations
+ } else {
+ $x = $this->H[$i][$i + 1];
+ $y = $this->H[$i + 1][$i];
+ $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
+ $t = ($x * $s - $z * $r) / $q;
+ $this->H[$i][$n] = $t;
+ if (abs($x) > abs($z)) {
+ $this->H[$i + 1][$n] = (-$r - $w * $t) / $x;
+ } else {
+ $this->H[$i + 1][$n] = (-$s - $y * $t) / $z;
+ }
+ }
+ // Overflow control
+ $t = abs($this->H[$i][$n]);
+ if (($eps * $t) * $t > 1) {
+ for ($j = $i; $j <= $n; ++$j) {
+ $this->H[$j][$n] = $this->H[$j][$n] / $t;
+ }
+ }
+ }
+ }
+ // Complex vector
+ } elseif ($q < 0) {
+ $l = $n - 1;
+ // Last vector component imaginary so matrix is triangular
+ if (abs($this->H[$n][$n - 1]) > abs($this->H[$n - 1][$n])) {
+ $this->H[$n - 1][$n - 1] = $q / $this->H[$n][$n - 1];
+ $this->H[$n - 1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n - 1];
+ } else {
+ $this->cdiv(0.0, -$this->H[$n - 1][$n], $this->H[$n - 1][$n - 1] - $p, $q);
+ $this->H[$n - 1][$n - 1] = $this->cdivr;
+ $this->H[$n - 1][$n] = $this->cdivi;
+ }
+ $this->H[$n][$n - 1] = 0.0;
+ $this->H[$n][$n] = 1.0;
+ for ($i = $n - 2; $i >= 0; --$i) {
+ // double ra,sa,vr,vi;
+ $ra = 0.0;
+ $sa = 0.0;
+ for ($j = $l; $j <= $n; ++$j) {
+ $ra = $ra + $this->H[$i][$j] * $this->H[$j][$n - 1];
+ $sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
+ }
+ $w = $this->H[$i][$i] - $p;
+ if ($this->e[$i] < 0.0) {
+ $z = $w;
+ $r = $ra;
+ $s = $sa;
+ } else {
+ $l = $i;
+ if ($this->e[$i] == 0) {
+ $this->cdiv(-$ra, -$sa, $w, $q);
+ $this->H[$i][$n - 1] = $this->cdivr;
+ $this->H[$i][$n] = $this->cdivi;
+ } else {
+ // Solve complex equations
+ $x = $this->H[$i][$i + 1];
+ $y = $this->H[$i + 1][$i];
+ $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
+ $vi = ($this->d[$i] - $p) * 2.0 * $q;
+ if ($vr == 0.0 & $vi == 0.0) {
+ $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
+ }
+ $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
+ $this->H[$i][$n - 1] = $this->cdivr;
+ $this->H[$i][$n] = $this->cdivi;
+ if (abs($x) > (abs($z) + abs($q))) {
+ $this->H[$i + 1][$n - 1] = (-$ra - $w * $this->H[$i][$n - 1] + $q * $this->H[$i][$n]) / $x;
+ $this->H[$i + 1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n - 1]) / $x;
+ } else {
+ $this->cdiv(-$r - $y * $this->H[$i][$n - 1], -$s - $y * $this->H[$i][$n], $z, $q);
+ $this->H[$i + 1][$n - 1] = $this->cdivr;
+ $this->H[$i + 1][$n] = $this->cdivi;
+ }
+ }
+ // Overflow control
+ $t = max(abs($this->H[$i][$n - 1]), abs($this->H[$i][$n]));
+ if (($eps * $t) * $t > 1) {
+ for ($j = $i; $j <= $n; ++$j) {
+ $this->H[$j][$n - 1] = $this->H[$j][$n - 1] / $t;
+ $this->H[$j][$n] = $this->H[$j][$n] / $t;
+ }
+ }
+ } // end else
+ } // end for
+ } // end else for complex case
+ } // end for
+
+ // Vectors of isolated roots
+ for ($i = 0; $i < $nn; ++$i) {
+ if ($i < $low | $i > $high) {
+ for ($j = $i; $j < $nn; ++$j) {
+ $this->V[$i][$j] = $this->H[$i][$j];
+ }
+ }
+ }
+
+ // Back transformation to get eigenvectors of original matrix
+ for ($j = $nn - 1; $j >= $low; --$j) {
+ for ($i = $low; $i <= $high; ++$i) {
+ $z = 0.0;
+ $kMax = min($j, $high);
+ for ($k = $low; $k <= $kMax; ++$k) {
+ $z = $z + $this->V[$i][$k] * $this->H[$k][$j];
+ }
+ $this->V[$i][$j] = $z;
+ }
+ }
+ }
+
+ // end hqr2
+
+ /**
+ * Constructor: Check for symmetry, then construct the eigenvalue decomposition.
+ *
+ * @param mixed $Arg A Square matrix
+ */
+ public function __construct($Arg)
+ {
+ $this->A = $Arg->getArray();
+ $this->n = $Arg->getColumnDimension();
+
+ $issymmetric = true;
+ for ($j = 0; ($j < $this->n) & $issymmetric; ++$j) {
+ for ($i = 0; ($i < $this->n) & $issymmetric; ++$i) {
+ $issymmetric = ($this->A[$i][$j] == $this->A[$j][$i]);
+ }
+ }
+
+ if ($issymmetric) {
+ $this->V = $this->A;
+ // Tridiagonalize.
+ $this->tred2();
+ // Diagonalize.
+ $this->tql2();
+ } else {
+ $this->H = $this->A;
+ $this->ort = [];
+ // Reduce to Hessenberg form.
+ $this->orthes();
+ // Reduce Hessenberg to real Schur form.
+ $this->hqr2();
+ }
+ }
+
+ /**
+ * Return the eigenvector matrix.
+ *
+ * @return Matrix V
+ */
+ public function getV()
+ {
+ return new Matrix($this->V, $this->n, $this->n);
+ }
+
+ /**
+ * Return the real parts of the eigenvalues.
+ *
+ * @return array real(diag(D))
+ */
+ public function getRealEigenvalues()
+ {
+ return $this->d;
+ }
+
+ /**
+ * Return the imaginary parts of the eigenvalues.
+ *
+ * @return array imag(diag(D))
+ */
+ public function getImagEigenvalues()
+ {
+ return $this->e;
+ }
+
+ /**
+ * Return the block diagonal eigenvalue matrix.
+ *
+ * @return Matrix D
+ */
+ public function getD()
+ {
+ for ($i = 0; $i < $this->n; ++$i) {
+ $D[$i] = array_fill(0, $this->n, 0.0);
+ $D[$i][$i] = $this->d[$i];
+ if ($this->e[$i] == 0) {
+ continue;
+ }
+ $o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
+ $D[$i][$o] = $this->e[$i];
+ }
+
+ return new Matrix($D);
+ }
+}
diff --git a/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/LUDecomposition.php b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/LUDecomposition.php new file mode 100644 index 0000000..3871895 --- /dev/null +++ b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/LUDecomposition.php @@ -0,0 +1,282 @@ +<?php
+
+namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
+
+use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
+
+/**
+ * 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
+ */
+class LUDecomposition
+{
+ const MATRIX_SINGULAR_EXCEPTION = 'Can only perform operation on singular matrix.';
+ const MATRIX_SQUARE_EXCEPTION = 'Mismatched Row dimension';
+
+ /**
+ * 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 = [];
+
+ /**
+ * LU Decomposition constructor.
+ *
+ * @param Matrix $A Rectangular matrix
+ */
+ public function __construct($A)
+ {
+ if ($A instanceof Matrix) {
+ // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
+ $this->LU = $A->getArray();
+ $this->m = $A->getRowDimension();
+ $this->n = $A->getColumnDimension();
+ for ($i = 0; $i < $this->m; ++$i) {
+ $this->piv[$i] = $i;
+ }
+ $this->pivsign = 1;
+ $LUrowi = $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];
+ }
+ }
+ }
+ } else {
+ throw new CalculationException(Matrix::ARGUMENT_TYPE_EXCEPTION);
+ }
+ }
+
+ // function __construct()
+
+ /**
+ * Get lower triangular factor.
+ *
+ * @return Matrix Lower triangular factor
+ */
+ public function getL()
+ {
+ 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);
+ }
+
+ // function getL()
+
+ /**
+ * Get upper triangular factor.
+ *
+ * @return Matrix Upper triangular factor
+ */
+ public function getU()
+ {
+ 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);
+ }
+
+ // function getU()
+
+ /**
+ * Return pivot permutation vector.
+ *
+ * @return array Pivot vector
+ */
+ public function getPivot()
+ {
+ return $this->piv;
+ }
+
+ // function getPivot()
+
+ /**
+ * Alias for getPivot.
+ *
+ * @see getPivot
+ */
+ public function getDoublePivot()
+ {
+ return $this->getPivot();
+ }
+
+ // function getDoublePivot()
+
+ /**
+ * Is the matrix nonsingular?
+ *
+ * @return bool true if U, and hence A, is nonsingular
+ */
+ public function isNonsingular()
+ {
+ for ($j = 0; $j < $this->n; ++$j) {
+ if ($this->LU[$j][$j] == 0) {
+ return false;
+ }
+ }
+
+ return true;
+ }
+
+ // function isNonsingular()
+
+ /**
+ * Count determinants.
+ *
+ * @return array d matrix deterninat
+ */
+ public function det()
+ {
+ if ($this->m == $this->n) {
+ $d = $this->pivsign;
+ for ($j = 0; $j < $this->n; ++$j) {
+ $d *= $this->LU[$j][$j];
+ }
+
+ return $d;
+ }
+
+ throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
+ }
+
+ // function det()
+
+ /**
+ * Solve A*X = B.
+ *
+ * @param mixed $B a Matrix with as many rows as A and any number of columns
+ *
+ * @return Matrix X so that L*U*X = B(piv,:)
+ */
+ public function solve($B)
+ {
+ if ($B->getRowDimension() == $this->m) {
+ if ($this->isNonsingular()) {
+ // Copy right hand side with pivoting
+ $nx = $B->getColumnDimension();
+ $X = $B->getMatrix($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->A[$i][$j] -= $X->A[$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->A[$k][$j] /= $this->LU[$k][$k];
+ }
+ for ($i = 0; $i < $k; ++$i) {
+ for ($j = 0; $j < $nx; ++$j) {
+ $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
+ }
+ }
+ }
+
+ return $X;
+ }
+
+ throw new CalculationException(self::MATRIX_SINGULAR_EXCEPTION);
+ }
+
+ throw new CalculationException(self::MATRIX_SQUARE_EXCEPTION);
+ }
+}
diff --git a/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/Matrix.php b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/Matrix.php new file mode 100644 index 0000000..fc815f6 --- /dev/null +++ b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/Matrix.php @@ -0,0 +1,1202 @@ +<?php
+
+namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
+
+use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
+use PhpOffice\PhpSpreadsheet\Calculation\Functions;
+use PhpOffice\PhpSpreadsheet\Shared\StringHelper;
+
+/**
+ * Matrix class.
+ *
+ * @author Paul Meagher
+ * @author Michael Bommarito
+ * @author Lukasz Karapuda
+ * @author Bartek Matosiuk
+ *
+ * @version 1.8
+ *
+ * @see https://math.nist.gov/javanumerics/jama/
+ */
+class Matrix
+{
+ const POLYMORPHIC_ARGUMENT_EXCEPTION = 'Invalid argument pattern for polymorphic function.';
+ const ARGUMENT_TYPE_EXCEPTION = 'Invalid argument type.';
+ const ARGUMENT_BOUNDS_EXCEPTION = 'Invalid argument range.';
+ const MATRIX_DIMENSION_EXCEPTION = 'Matrix dimensions are not equal.';
+ const ARRAY_LENGTH_EXCEPTION = 'Array length must be a multiple of m.';
+ const MATRIX_SPD_EXCEPTION = 'Can only perform operation on symmetric positive definite matrix.';
+
+ /**
+ * Matrix storage.
+ *
+ * @var array
+ */
+ public $A = [];
+
+ /**
+ * Matrix row dimension.
+ *
+ * @var int
+ */
+ private $m;
+
+ /**
+ * Matrix column dimension.
+ *
+ * @var int
+ */
+ private $n;
+
+ /**
+ * Polymorphic constructor.
+ *
+ * As PHP has no support for polymorphic constructors, we use tricks to make our own sort of polymorphism using func_num_args, func_get_arg, and gettype. In essence, we're just implementing a simple RTTI filter and calling the appropriate constructor.
+ */
+ public function __construct(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ //Rectangular matrix - m x n initialized from 2D array
+ case 'array':
+ $this->m = count($args[0]);
+ $this->n = count($args[0][0]);
+ $this->A = $args[0];
+
+ break;
+ //Square matrix - n x n
+ case 'integer':
+ $this->m = $args[0];
+ $this->n = $args[0];
+ $this->A = array_fill(0, $this->m, array_fill(0, $this->n, 0));
+
+ break;
+ //Rectangular matrix - m x n
+ case 'integer,integer':
+ $this->m = $args[0];
+ $this->n = $args[1];
+ $this->A = array_fill(0, $this->m, array_fill(0, $this->n, 0));
+
+ break;
+ //Rectangular matrix - m x n initialized from packed array
+ case 'array,integer':
+ $this->m = $args[1];
+ if ($this->m != 0) {
+ $this->n = count($args[0]) / $this->m;
+ } else {
+ $this->n = 0;
+ }
+ if (($this->m * $this->n) == count($args[0])) {
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $this->A[$i][$j] = $args[0][$i + $j * $this->m];
+ }
+ }
+ } else {
+ throw new CalculationException(self::ARRAY_LENGTH_EXCEPTION);
+ }
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ } else {
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+ }
+
+ /**
+ * getArray.
+ *
+ * @return array Matrix array
+ */
+ public function getArray()
+ {
+ return $this->A;
+ }
+
+ /**
+ * getRowDimension.
+ *
+ * @return int Row dimension
+ */
+ public function getRowDimension()
+ {
+ return $this->m;
+ }
+
+ /**
+ * getColumnDimension.
+ *
+ * @return int Column dimension
+ */
+ public function getColumnDimension()
+ {
+ return $this->n;
+ }
+
+ /**
+ * get.
+ *
+ * Get the i,j-th element of the matrix.
+ *
+ * @param int $i Row position
+ * @param int $j Column position
+ *
+ * @return mixed Element (int/float/double)
+ */
+ public function get($i = null, $j = null)
+ {
+ return $this->A[$i][$j];
+ }
+
+ /**
+ * getMatrix.
+ *
+ * Get a submatrix
+ *
+ * @return Matrix Submatrix
+ */
+ public function getMatrix(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ //A($i0...; $j0...)
+ case 'integer,integer':
+ [$i0, $j0] = $args;
+ if ($i0 >= 0) {
+ $m = $this->m - $i0;
+ } else {
+ throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
+ }
+ if ($j0 >= 0) {
+ $n = $this->n - $j0;
+ } else {
+ throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
+ }
+ $R = new self($m, $n);
+ for ($i = $i0; $i < $this->m; ++$i) {
+ for ($j = $j0; $j < $this->n; ++$j) {
+ $R->set($i, $j, $this->A[$i][$j]);
+ }
+ }
+
+ return $R;
+
+ break;
+ //A($i0...$iF; $j0...$jF)
+ case 'integer,integer,integer,integer':
+ [$i0, $iF, $j0, $jF] = $args;
+ if (($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0)) {
+ $m = $iF - $i0;
+ } else {
+ throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
+ }
+ if (($jF > $j0) && ($this->n >= $jF) && ($j0 >= 0)) {
+ $n = $jF - $j0;
+ } else {
+ throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
+ }
+ $R = new self($m + 1, $n + 1);
+ for ($i = $i0; $i <= $iF; ++$i) {
+ for ($j = $j0; $j <= $jF; ++$j) {
+ $R->set($i - $i0, $j - $j0, $this->A[$i][$j]);
+ }
+ }
+
+ return $R;
+
+ break;
+ //$R = array of row indices; $C = array of column indices
+ case 'array,array':
+ [$RL, $CL] = $args;
+ if (count($RL) > 0) {
+ $m = count($RL);
+ } else {
+ throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
+ }
+ if (count($CL) > 0) {
+ $n = count($CL);
+ } else {
+ throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
+ }
+ $R = new self($m, $n);
+ for ($i = 0; $i < $m; ++$i) {
+ for ($j = 0; $j < $n; ++$j) {
+ $R->set($i, $j, $this->A[$RL[$i]][$CL[$j]]);
+ }
+ }
+
+ return $R;
+
+ break;
+ //A($i0...$iF); $CL = array of column indices
+ case 'integer,integer,array':
+ [$i0, $iF, $CL] = $args;
+ if (($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0)) {
+ $m = $iF - $i0;
+ } else {
+ throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
+ }
+ if (count($CL) > 0) {
+ $n = count($CL);
+ } else {
+ throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
+ }
+ $R = new self($m, $n);
+ for ($i = $i0; $i < $iF; ++$i) {
+ for ($j = 0; $j < $n; ++$j) {
+ $R->set($i - $i0, $j, $this->A[$i][$CL[$j]]);
+ }
+ }
+
+ return $R;
+
+ break;
+ //$RL = array of row indices
+ case 'array,integer,integer':
+ [$RL, $j0, $jF] = $args;
+ if (count($RL) > 0) {
+ $m = count($RL);
+ } else {
+ throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
+ }
+ if (($jF >= $j0) && ($this->n >= $jF) && ($j0 >= 0)) {
+ $n = $jF - $j0;
+ } else {
+ throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
+ }
+ $R = new self($m, $n + 1);
+ for ($i = 0; $i < $m; ++$i) {
+ for ($j = $j0; $j <= $jF; ++$j) {
+ $R->set($i, $j - $j0, $this->A[$RL[$i]][$j]);
+ }
+ }
+
+ return $R;
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ } else {
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+ }
+
+ /**
+ * checkMatrixDimensions.
+ *
+ * Is matrix B the same size?
+ *
+ * @param Matrix $B Matrix B
+ *
+ * @return bool
+ */
+ public function checkMatrixDimensions($B = null)
+ {
+ if ($B instanceof self) {
+ if (($this->m == $B->getRowDimension()) && ($this->n == $B->getColumnDimension())) {
+ return true;
+ }
+
+ throw new CalculationException(self::MATRIX_DIMENSION_EXCEPTION);
+ }
+
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ // function checkMatrixDimensions()
+
+ /**
+ * set.
+ *
+ * Set the i,j-th element of the matrix.
+ *
+ * @param int $i Row position
+ * @param int $j Column position
+ * @param mixed $c Int/float/double value
+ *
+ * @return mixed Element (int/float/double)
+ */
+ public function set($i = null, $j = null, $c = null)
+ {
+ // Optimized set version just has this
+ $this->A[$i][$j] = $c;
+ }
+
+ // function set()
+
+ /**
+ * identity.
+ *
+ * Generate an identity matrix.
+ *
+ * @param int $m Row dimension
+ * @param int $n Column dimension
+ *
+ * @return Matrix Identity matrix
+ */
+ public function identity($m = null, $n = null)
+ {
+ return $this->diagonal($m, $n, 1);
+ }
+
+ /**
+ * diagonal.
+ *
+ * Generate a diagonal matrix
+ *
+ * @param int $m Row dimension
+ * @param int $n Column dimension
+ * @param mixed $c Diagonal value
+ *
+ * @return Matrix Diagonal matrix
+ */
+ public function diagonal($m = null, $n = null, $c = 1)
+ {
+ $R = new self($m, $n);
+ for ($i = 0; $i < $m; ++$i) {
+ $R->set($i, $i, $c);
+ }
+
+ return $R;
+ }
+
+ /**
+ * getMatrixByRow.
+ *
+ * Get a submatrix by row index/range
+ *
+ * @param int $i0 Initial row index
+ * @param int $iF Final row index
+ *
+ * @return Matrix Submatrix
+ */
+ public function getMatrixByRow($i0 = null, $iF = null)
+ {
+ if (is_int($i0)) {
+ if (is_int($iF)) {
+ return $this->getMatrix($i0, 0, $iF + 1, $this->n);
+ }
+
+ return $this->getMatrix($i0, 0, $i0 + 1, $this->n);
+ }
+
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ /**
+ * getMatrixByCol.
+ *
+ * Get a submatrix by column index/range
+ *
+ * @param int $j0 Initial column index
+ * @param int $jF Final column index
+ *
+ * @return Matrix Submatrix
+ */
+ public function getMatrixByCol($j0 = null, $jF = null)
+ {
+ if (is_int($j0)) {
+ if (is_int($jF)) {
+ return $this->getMatrix(0, $j0, $this->m, $jF + 1);
+ }
+
+ return $this->getMatrix(0, $j0, $this->m, $j0 + 1);
+ }
+
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ /**
+ * transpose.
+ *
+ * Tranpose matrix
+ *
+ * @return Matrix Transposed matrix
+ */
+ public function transpose()
+ {
+ $R = new self($this->n, $this->m);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $R->set($j, $i, $this->A[$i][$j]);
+ }
+ }
+
+ return $R;
+ }
+
+ // function transpose()
+
+ /**
+ * trace.
+ *
+ * Sum of diagonal elements
+ *
+ * @return float Sum of diagonal elements
+ */
+ public function trace()
+ {
+ $s = 0;
+ $n = min($this->m, $this->n);
+ for ($i = 0; $i < $n; ++$i) {
+ $s += $this->A[$i][$i];
+ }
+
+ return $s;
+ }
+
+ /**
+ * uminus.
+ *
+ * Unary minus matrix -A
+ *
+ * @return Matrix Unary minus matrix
+ */
+ public function uminus()
+ {
+ }
+
+ /**
+ * plus.
+ *
+ * A + B
+ *
+ * @return Matrix Sum
+ */
+ public function plus(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $M->set($i, $j, $M->get($i, $j) + $this->A[$i][$j]);
+ }
+ }
+
+ return $M;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * plusEquals.
+ *
+ * A = A + B
+ *
+ * @return $this
+ */
+ public function plusEquals(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $validValues = true;
+ $value = $M->get($i, $j);
+ if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
+ $this->A[$i][$j] = trim($this->A[$i][$j], '"');
+ $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
+ }
+ if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
+ $value = trim($value, '"');
+ $validValues &= StringHelper::convertToNumberIfFraction($value);
+ }
+ if ($validValues) {
+ $this->A[$i][$j] += $value;
+ } else {
+ $this->A[$i][$j] = Functions::NAN();
+ }
+ }
+ }
+
+ return $this;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * minus.
+ *
+ * A - B
+ *
+ * @return Matrix Sum
+ */
+ public function minus(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $M->set($i, $j, $M->get($i, $j) - $this->A[$i][$j]);
+ }
+ }
+
+ return $M;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * minusEquals.
+ *
+ * A = A - B
+ *
+ * @return $this
+ */
+ public function minusEquals(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $validValues = true;
+ $value = $M->get($i, $j);
+ if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
+ $this->A[$i][$j] = trim($this->A[$i][$j], '"');
+ $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
+ }
+ if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
+ $value = trim($value, '"');
+ $validValues &= StringHelper::convertToNumberIfFraction($value);
+ }
+ if ($validValues) {
+ $this->A[$i][$j] -= $value;
+ } else {
+ $this->A[$i][$j] = Functions::NAN();
+ }
+ }
+ }
+
+ return $this;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * arrayTimes.
+ *
+ * Element-by-element multiplication
+ * Cij = Aij * Bij
+ *
+ * @return Matrix Matrix Cij
+ */
+ public function arrayTimes(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $M->set($i, $j, $M->get($i, $j) * $this->A[$i][$j]);
+ }
+ }
+
+ return $M;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * arrayTimesEquals.
+ *
+ * Element-by-element multiplication
+ * Aij = Aij * Bij
+ *
+ * @return $this
+ */
+ public function arrayTimesEquals(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $validValues = true;
+ $value = $M->get($i, $j);
+ if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
+ $this->A[$i][$j] = trim($this->A[$i][$j], '"');
+ $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
+ }
+ if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
+ $value = trim($value, '"');
+ $validValues &= StringHelper::convertToNumberIfFraction($value);
+ }
+ if ($validValues) {
+ $this->A[$i][$j] *= $value;
+ } else {
+ $this->A[$i][$j] = Functions::NAN();
+ }
+ }
+ }
+
+ return $this;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * arrayRightDivide.
+ *
+ * Element-by-element right division
+ * A / B
+ *
+ * @return Matrix Division result
+ */
+ public function arrayRightDivide(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $validValues = true;
+ $value = $M->get($i, $j);
+ if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
+ $this->A[$i][$j] = trim($this->A[$i][$j], '"');
+ $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
+ }
+ if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
+ $value = trim($value, '"');
+ $validValues &= StringHelper::convertToNumberIfFraction($value);
+ }
+ if ($validValues) {
+ if ($value == 0) {
+ // Trap for Divide by Zero error
+ $M->set($i, $j, '#DIV/0!');
+ } else {
+ $M->set($i, $j, $this->A[$i][$j] / $value);
+ }
+ } else {
+ $M->set($i, $j, Functions::NAN());
+ }
+ }
+ }
+
+ return $M;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * arrayRightDivideEquals.
+ *
+ * Element-by-element right division
+ * Aij = Aij / Bij
+ *
+ * @return Matrix Matrix Aij
+ */
+ public function arrayRightDivideEquals(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $this->A[$i][$j] = $this->A[$i][$j] / $M->get($i, $j);
+ }
+ }
+
+ return $M;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * arrayLeftDivide.
+ *
+ * Element-by-element Left division
+ * A / B
+ *
+ * @return Matrix Division result
+ */
+ public function arrayLeftDivide(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $M->set($i, $j, $M->get($i, $j) / $this->A[$i][$j]);
+ }
+ }
+
+ return $M;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * arrayLeftDivideEquals.
+ *
+ * Element-by-element Left division
+ * Aij = Aij / Bij
+ *
+ * @return Matrix Matrix Aij
+ */
+ public function arrayLeftDivideEquals(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $this->A[$i][$j] = $M->get($i, $j) / $this->A[$i][$j];
+ }
+ }
+
+ return $M;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * times.
+ *
+ * Matrix multiplication
+ *
+ * @return Matrix Product
+ */
+ public function times(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $B = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+ if ($this->n == $B->m) {
+ $C = new self($this->m, $B->n);
+ for ($j = 0; $j < $B->n; ++$j) {
+ $Bcolj = [];
+ for ($k = 0; $k < $this->n; ++$k) {
+ $Bcolj[$k] = $B->A[$k][$j];
+ }
+ for ($i = 0; $i < $this->m; ++$i) {
+ $Arowi = $this->A[$i];
+ $s = 0;
+ for ($k = 0; $k < $this->n; ++$k) {
+ $s += $Arowi[$k] * $Bcolj[$k];
+ }
+ $C->A[$i][$j] = $s;
+ }
+ }
+
+ return $C;
+ }
+
+ throw new CalculationException(self::MATRIX_DIMENSION_EXCEPTION);
+ case 'array':
+ $B = new self($args[0]);
+ if ($this->n == $B->m) {
+ $C = new self($this->m, $B->n);
+ for ($i = 0; $i < $C->m; ++$i) {
+ for ($j = 0; $j < $C->n; ++$j) {
+ $s = '0';
+ for ($k = 0; $k < $C->n; ++$k) {
+ $s += $this->A[$i][$k] * $B->A[$k][$j];
+ }
+ $C->A[$i][$j] = $s;
+ }
+ }
+
+ return $C;
+ }
+
+ throw new CalculationException(self::MATRIX_DIMENSION_EXCEPTION);
+ case 'integer':
+ $C = new self($this->A);
+ for ($i = 0; $i < $C->m; ++$i) {
+ for ($j = 0; $j < $C->n; ++$j) {
+ $C->A[$i][$j] *= $args[0];
+ }
+ }
+
+ return $C;
+ case 'double':
+ $C = new self($this->m, $this->n);
+ for ($i = 0; $i < $C->m; ++$i) {
+ for ($j = 0; $j < $C->n; ++$j) {
+ $C->A[$i][$j] = $args[0] * $this->A[$i][$j];
+ }
+ }
+
+ return $C;
+ case 'float':
+ $C = new self($this->A);
+ for ($i = 0; $i < $C->m; ++$i) {
+ for ($j = 0; $j < $C->n; ++$j) {
+ $C->A[$i][$j] *= $args[0];
+ }
+ }
+
+ return $C;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+ } else {
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+ }
+
+ /**
+ * power.
+ *
+ * A = A ^ B
+ *
+ * @return $this
+ */
+ public function power(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $validValues = true;
+ $value = $M->get($i, $j);
+ if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
+ $this->A[$i][$j] = trim($this->A[$i][$j], '"');
+ $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
+ }
+ if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
+ $value = trim($value, '"');
+ $validValues &= StringHelper::convertToNumberIfFraction($value);
+ }
+ if ($validValues) {
+ $this->A[$i][$j] = $this->A[$i][$j] ** $value;
+ } else {
+ $this->A[$i][$j] = Functions::NAN();
+ }
+ }
+ }
+
+ return $this;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * concat.
+ *
+ * A = A & B
+ *
+ * @return $this
+ */
+ public function concat(...$args)
+ {
+ if (count($args) > 0) {
+ $match = implode(',', array_map('gettype', $args));
+
+ switch ($match) {
+ case 'object':
+ if ($args[0] instanceof self) {
+ $M = $args[0];
+ } else {
+ throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
+ }
+
+ break;
+ case 'array':
+ $M = new self($args[0]);
+
+ break;
+ default:
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+
+ break;
+ }
+ $this->checkMatrixDimensions($M);
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $this->A[$i][$j] = trim($this->A[$i][$j], '"') . trim($M->get($i, $j), '"');
+ }
+ }
+
+ return $this;
+ }
+
+ throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
+ }
+
+ /**
+ * Solve A*X = B.
+ *
+ * @param Matrix $B Right hand side
+ *
+ * @return Matrix ... Solution if A is square, least squares solution otherwise
+ */
+ public function solve($B)
+ {
+ if ($this->m == $this->n) {
+ $LU = new LUDecomposition($this);
+
+ return $LU->solve($B);
+ }
+ $QR = new QRDecomposition($this);
+
+ return $QR->solve($B);
+ }
+
+ /**
+ * Matrix inverse or pseudoinverse.
+ *
+ * @return Matrix ... Inverse(A) if A is square, pseudoinverse otherwise.
+ */
+ public function inverse()
+ {
+ return $this->solve($this->identity($this->m, $this->m));
+ }
+
+ /**
+ * det.
+ *
+ * Calculate determinant
+ *
+ * @return float Determinant
+ */
+ public function det()
+ {
+ $L = new LUDecomposition($this);
+
+ return $L->det();
+ }
+}
diff --git a/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/QRDecomposition.php b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/QRDecomposition.php new file mode 100644 index 0000000..373c074 --- /dev/null +++ b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/QRDecomposition.php @@ -0,0 +1,249 @@ +<?php
+
+namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
+
+use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
+
+/**
+ * For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
+ * orthogonal matrix Q and an n-by-n upper triangular matrix R so that
+ * A = Q*R.
+ *
+ * The QR decompostion always exists, even if the matrix does not have
+ * full rank, so the constructor will never fail. The primary use of the
+ * QR decomposition is in the least squares solution of nonsquare systems
+ * of simultaneous linear equations. This will fail if isFullRank()
+ * returns false.
+ *
+ * @author Paul Meagher
+ *
+ * @version 1.1
+ */
+class QRDecomposition
+{
+ const MATRIX_RANK_EXCEPTION = 'Can only perform operation on full-rank matrix.';
+
+ /**
+ * Array for internal storage of decomposition.
+ *
+ * @var array
+ */
+ private $QR = [];
+
+ /**
+ * Row dimension.
+ *
+ * @var int
+ */
+ private $m;
+
+ /**
+ * Column dimension.
+ *
+ * @var int
+ */
+ private $n;
+
+ /**
+ * Array for internal storage of diagonal of R.
+ *
+ * @var array
+ */
+ private $Rdiag = [];
+
+ /**
+ * QR Decomposition computed by Householder reflections.
+ *
+ * @param matrix $A Rectangular matrix
+ */
+ public function __construct($A)
+ {
+ if ($A instanceof Matrix) {
+ // Initialize.
+ $this->QR = $A->getArray();
+ $this->m = $A->getRowDimension();
+ $this->n = $A->getColumnDimension();
+ // Main loop.
+ for ($k = 0; $k < $this->n; ++$k) {
+ // Compute 2-norm of k-th column without under/overflow.
+ $nrm = 0.0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $nrm = hypo($nrm, $this->QR[$i][$k]);
+ }
+ if ($nrm != 0.0) {
+ // Form k-th Householder vector.
+ if ($this->QR[$k][$k] < 0) {
+ $nrm = -$nrm;
+ }
+ for ($i = $k; $i < $this->m; ++$i) {
+ $this->QR[$i][$k] /= $nrm;
+ }
+ $this->QR[$k][$k] += 1.0;
+ // Apply transformation to remaining columns.
+ for ($j = $k + 1; $j < $this->n; ++$j) {
+ $s = 0.0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $s += $this->QR[$i][$k] * $this->QR[$i][$j];
+ }
+ $s = -$s / $this->QR[$k][$k];
+ for ($i = $k; $i < $this->m; ++$i) {
+ $this->QR[$i][$j] += $s * $this->QR[$i][$k];
+ }
+ }
+ }
+ $this->Rdiag[$k] = -$nrm;
+ }
+ } else {
+ throw new CalculationException(Matrix::ARGUMENT_TYPE_EXCEPTION);
+ }
+ }
+
+ // function __construct()
+
+ /**
+ * Is the matrix full rank?
+ *
+ * @return bool true if R, and hence A, has full rank, else false
+ */
+ public function isFullRank()
+ {
+ for ($j = 0; $j < $this->n; ++$j) {
+ if ($this->Rdiag[$j] == 0) {
+ return false;
+ }
+ }
+
+ return true;
+ }
+
+ // function isFullRank()
+
+ /**
+ * Return the Householder vectors.
+ *
+ * @return Matrix Lower trapezoidal matrix whose columns define the reflections
+ */
+ public function getH()
+ {
+ $H = [];
+ for ($i = 0; $i < $this->m; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ if ($i >= $j) {
+ $H[$i][$j] = $this->QR[$i][$j];
+ } else {
+ $H[$i][$j] = 0.0;
+ }
+ }
+ }
+
+ return new Matrix($H);
+ }
+
+ // function getH()
+
+ /**
+ * Return the upper triangular factor.
+ *
+ * @return Matrix upper triangular factor
+ */
+ public function getR()
+ {
+ $R = [];
+ for ($i = 0; $i < $this->n; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ if ($i < $j) {
+ $R[$i][$j] = $this->QR[$i][$j];
+ } elseif ($i == $j) {
+ $R[$i][$j] = $this->Rdiag[$i];
+ } else {
+ $R[$i][$j] = 0.0;
+ }
+ }
+ }
+
+ return new Matrix($R);
+ }
+
+ // function getR()
+
+ /**
+ * Generate and return the (economy-sized) orthogonal factor.
+ *
+ * @return Matrix orthogonal factor
+ */
+ public function getQ()
+ {
+ $Q = [];
+ for ($k = $this->n - 1; $k >= 0; --$k) {
+ for ($i = 0; $i < $this->m; ++$i) {
+ $Q[$i][$k] = 0.0;
+ }
+ $Q[$k][$k] = 1.0;
+ for ($j = $k; $j < $this->n; ++$j) {
+ if ($this->QR[$k][$k] != 0) {
+ $s = 0.0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $s += $this->QR[$i][$k] * $Q[$i][$j];
+ }
+ $s = -$s / $this->QR[$k][$k];
+ for ($i = $k; $i < $this->m; ++$i) {
+ $Q[$i][$j] += $s * $this->QR[$i][$k];
+ }
+ }
+ }
+ }
+
+ return new Matrix($Q);
+ }
+
+ // function getQ()
+
+ /**
+ * Least squares solution of A*X = B.
+ *
+ * @param Matrix $B a Matrix with as many rows as A and any number of columns
+ *
+ * @return Matrix matrix that minimizes the two norm of Q*R*X-B
+ */
+ public function solve($B)
+ {
+ if ($B->getRowDimension() == $this->m) {
+ if ($this->isFullRank()) {
+ // Copy right hand side
+ $nx = $B->getColumnDimension();
+ $X = $B->getArrayCopy();
+ // Compute Y = transpose(Q)*B
+ for ($k = 0; $k < $this->n; ++$k) {
+ for ($j = 0; $j < $nx; ++$j) {
+ $s = 0.0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $s += $this->QR[$i][$k] * $X[$i][$j];
+ }
+ $s = -$s / $this->QR[$k][$k];
+ for ($i = $k; $i < $this->m; ++$i) {
+ $X[$i][$j] += $s * $this->QR[$i][$k];
+ }
+ }
+ }
+ // Solve R*X = Y;
+ for ($k = $this->n - 1; $k >= 0; --$k) {
+ for ($j = 0; $j < $nx; ++$j) {
+ $X[$k][$j] /= $this->Rdiag[$k];
+ }
+ for ($i = 0; $i < $k; ++$i) {
+ for ($j = 0; $j < $nx; ++$j) {
+ $X[$i][$j] -= $X[$k][$j] * $this->QR[$i][$k];
+ }
+ }
+ }
+ $X = new Matrix($X);
+
+ return $X->getMatrix(0, $this->n - 1, 0, $nx);
+ }
+
+ throw new CalculationException(self::MATRIX_RANK_EXCEPTION);
+ }
+
+ throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
+ }
+}
diff --git a/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/SingularValueDecomposition.php b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/SingularValueDecomposition.php new file mode 100644 index 0000000..a86bf08 --- /dev/null +++ b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/SingularValueDecomposition.php @@ -0,0 +1,528 @@ +<?php
+
+namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
+
+/**
+ * For an m-by-n matrix A with m >= n, the singular value decomposition is
+ * an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
+ * an n-by-n orthogonal matrix V so that A = U*S*V'.
+ *
+ * The singular values, sigma[$k] = S[$k][$k], are ordered so that
+ * sigma[0] >= sigma[1] >= ... >= sigma[n-1].
+ *
+ * The singular value decompostion always exists, so the constructor will
+ * never fail. The matrix condition number and the effective numerical
+ * rank can be computed from this decomposition.
+ *
+ * @author Paul Meagher
+ *
+ * @version 1.1
+ */
+class SingularValueDecomposition
+{
+ /**
+ * Internal storage of U.
+ *
+ * @var array
+ */
+ private $U = [];
+
+ /**
+ * Internal storage of V.
+ *
+ * @var array
+ */
+ private $V = [];
+
+ /**
+ * Internal storage of singular values.
+ *
+ * @var array
+ */
+ private $s = [];
+
+ /**
+ * Row dimension.
+ *
+ * @var int
+ */
+ private $m;
+
+ /**
+ * Column dimension.
+ *
+ * @var int
+ */
+ private $n;
+
+ /**
+ * Construct the singular value decomposition.
+ *
+ * Derived from LINPACK code.
+ *
+ * @param mixed $Arg Rectangular matrix
+ */
+ public function __construct($Arg)
+ {
+ // Initialize.
+ $A = $Arg->getArrayCopy();
+ $this->m = $Arg->getRowDimension();
+ $this->n = $Arg->getColumnDimension();
+ $nu = min($this->m, $this->n);
+ $e = [];
+ $work = [];
+ $wantu = true;
+ $wantv = true;
+ $nct = min($this->m - 1, $this->n);
+ $nrt = max(0, min($this->n - 2, $this->m));
+
+ // Reduce A to bidiagonal form, storing the diagonal elements
+ // in s and the super-diagonal elements in e.
+ $kMax = max($nct, $nrt);
+ for ($k = 0; $k < $kMax; ++$k) {
+ if ($k < $nct) {
+ // Compute the transformation for the k-th column and
+ // place the k-th diagonal in s[$k].
+ // Compute 2-norm of k-th column without under/overflow.
+ $this->s[$k] = 0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
+ }
+ if ($this->s[$k] != 0.0) {
+ if ($A[$k][$k] < 0.0) {
+ $this->s[$k] = -$this->s[$k];
+ }
+ for ($i = $k; $i < $this->m; ++$i) {
+ $A[$i][$k] /= $this->s[$k];
+ }
+ $A[$k][$k] += 1.0;
+ }
+ $this->s[$k] = -$this->s[$k];
+ }
+
+ for ($j = $k + 1; $j < $this->n; ++$j) {
+ if (($k < $nct) & ($this->s[$k] != 0.0)) {
+ // Apply the transformation.
+ $t = 0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $t += $A[$i][$k] * $A[$i][$j];
+ }
+ $t = -$t / $A[$k][$k];
+ for ($i = $k; $i < $this->m; ++$i) {
+ $A[$i][$j] += $t * $A[$i][$k];
+ }
+ // Place the k-th row of A into e for the
+ // subsequent calculation of the row transformation.
+ $e[$j] = $A[$k][$j];
+ }
+ }
+
+ if ($wantu && ($k < $nct)) {
+ // Place the transformation in U for subsequent back
+ // multiplication.
+ for ($i = $k; $i < $this->m; ++$i) {
+ $this->U[$i][$k] = $A[$i][$k];
+ }
+ }
+
+ if ($k < $nrt) {
+ // Compute the k-th row transformation and place the
+ // k-th super-diagonal in e[$k].
+ // Compute 2-norm without under/overflow.
+ $e[$k] = 0;
+ for ($i = $k + 1; $i < $this->n; ++$i) {
+ $e[$k] = hypo($e[$k], $e[$i]);
+ }
+ if ($e[$k] != 0.0) {
+ if ($e[$k + 1] < 0.0) {
+ $e[$k] = -$e[$k];
+ }
+ for ($i = $k + 1; $i < $this->n; ++$i) {
+ $e[$i] /= $e[$k];
+ }
+ $e[$k + 1] += 1.0;
+ }
+ $e[$k] = -$e[$k];
+ if (($k + 1 < $this->m) && ($e[$k] != 0.0)) {
+ // Apply the transformation.
+ for ($i = $k + 1; $i < $this->m; ++$i) {
+ $work[$i] = 0.0;
+ }
+ for ($j = $k + 1; $j < $this->n; ++$j) {
+ for ($i = $k + 1; $i < $this->m; ++$i) {
+ $work[$i] += $e[$j] * $A[$i][$j];
+ }
+ }
+ for ($j = $k + 1; $j < $this->n; ++$j) {
+ $t = -$e[$j] / $e[$k + 1];
+ for ($i = $k + 1; $i < $this->m; ++$i) {
+ $A[$i][$j] += $t * $work[$i];
+ }
+ }
+ }
+ if ($wantv) {
+ // Place the transformation in V for subsequent
+ // back multiplication.
+ for ($i = $k + 1; $i < $this->n; ++$i) {
+ $this->V[$i][$k] = $e[$i];
+ }
+ }
+ }
+ }
+
+ // Set up the final bidiagonal matrix or order p.
+ $p = min($this->n, $this->m + 1);
+ if ($nct < $this->n) {
+ $this->s[$nct] = $A[$nct][$nct];
+ }
+ if ($this->m < $p) {
+ $this->s[$p - 1] = 0.0;
+ }
+ if ($nrt + 1 < $p) {
+ $e[$nrt] = $A[$nrt][$p - 1];
+ }
+ $e[$p - 1] = 0.0;
+ // If required, generate U.
+ if ($wantu) {
+ for ($j = $nct; $j < $nu; ++$j) {
+ for ($i = 0; $i < $this->m; ++$i) {
+ $this->U[$i][$j] = 0.0;
+ }
+ $this->U[$j][$j] = 1.0;
+ }
+ for ($k = $nct - 1; $k >= 0; --$k) {
+ if ($this->s[$k] != 0.0) {
+ for ($j = $k + 1; $j < $nu; ++$j) {
+ $t = 0;
+ for ($i = $k; $i < $this->m; ++$i) {
+ $t += $this->U[$i][$k] * $this->U[$i][$j];
+ }
+ $t = -$t / $this->U[$k][$k];
+ for ($i = $k; $i < $this->m; ++$i) {
+ $this->U[$i][$j] += $t * $this->U[$i][$k];
+ }
+ }
+ for ($i = $k; $i < $this->m; ++$i) {
+ $this->U[$i][$k] = -$this->U[$i][$k];
+ }
+ $this->U[$k][$k] = 1.0 + $this->U[$k][$k];
+ for ($i = 0; $i < $k - 1; ++$i) {
+ $this->U[$i][$k] = 0.0;
+ }
+ } else {
+ for ($i = 0; $i < $this->m; ++$i) {
+ $this->U[$i][$k] = 0.0;
+ }
+ $this->U[$k][$k] = 1.0;
+ }
+ }
+ }
+
+ // If required, generate V.
+ if ($wantv) {
+ for ($k = $this->n - 1; $k >= 0; --$k) {
+ if (($k < $nrt) && ($e[$k] != 0.0)) {
+ for ($j = $k + 1; $j < $nu; ++$j) {
+ $t = 0;
+ for ($i = $k + 1; $i < $this->n; ++$i) {
+ $t += $this->V[$i][$k] * $this->V[$i][$j];
+ }
+ $t = -$t / $this->V[$k + 1][$k];
+ for ($i = $k + 1; $i < $this->n; ++$i) {
+ $this->V[$i][$j] += $t * $this->V[$i][$k];
+ }
+ }
+ }
+ for ($i = 0; $i < $this->n; ++$i) {
+ $this->V[$i][$k] = 0.0;
+ }
+ $this->V[$k][$k] = 1.0;
+ }
+ }
+
+ // Main iteration loop for the singular values.
+ $pp = $p - 1;
+ $iter = 0;
+ $eps = 2.0 ** (-52.0);
+
+ while ($p > 0) {
+ // Here is where a test for too many iterations would go.
+ // This section of the program inspects for negligible
+ // elements in the s and e arrays. On completion the
+ // variables kase and k are set as follows:
+ // kase = 1 if s(p) and e[k-1] are negligible and k<p
+ // kase = 2 if s(k) is negligible and k<p
+ // kase = 3 if e[k-1] is negligible, k<p, and
+ // s(k), ..., s(p) are not negligible (qr step).
+ // kase = 4 if e(p-1) is negligible (convergence).
+ for ($k = $p - 2; $k >= -1; --$k) {
+ if ($k == -1) {
+ break;
+ }
+ if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k + 1]))) {
+ $e[$k] = 0.0;
+
+ break;
+ }
+ }
+ if ($k == $p - 2) {
+ $kase = 4;
+ } else {
+ for ($ks = $p - 1; $ks >= $k; --$ks) {
+ if ($ks == $k) {
+ break;
+ }
+ $t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks - 1]) : 0.);
+ if (abs($this->s[$ks]) <= $eps * $t) {
+ $this->s[$ks] = 0.0;
+
+ break;
+ }
+ }
+ if ($ks == $k) {
+ $kase = 3;
+ } elseif ($ks == $p - 1) {
+ $kase = 1;
+ } else {
+ $kase = 2;
+ $k = $ks;
+ }
+ }
+ ++$k;
+
+ // Perform the task indicated by kase.
+ switch ($kase) {
+ // Deflate negligible s(p).
+ case 1:
+ $f = $e[$p - 2];
+ $e[$p - 2] = 0.0;
+ for ($j = $p - 2; $j >= $k; --$j) {
+ $t = hypo($this->s[$j], $f);
+ $cs = $this->s[$j] / $t;
+ $sn = $f / $t;
+ $this->s[$j] = $t;
+ if ($j != $k) {
+ $f = -$sn * $e[$j - 1];
+ $e[$j - 1] = $cs * $e[$j - 1];
+ }
+ if ($wantv) {
+ for ($i = 0; $i < $this->n; ++$i) {
+ $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p - 1];
+ $this->V[$i][$p - 1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p - 1];
+ $this->V[$i][$j] = $t;
+ }
+ }
+ }
+
+ break;
+ // Split at negligible s(k).
+ case 2:
+ $f = $e[$k - 1];
+ $e[$k - 1] = 0.0;
+ for ($j = $k; $j < $p; ++$j) {
+ $t = hypo($this->s[$j], $f);
+ $cs = $this->s[$j] / $t;
+ $sn = $f / $t;
+ $this->s[$j] = $t;
+ $f = -$sn * $e[$j];
+ $e[$j] = $cs * $e[$j];
+ if ($wantu) {
+ for ($i = 0; $i < $this->m; ++$i) {
+ $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k - 1];
+ $this->U[$i][$k - 1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k - 1];
+ $this->U[$i][$j] = $t;
+ }
+ }
+ }
+
+ break;
+ // Perform one qr step.
+ case 3:
+ // Calculate the shift.
+ $scale = max(max(max(max(abs($this->s[$p - 1]), abs($this->s[$p - 2])), abs($e[$p - 2])), abs($this->s[$k])), abs($e[$k]));
+ $sp = $this->s[$p - 1] / $scale;
+ $spm1 = $this->s[$p - 2] / $scale;
+ $epm1 = $e[$p - 2] / $scale;
+ $sk = $this->s[$k] / $scale;
+ $ek = $e[$k] / $scale;
+ $b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
+ $c = ($sp * $epm1) * ($sp * $epm1);
+ $shift = 0.0;
+ if (($b != 0.0) || ($c != 0.0)) {
+ $shift = sqrt($b * $b + $c);
+ if ($b < 0.0) {
+ $shift = -$shift;
+ }
+ $shift = $c / ($b + $shift);
+ }
+ $f = ($sk + $sp) * ($sk - $sp) + $shift;
+ $g = $sk * $ek;
+ // Chase zeros.
+ for ($j = $k; $j < $p - 1; ++$j) {
+ $t = hypo($f, $g);
+ $cs = $f / $t;
+ $sn = $g / $t;
+ if ($j != $k) {
+ $e[$j - 1] = $t;
+ }
+ $f = $cs * $this->s[$j] + $sn * $e[$j];
+ $e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
+ $g = $sn * $this->s[$j + 1];
+ $this->s[$j + 1] = $cs * $this->s[$j + 1];
+ if ($wantv) {
+ for ($i = 0; $i < $this->n; ++$i) {
+ $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j + 1];
+ $this->V[$i][$j + 1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j + 1];
+ $this->V[$i][$j] = $t;
+ }
+ }
+ $t = hypo($f, $g);
+ $cs = $f / $t;
+ $sn = $g / $t;
+ $this->s[$j] = $t;
+ $f = $cs * $e[$j] + $sn * $this->s[$j + 1];
+ $this->s[$j + 1] = -$sn * $e[$j] + $cs * $this->s[$j + 1];
+ $g = $sn * $e[$j + 1];
+ $e[$j + 1] = $cs * $e[$j + 1];
+ if ($wantu && ($j < $this->m - 1)) {
+ for ($i = 0; $i < $this->m; ++$i) {
+ $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j + 1];
+ $this->U[$i][$j + 1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j + 1];
+ $this->U[$i][$j] = $t;
+ }
+ }
+ }
+ $e[$p - 2] = $f;
+ $iter = $iter + 1;
+
+ break;
+ // Convergence.
+ case 4:
+ // Make the singular values positive.
+ if ($this->s[$k] <= 0.0) {
+ $this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);
+ if ($wantv) {
+ for ($i = 0; $i <= $pp; ++$i) {
+ $this->V[$i][$k] = -$this->V[$i][$k];
+ }
+ }
+ }
+ // Order the singular values.
+ while ($k < $pp) {
+ if ($this->s[$k] >= $this->s[$k + 1]) {
+ break;
+ }
+ $t = $this->s[$k];
+ $this->s[$k] = $this->s[$k + 1];
+ $this->s[$k + 1] = $t;
+ if ($wantv && ($k < $this->n - 1)) {
+ for ($i = 0; $i < $this->n; ++$i) {
+ $t = $this->V[$i][$k + 1];
+ $this->V[$i][$k + 1] = $this->V[$i][$k];
+ $this->V[$i][$k] = $t;
+ }
+ }
+ if ($wantu && ($k < $this->m - 1)) {
+ for ($i = 0; $i < $this->m; ++$i) {
+ $t = $this->U[$i][$k + 1];
+ $this->U[$i][$k + 1] = $this->U[$i][$k];
+ $this->U[$i][$k] = $t;
+ }
+ }
+ ++$k;
+ }
+ $iter = 0;
+ --$p;
+
+ break;
+ } // end switch
+ } // end while
+ }
+
+ /**
+ * Return the left singular vectors.
+ *
+ * @return Matrix U
+ */
+ public function getU()
+ {
+ return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
+ }
+
+ /**
+ * Return the right singular vectors.
+ *
+ * @return Matrix V
+ */
+ public function getV()
+ {
+ return new Matrix($this->V);
+ }
+
+ /**
+ * Return the one-dimensional array of singular values.
+ *
+ * @return array diagonal of S
+ */
+ public function getSingularValues()
+ {
+ return $this->s;
+ }
+
+ /**
+ * Return the diagonal matrix of singular values.
+ *
+ * @return Matrix S
+ */
+ public function getS()
+ {
+ for ($i = 0; $i < $this->n; ++$i) {
+ for ($j = 0; $j < $this->n; ++$j) {
+ $S[$i][$j] = 0.0;
+ }
+ $S[$i][$i] = $this->s[$i];
+ }
+
+ return new Matrix($S);
+ }
+
+ /**
+ * Two norm.
+ *
+ * @return float max(S)
+ */
+ public function norm2()
+ {
+ return $this->s[0];
+ }
+
+ /**
+ * Two norm condition number.
+ *
+ * @return float max(S)/min(S)
+ */
+ public function cond()
+ {
+ return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
+ }
+
+ /**
+ * Effective numerical matrix rank.
+ *
+ * @return int Number of nonnegligible singular values
+ */
+ public function rank()
+ {
+ $eps = 2.0 ** (-52.0);
+ $tol = max($this->m, $this->n) * $this->s[0] * $eps;
+ $r = 0;
+ $iMax = count($this->s);
+ for ($i = 0; $i < $iMax; ++$i) {
+ if ($this->s[$i] > $tol) {
+ ++$r;
+ }
+ }
+
+ return $r;
+ }
+}
diff --git a/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/utils/Maths.php b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/utils/Maths.php new file mode 100644 index 0000000..8e7c92b --- /dev/null +++ b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/utils/Maths.php @@ -0,0 +1,31 @@ +<?php
+
+/**
+ * Pythagorean Theorem:.
+ *
+ * a = 3
+ * b = 4
+ * r = sqrt(square(a) + square(b))
+ * r = 5
+ *
+ * r = sqrt(a^2 + b^2) without under/overflow.
+ *
+ * @param mixed $a
+ * @param mixed $b
+ *
+ * @return float
+ */
+function hypo($a, $b)
+{
+ if (abs($a) > abs($b)) {
+ $r = $b / $a;
+ $r = abs($a) * sqrt(1 + $r * $r);
+ } elseif ($b != 0) {
+ $r = $a / $b;
+ $r = abs($b) * sqrt(1 + $r * $r);
+ } else {
+ $r = 0.0;
+ }
+
+ return $r;
+}
|