[Summary view]

Differences Between: [Versions 310 and 311] [Versions 311 and 400] [Versions 311 and 401] [Versions 39 and 311]

   1  <?php
   2  
   3  namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
   4  
   5  use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
   6  use PhpOffice\PhpSpreadsheet\Calculation\Functions;
   7  use PhpOffice\PhpSpreadsheet\Shared\StringHelper;
   8  
   9  /**
  10   * Matrix class.
  11   *
  12   * @author Paul Meagher
  13   * @author Michael Bommarito
  14   * @author Lukasz Karapuda
  15   * @author Bartek Matosiuk
  16   *
  17   * @version 1.8
  18   *
  19   * @see https://math.nist.gov/javanumerics/jama/
  20   */
  21  class Matrix
  22  {
  23      const POLYMORPHIC_ARGUMENT_EXCEPTION = 'Invalid argument pattern for polymorphic function.';
  24      const ARGUMENT_TYPE_EXCEPTION = 'Invalid argument type.';
  25      const ARGUMENT_BOUNDS_EXCEPTION = 'Invalid argument range.';
  26      const MATRIX_DIMENSION_EXCEPTION = 'Matrix dimensions are not equal.';
  27      const ARRAY_LENGTH_EXCEPTION = 'Array length must be a multiple of m.';
  28      const MATRIX_SPD_EXCEPTION = 'Can only perform operation on symmetric positive definite matrix.';
  29  
  30      /**
  31       * Matrix storage.
  32       *
  33       * @var array
  34       */
  35      public $A = [];
  36  
  37      /**
  38       * Matrix row dimension.
  39       *
  40       * @var int
  41       */
  42      private $m;
  43  
  44      /**
  45       * Matrix column dimension.
  46       *
  47       * @var int
  48       */
  49      private $n;
  50  
  51      /**
  52       * Polymorphic constructor.
  53       *
  54       * 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.
  55       */
  56      public function __construct(...$args)
  57      {
  58          if (count($args) > 0) {
  59              $match = implode(',', array_map('gettype', $args));
  60  
  61              switch ($match) {
  62                  //Rectangular matrix - m x n initialized from 2D array
  63                  case 'array':
  64                      $this->m = count($args[0]);
  65                      $this->n = count($args[0][0]);
  66                      $this->A = $args[0];
  67  
  68                      break;
  69                  //Square matrix - n x n
  70                  case 'integer':
  71                      $this->m = $args[0];
  72                      $this->n = $args[0];
  73                      $this->A = array_fill(0, $this->m, array_fill(0, $this->n, 0));
  74  
  75                      break;
  76                  //Rectangular matrix - m x n
  77                  case 'integer,integer':
  78                      $this->m = $args[0];
  79                      $this->n = $args[1];
  80                      $this->A = array_fill(0, $this->m, array_fill(0, $this->n, 0));
  81  
  82                      break;
  83                  //Rectangular matrix - m x n initialized from packed array
  84                  case 'array,integer':
  85                      $this->m = $args[1];
  86                      if ($this->m != 0) {
  87                          $this->n = count($args[0]) / $this->m;
  88                      } else {
  89                          $this->n = 0;
  90                      }
  91                      if (($this->m * $this->n) == count($args[0])) {
  92                          for ($i = 0; $i < $this->m; ++$i) {
  93                              for ($j = 0; $j < $this->n; ++$j) {
  94                                  $this->A[$i][$j] = $args[0][$i + $j * $this->m];
  95                              }
  96                          }
  97                      } else {
  98                          throw new CalculationException(self::ARRAY_LENGTH_EXCEPTION);
  99                      }
 100  
 101                      break;
 102                  default:
 103                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 104  
 105                      break;
 106              }
 107          } else {
 108              throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 109          }
 110      }
 111  
 112      /**
 113       * getArray.
 114       *
 115       * @return array Matrix array
 116       */
 117      public function getArray()
 118      {
 119          return $this->A;
 120      }
 121  
 122      /**
 123       * getRowDimension.
 124       *
 125       * @return int Row dimension
 126       */
 127      public function getRowDimension()
 128      {
 129          return $this->m;
 130      }
 131  
 132      /**
 133       * getColumnDimension.
 134       *
 135       * @return int Column dimension
 136       */
 137      public function getColumnDimension()
 138      {
 139          return $this->n;
 140      }
 141  
 142      /**
 143       * get.
 144       *
 145       * Get the i,j-th element of the matrix.
 146       *
 147       * @param int $i Row position
 148       * @param int $j Column position
 149       *
 150       * @return mixed Element (int/float/double)
 151       */
 152      public function get($i = null, $j = null)
 153      {
 154          return $this->A[$i][$j];
 155      }
 156  
 157      /**
 158       * getMatrix.
 159       *
 160       *    Get a submatrix
 161       *
 162       * @return Matrix Submatrix
 163       */
 164      public function getMatrix(...$args)
 165      {
 166          if (count($args) > 0) {
 167              $match = implode(',', array_map('gettype', $args));
 168  
 169              switch ($match) {
 170                  //A($i0...; $j0...)
 171                  case 'integer,integer':
 172                      [$i0, $j0] = $args;
 173                      if ($i0 >= 0) {
 174                          $m = $this->m - $i0;
 175                      } else {
 176                          throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
 177                      }
 178                      if ($j0 >= 0) {
 179                          $n = $this->n - $j0;
 180                      } else {
 181                          throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
 182                      }
 183                      $R = new self($m, $n);
 184                      for ($i = $i0; $i < $this->m; ++$i) {
 185                          for ($j = $j0; $j < $this->n; ++$j) {
 186                              $R->set($i, $j, $this->A[$i][$j]);
 187                          }
 188                      }
 189  
 190                      return $R;
 191  
 192                      break;
 193                  //A($i0...$iF; $j0...$jF)
 194                  case 'integer,integer,integer,integer':
 195                      [$i0, $iF, $j0, $jF] = $args;
 196                      if (($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0)) {
 197                          $m = $iF - $i0;
 198                      } else {
 199                          throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
 200                      }
 201                      if (($jF > $j0) && ($this->n >= $jF) && ($j0 >= 0)) {
 202                          $n = $jF - $j0;
 203                      } else {
 204                          throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
 205                      }
 206                      $R = new self($m + 1, $n + 1);
 207                      for ($i = $i0; $i <= $iF; ++$i) {
 208                          for ($j = $j0; $j <= $jF; ++$j) {
 209                              $R->set($i - $i0, $j - $j0, $this->A[$i][$j]);
 210                          }
 211                      }
 212  
 213                      return $R;
 214  
 215                      break;
 216                  //$R = array of row indices; $C = array of column indices
 217                  case 'array,array':
 218                      [$RL, $CL] = $args;
 219                      if (count($RL) > 0) {
 220                          $m = count($RL);
 221                      } else {
 222                          throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
 223                      }
 224                      if (count($CL) > 0) {
 225                          $n = count($CL);
 226                      } else {
 227                          throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
 228                      }
 229                      $R = new self($m, $n);
 230                      for ($i = 0; $i < $m; ++$i) {
 231                          for ($j = 0; $j < $n; ++$j) {
 232                              $R->set($i, $j, $this->A[$RL[$i]][$CL[$j]]);
 233                          }
 234                      }
 235  
 236                      return $R;
 237  
 238                      break;
 239                  //A($i0...$iF); $CL = array of column indices
 240                  case 'integer,integer,array':
 241                      [$i0, $iF, $CL] = $args;
 242                      if (($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0)) {
 243                          $m = $iF - $i0;
 244                      } else {
 245                          throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
 246                      }
 247                      if (count($CL) > 0) {
 248                          $n = count($CL);
 249                      } else {
 250                          throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
 251                      }
 252                      $R = new self($m, $n);
 253                      for ($i = $i0; $i < $iF; ++$i) {
 254                          for ($j = 0; $j < $n; ++$j) {
 255                              $R->set($i - $i0, $j, $this->A[$i][$CL[$j]]);
 256                          }
 257                      }
 258  
 259                      return $R;
 260  
 261                      break;
 262                  //$RL = array of row indices
 263                  case 'array,integer,integer':
 264                      [$RL, $j0, $jF] = $args;
 265                      if (count($RL) > 0) {
 266                          $m = count($RL);
 267                      } else {
 268                          throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
 269                      }
 270                      if (($jF >= $j0) && ($this->n >= $jF) && ($j0 >= 0)) {
 271                          $n = $jF - $j0;
 272                      } else {
 273                          throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
 274                      }
 275                      $R = new self($m, $n + 1);
 276                      for ($i = 0; $i < $m; ++$i) {
 277                          for ($j = $j0; $j <= $jF; ++$j) {
 278                              $R->set($i, $j - $j0, $this->A[$RL[$i]][$j]);
 279                          }
 280                      }
 281  
 282                      return $R;
 283  
 284                      break;
 285                  default:
 286                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 287  
 288                      break;
 289              }
 290          } else {
 291              throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 292          }
 293      }
 294  
 295      /**
 296       * checkMatrixDimensions.
 297       *
 298       *    Is matrix B the same size?
 299       *
 300       * @param Matrix $B Matrix B
 301       *
 302       * @return bool
 303       */
 304      public function checkMatrixDimensions($B = null)
 305      {
 306          if ($B instanceof self) {
 307              if (($this->m == $B->getRowDimension()) && ($this->n == $B->getColumnDimension())) {
 308                  return true;
 309              }
 310  
 311              throw new CalculationException(self::MATRIX_DIMENSION_EXCEPTION);
 312          }
 313  
 314          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 315      }
 316  
 317      //    function checkMatrixDimensions()
 318  
 319      /**
 320       * set.
 321       *
 322       * Set the i,j-th element of the matrix.
 323       *
 324       * @param int $i Row position
 325       * @param int $j Column position
 326       * @param mixed $c Int/float/double value
 327       *
 328       * @return mixed Element (int/float/double)
 329       */
 330      public function set($i = null, $j = null, $c = null)
 331      {
 332          // Optimized set version just has this
 333          $this->A[$i][$j] = $c;
 334      }
 335  
 336      //    function set()
 337  
 338      /**
 339       * identity.
 340       *
 341       * Generate an identity matrix.
 342       *
 343       * @param int $m Row dimension
 344       * @param int $n Column dimension
 345       *
 346       * @return Matrix Identity matrix
 347       */
 348      public function identity($m = null, $n = null)
 349      {
 350          return $this->diagonal($m, $n, 1);
 351      }
 352  
 353      /**
 354       * diagonal.
 355       *
 356       *    Generate a diagonal matrix
 357       *
 358       * @param int $m Row dimension
 359       * @param int $n Column dimension
 360       * @param mixed $c Diagonal value
 361       *
 362       * @return Matrix Diagonal matrix
 363       */
 364      public function diagonal($m = null, $n = null, $c = 1)
 365      {
 366          $R = new self($m, $n);
 367          for ($i = 0; $i < $m; ++$i) {
 368              $R->set($i, $i, $c);
 369          }
 370  
 371          return $R;
 372      }
 373  
 374      /**
 375       * getMatrixByRow.
 376       *
 377       *    Get a submatrix by row index/range
 378       *
 379       * @param int $i0 Initial row index
 380       * @param int $iF Final row index
 381       *
 382       * @return Matrix Submatrix
 383       */
 384      public function getMatrixByRow($i0 = null, $iF = null)
 385      {
 386          if (is_int($i0)) {
 387              if (is_int($iF)) {
 388                  return $this->getMatrix($i0, 0, $iF + 1, $this->n);
 389              }
 390  
 391              return $this->getMatrix($i0, 0, $i0 + 1, $this->n);
 392          }
 393  
 394          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 395      }
 396  
 397      /**
 398       * getMatrixByCol.
 399       *
 400       *    Get a submatrix by column index/range
 401       *
 402       * @param int $j0 Initial column index
 403       * @param int $jF Final column index
 404       *
 405       * @return Matrix Submatrix
 406       */
 407      public function getMatrixByCol($j0 = null, $jF = null)
 408      {
 409          if (is_int($j0)) {
 410              if (is_int($jF)) {
 411                  return $this->getMatrix(0, $j0, $this->m, $jF + 1);
 412              }
 413  
 414              return $this->getMatrix(0, $j0, $this->m, $j0 + 1);
 415          }
 416  
 417          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 418      }
 419  
 420      /**
 421       * transpose.
 422       *
 423       *    Tranpose matrix
 424       *
 425       * @return Matrix Transposed matrix
 426       */
 427      public function transpose()
 428      {
 429          $R = new self($this->n, $this->m);
 430          for ($i = 0; $i < $this->m; ++$i) {
 431              for ($j = 0; $j < $this->n; ++$j) {
 432                  $R->set($j, $i, $this->A[$i][$j]);
 433              }
 434          }
 435  
 436          return $R;
 437      }
 438  
 439      //    function transpose()
 440  
 441      /**
 442       * trace.
 443       *
 444       *    Sum of diagonal elements
 445       *
 446       * @return float Sum of diagonal elements
 447       */
 448      public function trace()
 449      {
 450          $s = 0;
 451          $n = min($this->m, $this->n);
 452          for ($i = 0; $i < $n; ++$i) {
 453              $s += $this->A[$i][$i];
 454          }
 455  
 456          return $s;
 457      }
 458  
 459      /**
 460       * uminus.
 461       *
 462       *    Unary minus matrix -A
 463       *
 464       * @return Matrix Unary minus matrix
 465       */
 466      public function uminus()
 467      {
 468      }
 469  
 470      /**
 471       * plus.
 472       *
 473       *    A + B
 474       *
 475       * @return Matrix Sum
 476       */
 477      public function plus(...$args)
 478      {
 479          if (count($args) > 0) {
 480              $match = implode(',', array_map('gettype', $args));
 481  
 482              switch ($match) {
 483                  case 'object':
 484                      if ($args[0] instanceof self) {
 485                          $M = $args[0];
 486                      } else {
 487                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 488                      }
 489  
 490                      break;
 491                  case 'array':
 492                      $M = new self($args[0]);
 493  
 494                      break;
 495                  default:
 496                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 497  
 498                      break;
 499              }
 500              $this->checkMatrixDimensions($M);
 501              for ($i = 0; $i < $this->m; ++$i) {
 502                  for ($j = 0; $j < $this->n; ++$j) {
 503                      $M->set($i, $j, $M->get($i, $j) + $this->A[$i][$j]);
 504                  }
 505              }
 506  
 507              return $M;
 508          }
 509  
 510          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 511      }
 512  
 513      /**
 514       * plusEquals.
 515       *
 516       *    A = A + B
 517       *
 518       * @return $this
 519       */
 520      public function plusEquals(...$args)
 521      {
 522          if (count($args) > 0) {
 523              $match = implode(',', array_map('gettype', $args));
 524  
 525              switch ($match) {
 526                  case 'object':
 527                      if ($args[0] instanceof self) {
 528                          $M = $args[0];
 529                      } else {
 530                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 531                      }
 532  
 533                      break;
 534                  case 'array':
 535                      $M = new self($args[0]);
 536  
 537                      break;
 538                  default:
 539                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 540  
 541                      break;
 542              }
 543              $this->checkMatrixDimensions($M);
 544              for ($i = 0; $i < $this->m; ++$i) {
 545                  for ($j = 0; $j < $this->n; ++$j) {
 546                      $validValues = true;
 547                      $value = $M->get($i, $j);
 548                      if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
 549                          $this->A[$i][$j] = trim($this->A[$i][$j], '"');
 550                          $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
 551                      }
 552                      if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
 553                          $value = trim($value, '"');
 554                          $validValues &= StringHelper::convertToNumberIfFraction($value);
 555                      }
 556                      if ($validValues) {
 557                          $this->A[$i][$j] += $value;
 558                      } else {
 559                          $this->A[$i][$j] = Functions::NAN();
 560                      }
 561                  }
 562              }
 563  
 564              return $this;
 565          }
 566  
 567          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 568      }
 569  
 570      /**
 571       * minus.
 572       *
 573       *    A - B
 574       *
 575       * @return Matrix Sum
 576       */
 577      public function minus(...$args)
 578      {
 579          if (count($args) > 0) {
 580              $match = implode(',', array_map('gettype', $args));
 581  
 582              switch ($match) {
 583                  case 'object':
 584                      if ($args[0] instanceof self) {
 585                          $M = $args[0];
 586                      } else {
 587                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 588                      }
 589  
 590                      break;
 591                  case 'array':
 592                      $M = new self($args[0]);
 593  
 594                      break;
 595                  default:
 596                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 597  
 598                      break;
 599              }
 600              $this->checkMatrixDimensions($M);
 601              for ($i = 0; $i < $this->m; ++$i) {
 602                  for ($j = 0; $j < $this->n; ++$j) {
 603                      $M->set($i, $j, $M->get($i, $j) - $this->A[$i][$j]);
 604                  }
 605              }
 606  
 607              return $M;
 608          }
 609  
 610          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 611      }
 612  
 613      /**
 614       * minusEquals.
 615       *
 616       *    A = A - B
 617       *
 618       * @return $this
 619       */
 620      public function minusEquals(...$args)
 621      {
 622          if (count($args) > 0) {
 623              $match = implode(',', array_map('gettype', $args));
 624  
 625              switch ($match) {
 626                  case 'object':
 627                      if ($args[0] instanceof self) {
 628                          $M = $args[0];
 629                      } else {
 630                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 631                      }
 632  
 633                      break;
 634                  case 'array':
 635                      $M = new self($args[0]);
 636  
 637                      break;
 638                  default:
 639                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 640  
 641                      break;
 642              }
 643              $this->checkMatrixDimensions($M);
 644              for ($i = 0; $i < $this->m; ++$i) {
 645                  for ($j = 0; $j < $this->n; ++$j) {
 646                      $validValues = true;
 647                      $value = $M->get($i, $j);
 648                      if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
 649                          $this->A[$i][$j] = trim($this->A[$i][$j], '"');
 650                          $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
 651                      }
 652                      if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
 653                          $value = trim($value, '"');
 654                          $validValues &= StringHelper::convertToNumberIfFraction($value);
 655                      }
 656                      if ($validValues) {
 657                          $this->A[$i][$j] -= $value;
 658                      } else {
 659                          $this->A[$i][$j] = Functions::NAN();
 660                      }
 661                  }
 662              }
 663  
 664              return $this;
 665          }
 666  
 667          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 668      }
 669  
 670      /**
 671       * arrayTimes.
 672       *
 673       *    Element-by-element multiplication
 674       *    Cij = Aij * Bij
 675       *
 676       * @return Matrix Matrix Cij
 677       */
 678      public function arrayTimes(...$args)
 679      {
 680          if (count($args) > 0) {
 681              $match = implode(',', array_map('gettype', $args));
 682  
 683              switch ($match) {
 684                  case 'object':
 685                      if ($args[0] instanceof self) {
 686                          $M = $args[0];
 687                      } else {
 688                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 689                      }
 690  
 691                      break;
 692                  case 'array':
 693                      $M = new self($args[0]);
 694  
 695                      break;
 696                  default:
 697                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 698  
 699                      break;
 700              }
 701              $this->checkMatrixDimensions($M);
 702              for ($i = 0; $i < $this->m; ++$i) {
 703                  for ($j = 0; $j < $this->n; ++$j) {
 704                      $M->set($i, $j, $M->get($i, $j) * $this->A[$i][$j]);
 705                  }
 706              }
 707  
 708              return $M;
 709          }
 710  
 711          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 712      }
 713  
 714      /**
 715       * arrayTimesEquals.
 716       *
 717       *    Element-by-element multiplication
 718       *    Aij = Aij * Bij
 719       *
 720       * @return $this
 721       */
 722      public function arrayTimesEquals(...$args)
 723      {
 724          if (count($args) > 0) {
 725              $match = implode(',', array_map('gettype', $args));
 726  
 727              switch ($match) {
 728                  case 'object':
 729                      if ($args[0] instanceof self) {
 730                          $M = $args[0];
 731                      } else {
 732                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 733                      }
 734  
 735                      break;
 736                  case 'array':
 737                      $M = new self($args[0]);
 738  
 739                      break;
 740                  default:
 741                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 742  
 743                      break;
 744              }
 745              $this->checkMatrixDimensions($M);
 746              for ($i = 0; $i < $this->m; ++$i) {
 747                  for ($j = 0; $j < $this->n; ++$j) {
 748                      $validValues = true;
 749                      $value = $M->get($i, $j);
 750                      if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
 751                          $this->A[$i][$j] = trim($this->A[$i][$j], '"');
 752                          $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
 753                      }
 754                      if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
 755                          $value = trim($value, '"');
 756                          $validValues &= StringHelper::convertToNumberIfFraction($value);
 757                      }
 758                      if ($validValues) {
 759                          $this->A[$i][$j] *= $value;
 760                      } else {
 761                          $this->A[$i][$j] = Functions::NAN();
 762                      }
 763                  }
 764              }
 765  
 766              return $this;
 767          }
 768  
 769          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 770      }
 771  
 772      /**
 773       * arrayRightDivide.
 774       *
 775       *    Element-by-element right division
 776       *    A / B
 777       *
 778       * @return Matrix Division result
 779       */
 780      public function arrayRightDivide(...$args)
 781      {
 782          if (count($args) > 0) {
 783              $match = implode(',', array_map('gettype', $args));
 784  
 785              switch ($match) {
 786                  case 'object':
 787                      if ($args[0] instanceof self) {
 788                          $M = $args[0];
 789                      } else {
 790                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 791                      }
 792  
 793                      break;
 794                  case 'array':
 795                      $M = new self($args[0]);
 796  
 797                      break;
 798                  default:
 799                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 800  
 801                      break;
 802              }
 803              $this->checkMatrixDimensions($M);
 804              for ($i = 0; $i < $this->m; ++$i) {
 805                  for ($j = 0; $j < $this->n; ++$j) {
 806                      $validValues = true;
 807                      $value = $M->get($i, $j);
 808                      if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
 809                          $this->A[$i][$j] = trim($this->A[$i][$j], '"');
 810                          $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
 811                      }
 812                      if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
 813                          $value = trim($value, '"');
 814                          $validValues &= StringHelper::convertToNumberIfFraction($value);
 815                      }
 816                      if ($validValues) {
 817                          if ($value == 0) {
 818                              //    Trap for Divide by Zero error
 819                              $M->set($i, $j, '#DIV/0!');
 820                          } else {
 821                              $M->set($i, $j, $this->A[$i][$j] / $value);
 822                          }
 823                      } else {
 824                          $M->set($i, $j, Functions::NAN());
 825                      }
 826                  }
 827              }
 828  
 829              return $M;
 830          }
 831  
 832          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 833      }
 834  
 835      /**
 836       * arrayRightDivideEquals.
 837       *
 838       *    Element-by-element right division
 839       *    Aij = Aij / Bij
 840       *
 841       * @return Matrix Matrix Aij
 842       */
 843      public function arrayRightDivideEquals(...$args)
 844      {
 845          if (count($args) > 0) {
 846              $match = implode(',', array_map('gettype', $args));
 847  
 848              switch ($match) {
 849                  case 'object':
 850                      if ($args[0] instanceof self) {
 851                          $M = $args[0];
 852                      } else {
 853                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 854                      }
 855  
 856                      break;
 857                  case 'array':
 858                      $M = new self($args[0]);
 859  
 860                      break;
 861                  default:
 862                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 863  
 864                      break;
 865              }
 866              $this->checkMatrixDimensions($M);
 867              for ($i = 0; $i < $this->m; ++$i) {
 868                  for ($j = 0; $j < $this->n; ++$j) {
 869                      $this->A[$i][$j] = $this->A[$i][$j] / $M->get($i, $j);
 870                  }
 871              }
 872  
 873              return $M;
 874          }
 875  
 876          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 877      }
 878  
 879      /**
 880       * arrayLeftDivide.
 881       *
 882       *    Element-by-element Left division
 883       *    A / B
 884       *
 885       * @return Matrix Division result
 886       */
 887      public function arrayLeftDivide(...$args)
 888      {
 889          if (count($args) > 0) {
 890              $match = implode(',', array_map('gettype', $args));
 891  
 892              switch ($match) {
 893                  case 'object':
 894                      if ($args[0] instanceof self) {
 895                          $M = $args[0];
 896                      } else {
 897                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 898                      }
 899  
 900                      break;
 901                  case 'array':
 902                      $M = new self($args[0]);
 903  
 904                      break;
 905                  default:
 906                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 907  
 908                      break;
 909              }
 910              $this->checkMatrixDimensions($M);
 911              for ($i = 0; $i < $this->m; ++$i) {
 912                  for ($j = 0; $j < $this->n; ++$j) {
 913                      $M->set($i, $j, $M->get($i, $j) / $this->A[$i][$j]);
 914                  }
 915              }
 916  
 917              return $M;
 918          }
 919  
 920          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 921      }
 922  
 923      /**
 924       * arrayLeftDivideEquals.
 925       *
 926       *    Element-by-element Left division
 927       *    Aij = Aij / Bij
 928       *
 929       * @return Matrix Matrix Aij
 930       */
 931      public function arrayLeftDivideEquals(...$args)
 932      {
 933          if (count($args) > 0) {
 934              $match = implode(',', array_map('gettype', $args));
 935  
 936              switch ($match) {
 937                  case 'object':
 938                      if ($args[0] instanceof self) {
 939                          $M = $args[0];
 940                      } else {
 941                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 942                      }
 943  
 944                      break;
 945                  case 'array':
 946                      $M = new self($args[0]);
 947  
 948                      break;
 949                  default:
 950                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 951  
 952                      break;
 953              }
 954              $this->checkMatrixDimensions($M);
 955              for ($i = 0; $i < $this->m; ++$i) {
 956                  for ($j = 0; $j < $this->n; ++$j) {
 957                      $this->A[$i][$j] = $M->get($i, $j) / $this->A[$i][$j];
 958                  }
 959              }
 960  
 961              return $M;
 962          }
 963  
 964          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
 965      }
 966  
 967      /**
 968       * times.
 969       *
 970       *    Matrix multiplication
 971       *
 972       * @return Matrix Product
 973       */
 974      public function times(...$args)
 975      {
 976          if (count($args) > 0) {
 977              $match = implode(',', array_map('gettype', $args));
 978  
 979              switch ($match) {
 980                  case 'object':
 981                      if ($args[0] instanceof self) {
 982                          $B = $args[0];
 983                      } else {
 984                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
 985                      }
 986                      if ($this->n == $B->m) {
 987                          $C = new self($this->m, $B->n);
 988                          for ($j = 0; $j < $B->n; ++$j) {
 989                              $Bcolj = [];
 990                              for ($k = 0; $k < $this->n; ++$k) {
 991                                  $Bcolj[$k] = $B->A[$k][$j];
 992                              }
 993                              for ($i = 0; $i < $this->m; ++$i) {
 994                                  $Arowi = $this->A[$i];
 995                                  $s = 0;
 996                                  for ($k = 0; $k < $this->n; ++$k) {
 997                                      $s += $Arowi[$k] * $Bcolj[$k];
 998                                  }
 999                                  $C->A[$i][$j] = $s;
1000                              }
1001                          }
1002  
1003                          return $C;
1004                      }
1005  
1006                      throw new CalculationException(self::MATRIX_DIMENSION_EXCEPTION);
1007                  case 'array':
1008                      $B = new self($args[0]);
1009                      if ($this->n == $B->m) {
1010                          $C = new self($this->m, $B->n);
1011                          for ($i = 0; $i < $C->m; ++$i) {
1012                              for ($j = 0; $j < $C->n; ++$j) {
1013                                  $s = '0';
1014                                  for ($k = 0; $k < $C->n; ++$k) {
1015                                      $s += $this->A[$i][$k] * $B->A[$k][$j];
1016                                  }
1017                                  $C->A[$i][$j] = $s;
1018                              }
1019                          }
1020  
1021                          return $C;
1022                      }
1023  
1024                      throw new CalculationException(self::MATRIX_DIMENSION_EXCEPTION);
1025                  case 'integer':
1026                      $C = new self($this->A);
1027                      for ($i = 0; $i < $C->m; ++$i) {
1028                          for ($j = 0; $j < $C->n; ++$j) {
1029                              $C->A[$i][$j] *= $args[0];
1030                          }
1031                      }
1032  
1033                      return $C;
1034                  case 'double':
1035                      $C = new self($this->m, $this->n);
1036                      for ($i = 0; $i < $C->m; ++$i) {
1037                          for ($j = 0; $j < $C->n; ++$j) {
1038                              $C->A[$i][$j] = $args[0] * $this->A[$i][$j];
1039                          }
1040                      }
1041  
1042                      return $C;
1043                  case 'float':
1044                      $C = new self($this->A);
1045                      for ($i = 0; $i < $C->m; ++$i) {
1046                          for ($j = 0; $j < $C->n; ++$j) {
1047                              $C->A[$i][$j] *= $args[0];
1048                          }
1049                      }
1050  
1051                      return $C;
1052                  default:
1053                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
1054              }
1055          } else {
1056              throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
1057          }
1058      }
1059  
1060      /**
1061       * power.
1062       *
1063       *    A = A ^ B
1064       *
1065       * @return $this
1066       */
1067      public function power(...$args)
1068      {
1069          if (count($args) > 0) {
1070              $match = implode(',', array_map('gettype', $args));
1071  
1072              switch ($match) {
1073                  case 'object':
1074                      if ($args[0] instanceof self) {
1075                          $M = $args[0];
1076                      } else {
1077                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
1078                      }
1079  
1080                      break;
1081                  case 'array':
1082                      $M = new self($args[0]);
1083  
1084                      break;
1085                  default:
1086                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
1087  
1088                      break;
1089              }
1090              $this->checkMatrixDimensions($M);
1091              for ($i = 0; $i < $this->m; ++$i) {
1092                  for ($j = 0; $j < $this->n; ++$j) {
1093                      $validValues = true;
1094                      $value = $M->get($i, $j);
1095                      if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
1096                          $this->A[$i][$j] = trim($this->A[$i][$j], '"');
1097                          $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
1098                      }
1099                      if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
1100                          $value = trim($value, '"');
1101                          $validValues &= StringHelper::convertToNumberIfFraction($value);
1102                      }
1103                      if ($validValues) {
1104                          $this->A[$i][$j] = $this->A[$i][$j] ** $value;
1105                      } else {
1106                          $this->A[$i][$j] = Functions::NAN();
1107                      }
1108                  }
1109              }
1110  
1111              return $this;
1112          }
1113  
1114          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
1115      }
1116  
1117      /**
1118       * concat.
1119       *
1120       *    A = A & B
1121       *
1122       * @return $this
1123       */
1124      public function concat(...$args)
1125      {
1126          if (count($args) > 0) {
1127              $match = implode(',', array_map('gettype', $args));
1128  
1129              switch ($match) {
1130                  case 'object':
1131                      if ($args[0] instanceof self) {
1132                          $M = $args[0];
1133                      } else {
1134                          throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
1135                      }
1136  
1137                      break;
1138                  case 'array':
1139                      $M = new self($args[0]);
1140  
1141                      break;
1142                  default:
1143                      throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
1144  
1145                      break;
1146              }
1147              $this->checkMatrixDimensions($M);
1148              for ($i = 0; $i < $this->m; ++$i) {
1149                  for ($j = 0; $j < $this->n; ++$j) {
1150                      $this->A[$i][$j] = trim($this->A[$i][$j], '"') . trim($M->get($i, $j), '"');
1151                  }
1152              }
1153  
1154              return $this;
1155          }
1156  
1157          throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
1158      }
1159  
1160      /**
1161       * Solve A*X = B.
1162       *
1163       * @param Matrix $B Right hand side
1164       *
1165       * @return Matrix ... Solution if A is square, least squares solution otherwise
1166       */
1167      public function solve($B)
1168      {
1169          if ($this->m == $this->n) {
1170              $LU = new LUDecomposition($this);
1171  
1172              return $LU->solve($B);
1173          }
1174          $QR = new QRDecomposition($this);
1175  
1176          return $QR->solve($B);
1177      }
1178  
1179      /**
1180       * Matrix inverse or pseudoinverse.
1181       *
1182       * @return Matrix ... Inverse(A) if A is square, pseudoinverse otherwise.
1183       */
1184      public function inverse()
1185      {
1186          return $this->solve($this->identity($this->m, $this->m));
1187      }
1188  
1189      /**
1190       * det.
1191       *
1192       *    Calculate determinant
1193       *
1194       * @return float Determinant
1195       */
1196      public function det()
1197      {
1198          $L = new LUDecomposition($this);
1199  
1200          return $L->det();
1201      }
1202  }