Matrix.php 37 KB

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