| 1 |
<?php |
| 2 |
/** |
| 3 |
* Copyright (C) 2023 Graham Breach |
| 4 |
* |
| 5 |
* This program is free software: you can redistribute it and/or modify |
| 6 |
* it under the terms of the GNU Lesser General Public License as published by |
| 7 |
* the Free Software Foundation, either version 3 of the License, or |
| 8 |
* (at your option) any later version. |
| 9 |
* |
| 10 |
* This program is distributed in the hope that it will be useful, |
| 11 |
* but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 12 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 13 |
* GNU Lesser General Public License for more details. |
| 14 |
* |
| 15 |
* You should have received a copy of the GNU Lesser General Public License |
| 16 |
* along with this program. If not, see <http://www.gnu.org/licenses/>. |
| 17 |
*/ |
| 18 |
/** |
| 19 |
* For more information, please contact <graham@goat1000.com> |
| 20 |
*/ |
| 21 |
|
| 22 |
namespace Goat1000\SVGGraph; |
| 23 |
|
| 24 |
/** |
| 25 |
* Class for matrix maths |
| 26 |
*/ |
| 27 |
class Matrix { |
| 28 |
private $rows = 0; |
| 29 |
private $cols = 0; |
| 30 |
private $data = []; |
| 31 |
const MAX_SIZE = 100000; |
| 32 |
|
| 33 |
public function __construct($r, $c) |
| 34 |
{ |
| 35 |
if(!is_int($c) || !is_int($r) || $c < 1 || $r < 1) |
| 36 |
throw new \InvalidArgumentException("{$r}\u{00d7}{$c} matrix size invalid"); |
| 37 |
if($c * $r > Matrix::MAX_SIZE) |
| 38 |
throw new \InvalidArgumentException("{$r}\u{00d7}{$c} matrix too large"); |
| 39 |
|
| 40 |
$this->rows = $r; |
| 41 |
$this->cols = $c; |
| 42 |
$this->data = array_fill(0, $c * $r, 0); |
| 43 |
} |
| 44 |
|
| 45 |
/** |
| 46 |
* Returns the number of rows and cols in an array |
| 47 |
*/ |
| 48 |
public function dimensions() |
| 49 |
{ |
| 50 |
return [$this->rows, $this->cols]; |
| 51 |
} |
| 52 |
|
| 53 |
/** |
| 54 |
* Returns the size as a string |
| 55 |
*/ |
| 56 |
public function size() |
| 57 |
{ |
| 58 |
return $this->rows . "\u{00d7}" . $this->cols; |
| 59 |
} |
| 60 |
|
| 61 |
/** |
| 62 |
* Element access |
| 63 |
*/ |
| 64 |
public function &__invoke($r, $c, $v = null) |
| 65 |
{ |
| 66 |
if($c < 0 || $r < 0 || $c >= $this->cols || $r >= $this->rows) |
| 67 |
throw new \InvalidArgumentException("({$r},{$c}) out of range of " . |
| 68 |
$this->size() . " matrix"); |
| 69 |
|
| 70 |
$item = $this->cols * $r + $c; |
| 71 |
if($v !== null) { |
| 72 |
if(!is_numeric($v)) |
| 73 |
throw new \InvalidArgumentException("Matrix values must be numeric"); |
| 74 |
if(is_string($v)) |
| 75 |
$this->data[$item] = $v; |
| 76 |
elseif(is_int($v)) |
| 77 |
$this->data[$item] = (string)$v; |
| 78 |
else |
| 79 |
$this->data[$item] = sprintf("%.40F", $v); |
| 80 |
} |
| 81 |
return $this->data[$item]; |
| 82 |
} |
| 83 |
|
| 84 |
/** |
| 85 |
* Fill the array with values (row-major order) |
| 86 |
*/ |
| 87 |
public function load(array $values) |
| 88 |
{ |
| 89 |
for($i = 0; $i < $this->rows * $this->cols; ++$i) { |
| 90 |
$value = $values[$i]; |
| 91 |
if(isset($value) && is_numeric($value)) { |
| 92 |
$this->data[$i] = is_int($value) ? (string)$value : |
| 93 |
sprintf("%.40F", $value); |
| 94 |
} |
| 95 |
} |
| 96 |
} |
| 97 |
|
| 98 |
/** |
| 99 |
* Sets the identity matrix |
| 100 |
*/ |
| 101 |
public function identity() |
| 102 |
{ |
| 103 |
if($this->rows != $this->cols) |
| 104 |
throw new \Exception($this->size() . " not a square matrix"); |
| 105 |
$this->data = array_fill(0, $this->rows * $this->cols, 0); |
| 106 |
for($i = 0; $i < $this->rows; ++$i) |
| 107 |
$this->data[($i * $this->rows) + $i] = 1; |
| 108 |
} |
| 109 |
|
| 110 |
/** |
| 111 |
* Return the transpose of the matrix |
| 112 |
*/ |
| 113 |
public function transpose() |
| 114 |
{ |
| 115 |
$tc = $this->rows; |
| 116 |
$tr = $this->cols; |
| 117 |
$result = new Matrix($tr, $tc); |
| 118 |
|
| 119 |
// row or column vector = same data when transposed |
| 120 |
if($tc == 1 || $tr == 1) { |
| 121 |
$result->data = $this->data; |
| 122 |
return $result; |
| 123 |
} |
| 124 |
|
| 125 |
for($c = 0; $c < $tc; ++$c) { |
| 126 |
for($r = 0; $r < $tr; ++$r) { |
| 127 |
$result($r, $c, $this($c, $r)); |
| 128 |
} |
| 129 |
} |
| 130 |
return $result; |
| 131 |
} |
| 132 |
|
| 133 |
/** |
| 134 |
* Add two matrices |
| 135 |
*/ |
| 136 |
public function add(Matrix $m, $bcscale = 50) |
| 137 |
{ |
| 138 |
if($m->rows != $this->rows || $m->cols != $this->cols) |
| 139 |
throw new \InvalidArgumentException("Cannot add " . $this->size() . |
| 140 |
" and " . $m->size() . " matrices."); |
| 141 |
|
| 142 |
$old_scale = bcscale($bcscale); |
| 143 |
$result = new Matrix($this->cols, $this->rows); |
| 144 |
foreach($this->data as $k => $value) |
| 145 |
$result->data[$k] = bcadd($value, $m->data[$k]); |
| 146 |
|
| 147 |
bcscale($old_scale); |
| 148 |
return $result; |
| 149 |
} |
| 150 |
|
| 151 |
/** |
| 152 |
* Subtract a matrix |
| 153 |
*/ |
| 154 |
public function subtract(Matrix $m, $bcscale = 50) |
| 155 |
{ |
| 156 |
if($m->rows != $this->rows || $m->cols != $this->cols) |
| 157 |
throw new \InvalidArgumentException("Cannot subtract " . $m->size() . |
| 158 |
" matrix from " . $this->size() . " matrix."); |
| 159 |
|
| 160 |
$old_scale = bcscale($bcscale); |
| 161 |
$result = new Matrix($this->cols, $this->rows); |
| 162 |
foreach($this->data as $k => $value) |
| 163 |
$result->data[$k] = bcsub($value, $m->data[$k]); |
| 164 |
|
| 165 |
bcscale($old_scale); |
| 166 |
return $result; |
| 167 |
} |
| 168 |
|
| 169 |
/** |
| 170 |
* Multiplication by a matrix |
| 171 |
*/ |
| 172 |
public function multiply(Matrix $m1, $bcscale = 50) |
| 173 |
{ |
| 174 |
if($this->cols != $m1->rows) |
| 175 |
throw new \InvalidArgumentException("Cannot multiply " . $this->size() . |
| 176 |
" matrix by " . $m1->size() . " matrix."); |
| 177 |
$m = $this->rows; |
| 178 |
$n = $this->cols; |
| 179 |
$p = $m1->cols; |
| 180 |
|
| 181 |
$old_scale = bcscale($bcscale); |
| 182 |
$result = new Matrix($m, $p); |
| 183 |
for($i = 0; $i < $m; ++$i) { |
| 184 |
for($j = 0; $j < $p; ++$j) { |
| 185 |
$value = "0"; |
| 186 |
for($k = 0; $k < $n; ++$k) { |
| 187 |
$value = bcadd($value, bcmul($this($i, $k), $m1($k, $j))); |
| 188 |
} |
| 189 |
$result($i, $j, $value); |
| 190 |
} |
| 191 |
} |
| 192 |
|
| 193 |
bcscale($old_scale); |
| 194 |
return $result; |
| 195 |
} |
| 196 |
|
| 197 |
/** |
| 198 |
* Gaussian elimination |
| 199 |
*/ |
| 200 |
public function gaussian($bcscale = 50) |
| 201 |
{ |
| 202 |
$old_scale = bcscale($bcscale); |
| 203 |
|
| 204 |
$argmax = function($a, $b, $m, $col) { |
| 205 |
$max = 0; |
| 206 |
$max_i = 0; |
| 207 |
for($i = $a; $i < $b; ++$i) { |
| 208 |
$value = $m($a, $col); |
| 209 |
if(bccomp($value, "0") == -1) |
| 210 |
$value = bcmul($value, "-1"); |
| 211 |
if(bccomp($value, $max) == 1) { |
| 212 |
$max_i = $i; |
| 213 |
$max = $value; |
| 214 |
} |
| 215 |
return $max_i; |
| 216 |
} |
| 217 |
}; |
| 218 |
|
| 219 |
$m = $this->rows; |
| 220 |
$n = $this->cols; |
| 221 |
$h = $k = 0; |
| 222 |
while($h < $m && $k < $n) { |
| 223 |
$i_max = $argmax($h, $m, $this, $k); |
| 224 |
if($this($i_max, $k) == 0) { |
| 225 |
++$k; |
| 226 |
} else { |
| 227 |
$this->rowSwap($h, $i_max); |
| 228 |
for($i = $h + 1; $i < $m; ++$i) { |
| 229 |
$f = bcdiv($this($i, $k), $this($h, $k)); |
| 230 |
$this($i, $k, 0); |
| 231 |
for($j = $k + 1; $j < $n; ++$j) { |
| 232 |
$val = bcsub($this($i, $j), bcmul($this($h, $j), $f)); |
| 233 |
$this($i, $j, $val); |
| 234 |
} |
| 235 |
} |
| 236 |
++$h; |
| 237 |
++$k; |
| 238 |
} |
| 239 |
} |
| 240 |
|
| 241 |
bcscale($old_scale); |
| 242 |
} |
| 243 |
|
| 244 |
/** |
| 245 |
* Use Gaussian elimination to solve the equations with given RHS |
| 246 |
*/ |
| 247 |
public function gaussian_solve(Matrix $rhs, $bcscale = 50) |
| 248 |
{ |
| 249 |
$a = $this->augment($rhs); |
| 250 |
$a->gaussian($bcscale); |
| 251 |
return $this->solve($a, $bcscale); |
| 252 |
} |
| 253 |
|
| 254 |
/** |
| 255 |
* Creates a new matrix with $this on left and $rhs on right |
| 256 |
*/ |
| 257 |
public function augment(Matrix $rhs) |
| 258 |
{ |
| 259 |
$m = $this->rows; |
| 260 |
$n = $this->cols + $rhs->cols; |
| 261 |
$aug = new Matrix($m, $n); |
| 262 |
|
| 263 |
$c = 0; |
| 264 |
for($i = 0; $i < $m; ++$i) { |
| 265 |
for($j = 0; $j < $this->cols; ++$j) |
| 266 |
$aug->data[$c++] = $this($i, $j); |
| 267 |
for($j = 0; $j < $rhs->cols; ++$j) |
| 268 |
$aug->data[$c++] = $rhs($i, $j); |
| 269 |
} |
| 270 |
return $aug; |
| 271 |
} |
| 272 |
|
| 273 |
/** |
| 274 |
* Solves simultaneous equations using Gaussian elimination |
| 275 |
*/ |
| 276 |
public function solve(Matrix $a, $bcscale) |
| 277 |
{ |
| 278 |
$result = new Matrix(1, $a->rows); |
| 279 |
$old_scale = bcscale($bcscale); |
| 280 |
|
| 281 |
// back substitution |
| 282 |
$m = $a->rows; |
| 283 |
$n = $a->cols; |
| 284 |
for($i = $m - 1; $i >= 0; --$i) { |
| 285 |
for($j = $n - 2; $j > $i; --$j) { |
| 286 |
$value = bcsub($a($i, $n - 1), |
| 287 |
bcmul($a($i, $j), $result(0, $j))); |
| 288 |
$a($i, $n - 1, $value); |
| 289 |
$a($i, $j, 0); |
| 290 |
} |
| 291 |
$d = $a($i, $i); |
| 292 |
if($d == 0) |
| 293 |
return null; |
| 294 |
$value = bcdiv($a($i, $n - 1), $a($i, $i)); |
| 295 |
$result(0, $i, $value); |
| 296 |
} |
| 297 |
|
| 298 |
bcscale($old_scale); |
| 299 |
return $result; |
| 300 |
} |
| 301 |
|
| 302 |
/** |
| 303 |
* Swaps two rows |
| 304 |
*/ |
| 305 |
public function rowSwap($r1, $r2) |
| 306 |
{ |
| 307 |
for($i = 0; $i < $this->cols; ++$i) { |
| 308 |
$c = $this($r1, $i); |
| 309 |
$this($r1, $i, $this($r2, $i)); |
| 310 |
$this($r2, $i, $c); |
| 311 |
} |
| 312 |
} |
| 313 |
|
| 314 |
/** |
| 315 |
* Output as string for debugging |
| 316 |
*/ |
| 317 |
public function __toString() |
| 318 |
{ |
| 319 |
$str = ''; |
| 320 |
$m = 0; |
| 321 |
foreach($this->data as $v) { |
| 322 |
$m1 = abs($v); |
| 323 |
if($m1 > $m) |
| 324 |
$m = $m1; |
| 325 |
} |
| 326 |
|
| 327 |
$digits = max(9,(int)log($m, 10) + 6); |
| 328 |
for($r = 0; $r < $this->rows; ++$r) { |
| 329 |
$str .= "\t"; |
| 330 |
$r_offset = $r * $this->cols; |
| 331 |
for($c = 0; $c < $this->cols; ++$c) { |
| 332 |
$str .= sprintf(" %{$digits}.4f", $this->data[$r_offset + $c]); |
| 333 |
} |
| 334 |
$str .= "\n"; |
| 335 |
} |
| 336 |
return $str; |
| 337 |
} |
| 338 |
|
| 339 |
/** |
| 340 |
* Returns the data array |
| 341 |
*/ |
| 342 |
public function asArray() |
| 343 |
{ |
| 344 |
return $this->data; |
| 345 |
} |
| 346 |
} |
| 347 |
|
