Passed
Push — main ( 71dcec...306c22 )
by Shubham
01:53
created

matrix::clip()   A

Complexity

Conditions 4
Paths 4

Size

Total Lines 14
Code Lines 10

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
cc 4
eloc 10
nc 4
nop 2
dl 0
loc 14
rs 9.9332
c 0
b 0
f 0
1
<?php
2
3
declare(strict_types=1);
4
5
namespace Np;
6
7
use Np\core\nd;
8
use Np\core\blas;
9
use Np\core\lapack;
10
use Np\reductions\ref;
11
use Np\reductions\rref;
12
use Np\decompositions\lu;
13
use Np\decompositions\svd;
14
use Np\decompositions\eigen;
15
use Np\decompositions\cholesky;
16
17
/**
18
 * Matrix
19
 * A fast lite memory efficient Scientific Computing for php
20
 * 
21
 * @package   NumPhp
22
 * @category  Scientific Computing
23
 * @author    ghost (Shubham Chaudhary)
24
 * @email     [email protected]
25
 * @copyright (c) 2020-2021, Shubham Chaudhary
26
 */
27
class matrix extends nd{
28
    
29
    use ops,linAlg;
30
    /**
31
     * create empty 2d matrix for given data type
32
     * @param int $row  num of rows 
33
     * @param int $col  num of cols
34
     * @param int $dtype matrix data type float|double
35
     * @return \Np\matrix
36
     */
37
    public static function factory(int $row, int $col, int $dtype = self::FLOAT): matrix {
38
        return new self($row, $col, $dtype);
39
    }
40
41
    /**
42
     * create 2d matrix using php array
43
     * @param array $data
44
     * @param int $dtype matrix data type float|double
45
     * @return \Np\matrix
46
     */
47
    public static function ar(array $data, int $dtype = self::FLOAT): matrix {
48
        if (is_array($data) && is_array($data[0])) {
49
            $ar = self::factory(count($data), count($data[0]), $dtype);
50
            $ar->setData($data);
51
            unset($data);
52
            return $ar;
53
        } else {
54
            self::_err('given array is not rank-2 or given is not an array');
55
        }
56
    }
57
58
    /**
59
     * create one like 2d matrix
60
     * @param int $row
61
     * @param int $col
62
     * @return \Np\matrix
63
     */
64
    public static function ones(int $row, int $col, int $dtype = self::FLOAT): matrix {
65
        $ar = self::factory($row, $col, $dtype);
66
        for ($i = 0; $i < $ar->ndim; ++$i) {
67
            $ar->data[$i] = 1;
68
        }
69
        return $ar;
70
    }
71
72
    /**
73
     * Create Matrix with random values
74
     * @param int $row
75
     * @param int $col
76
     * @param int $dtype  Float|Double
77
     * @return \Np\matrix
78
     */
79
    public static function randn(int $row, int $col, int $dtype = self::FLOAT): matrix {
80
        $ar = self::factory($row, $col, $dtype);
81
        $max = getrandmax();
82
        for ($i = 0; $i < $ar->ndim; ++$i) {
83
            $ar->data[$i] = rand() / $max;
84
        }
85
        return $ar;
86
    }
87
88
    /**
89
     * Return 2d matrix with uniform values
90
     * @param int $row
91
     * @param int $col
92
     * @param int $dtype
93
     * @return \Np\matrix
94
     */
95
    public static function uniform(int $row, int $col, int $dtype = self::FLOAT): matrix {
96
        $ar = self::factory($row, $col, $dtype);
97
        $max = getrandmax();
98
        for ($i = 0; $i < $ar->ndim; ++$i) {
99
            $ar->data[$i] = rand(-$max, $max) / $max;
100
        }
101
        return $ar;
102
    }
103
104
    /**
105
     * Return a zero matrix with the given dimensions.
106
     * @param int $row
107
     * @param int $col
108
     * @param int $dtype
109
     * @return \Np\matrix
110
     */
111
    public static function zeros(int $row, int $col, int $dtype = self::FLOAT): matrix {
112
        $ar = self::factory($row, $col, $dtype);
113
        for ($i = 0; $i < $ar->ndim; ++$i) {
114
            $ar->data[$i] = 0.0;
115
        }
116
        return $ar;
117
    }
118
119
    /**
120
     * create a null like 2d matrix
121
     * @param int $row
122
     * @param int $col
123
     * @return \Np\matrix
124
     */
125
    public static function null(int $row, int $col, int $dtype = self::FLOAT): matrix {
126
        $ar = self::factory($row, $col, $dtype);
127
        for ($i = 0; $i < $ar->ndim; ++$i) {
128
            $ar->data[$i] = null;
129
        }
130
        return $ar;
131
    }
132
133
    /**
134
     * create a 2d matrix with given scalar value
135
     * @param int $row
136
     * @param int $col
137
     * @param int|float|double $val
138
     * @return \Np\matrix
139
     */
140
    public static function full(int $row, int $col, $val, int $dtype = self::FLOAT): matrix {
141
        $ar = self::factory($row, $col, $dtype);
142
        for ($i = 0; $i < $ar->ndim; ++$i) {
143
            $ar->data[$i] = $val;
144
        }
145
        return $ar;
146
    }
147
148
    /**
149
     * create a diagonal 2d matrix with given 1d array;
150
     * @param array $elements
151
     * @return \Np\matrix
152
     */
153
    public static function diagonal(array $elements, int $dtype = self::FLOAT): matrix {
154
        $n = count($elements);
155
        $ar = self::factory($n, $n, $dtype);
156
        for ($i = 0; $i < $n; ++$i) {
157
            $ar->data[$i * $n + $i] = $elements[$i]; #for ($j = 0; $j < $n; ++$j) {$i === $j ? $elements[$i] : 0;#} 
158
        }
159
        return $ar;
160
    }
161
162
    /**
163
     * Generate a m x n matrix with elements from a Poisson distribution.
164
     * @param int $row
165
     * @param int $col
166
     * @param float $lambda
167
     * @param int $dtype 
168
     * @return \Np\matrix
169
     */
170
    public static function poisson(int $row, int $col, float $lambda = 1.0, int $dtype = self::FLOAT): matrix {
171
        $ar = self::factory($row, $col, $dtype);
172
        $max = getrandmax();
173
        $l = exp(-$lambda);
174
        for ($i = 0; $i < $row; ++$i) {
175
            for ($j = 0; $j < $col; ++$j) {
176
                $k = 0;
177
                $p = 1.0;
178
                while ($p > $l) {
179
                    ++$k;
180
                    $p = $p * rand() / $max;
181
                }
182
                $ar->data[$i * $col + $j] = $k - 1;
183
            }
184
        }
185
        return $ar;
186
    }
187
188
    /**
189
     * Return a standard normally distributed random matrix i.e values
190
     * between -1 and 1.
191
     * @param int $row
192
     * @param int $col
193
     * @param int $dtype Description
194
     * @return \Np\matrix
195
     */
196
    public static function gaussian(int $row, int $col, int $dtype = self::FLOAT): matrix {
197
        $max = getrandmax();
198
        $a = $extras = [];
199
200
        while (count($a) < $row) {
201
            $rowA = [];
202
203
            if (!empty($extras)) {
204
                $rowA[] = array_pop($extras);
205
            }
206
207
            while (count($rowA) < $col) {
208
                $r = sqrt(-2.0 * log(rand() / $max));
209
210
                $phi = rand() / $max * self::TWO_PI;
211
212
                $rowA[] = $r * sin($phi);
213
                $rowA[] = $r * cos($phi);
214
            }
215
216
            if (count($rowA) > $col) {
217
                $extras[] = array_pop($rowA);
218
            }
219
220
            $a[] = $rowA;
221
        }
222
223
        return self::ar($a, $dtype);
224
    }
225
226
    /**
227
     * create an identity matrix with the given dimensions.
228
     * @param int $n
229
     * @param int $dtype
230
     * @return matrix
231
     * @throws \InvalidArgumentException
232
     */
233
    public static function identity(int $n, int $dtype = self::FLOAT): matrix {
234
        if ($n < 1) {
235
            self::_dimensionaMisMatchErr('dimensionality must be greater than 0 on all axes.');
236
        }
237
238
        $ar = self::factory($n, $n, $dtype);
239
        for ($i = 0; $i < $n; ++$i) {
240
            for ($j = 0; $j < $n; ++$j) {
241
                $ar->data[$i * $n + $j] = $i === $j ? 1 : 0;
242
            }
243
        }
244
        return $ar;
245
    }
246
247
    /**
248
     * 2D convolution between a matrix ma and kernel kb, with a given stride.
249
     * @param \Np\matrix $m
250
     * @param int $stride
251
     * @return matrix
252
     */
253
    public function convolve(matrix $m, int $stride = 1): matrix {
254
        return convolve::conv2D($this, $m, $stride);
255
    }
256
257
    /**
258
     * Calculate the determinant of the matrix.
259
     * @return float
260
     */
261
    public function det(): float {
262
        if (!$this->isSquare()) {
263
            self::_err('determinant is undefined for a non square matrix');
264
        }
265
        $lu = $this->lu();
266
        $nSwaps = $lu->p()->diagonalAsVector()->subtract($lu->p()->diagonalAsVector()->sum())->col - 1;
267
        $detP = (-1) ** $nSwaps;
268
        $detL = $lu->l()->diagonalAsVector()->product();
269
        $detU = $lu->u()->diagonalAsVector()->product();
270
        unset($lu);
271
        return ($detP * $detL * $detU);
272
    }
273
274
    /**
275
     * Return the trace of the matrix i.e the sum of all diagonal elements of a square matrix.
276
     * @return float
277
     */
278
    public function trace(): float {
279
        if (!$this->isSquare()) {
280
            self::_err('Error::matrix is not a squared can not Trace!');
281
        }
282
        $trace = 0.0;
283
        for ($i = 0; $i < $this->row; ++$i) {
284
            for ($j = 0; $j < $this->col; ++$j) {
285
                if ($i == $j) {
286
                    $trace += $this->data[$i * $this->col + $i];
287
                }
288
            }
289
        }
290
        return $trace;
291
    }
292
293
    /**
294
     * dignoalInterChange
295
     */
296
    public function dignoalInterChange() {
297
        for ($i = 0; $i < $this->row; ++$i) {
298
            for ($j = 0; $j < $this->col; ++$j) {
299
                $tmp = $this->data[$i * $this->col - $j];
300
                $this->data[$i * $this->col - $j] = $tmp;
301
            }
302
        }
303
    }
304
    
305
    //---------------Arthmetic Opreations-----------------------------------
306
307
    /**
308
     * multiply this matrix with another matrix|scalar element-wise
309
     * Matrix Scalar\Matrix multiplication
310
     * @param int|float|matrix|vector $m
311
     * @return matrix|vector
312
     */
313
    public function multiply(int|float|matrix|vector $m): matrix|vector {
314
        if ($m instanceof self) {
315
            return $this->multiplyMatrix($m);
316
        } else if ($m instanceof vector) {
317
            return $this->multiplyVector($m);
318
        } else {
319
            return $this->scale($m);
320
        }
321
    }
322
323
    /**
324
     * 
325
     * @param \Np\vector $v
326
     * @return matrix
327
     */
328
    protected function multiplyVector(vector $v): matrix {
329
        if ($this->checkDimensions($v, $this) && $this->checkDtype($this, $v)) {
330
            $ar = matrix::factory($this->row, $this->col, $this->dtype);
331
            for ($i = 0; $i < $this->row; ++$i) {
332
                for ($j = 0; $j < $this->col; ++$j) {
333
                    $ar->data[$i * $this->col + $j] = $v->data[$j] * $this->data[$i * $this->col + $j];
334
                }
335
            }
336
            return $ar;
337
        }
338
    }
339
340
    /**
341
     * 
342
     * @param \Np\matrix $m
343
     * @return matrix
344
     */
345
    protected function multiplyMatrix(matrix $m): matrix {
346
        if ($this->checkDtype($this, $m) && $this->checkShape($this, $m)) {
347
            $ar = self::factory($this->row, $this->col, $this->dtype);
348
            for ($i = 0; $i < $this->row; ++$i) {
349
                for ($j = 0; $j < $this->col; ++$j) {
350
                    $ar->data[$i * $this->col + $j] = $this->data[$i * $this->col + $j] * $m->data[$i * $this->col + $j];
351
                }
352
            }
353
            return $ar;
354
        }
355
    }
356
357
    /**
358
     * Sum of Rows of matrix
359
     * @return vector
360
     */
361
    public function sumRows(): vector {
362
        $vr = vector::factory($this->row, $this->dtype);
363
        for ($i = 0; $i < $this->row; ++$i) {
364
            $sum = 0.0;
365
            for ($j = 0; $j < $this->col; ++$j) {
366
                $sum += $this->data[$i * $this->col + $j];
367
            }
368
            $vr->data[$i] = $sum;
369
        }
370
        return $vr;
371
    }
372
373
    /**
374
     * Sum of two matrix, vector or a scalar to current matrix
375
     * 
376
     * @param int|float|matrix|vector $m
377
     * @return matrix
378
     */
379
    public function sum(int|float|matrix|vector $m): matrix {
380
        if ($m instanceof self) {
381
            return $this->sumMatrix($m);
382
        } elseif ($m instanceof vector) {
383
            return $this->sumVector($m);
384
        } else {
385
            return $this->sumScalar($m);
386
        }
387
    }
388
389
    protected function sumScalar(int|float $s): matrix {
390
        $ar = self::factory($this->row, $this->col, $this->dtype);
391
        for ($i = 0; $i < $this->ndim; ++$i) {
392
            $ar->data[$i] = $this->data[$i] + $s;
393
        }
394
        return $ar;
395
    }
396
397
    protected function sumMatrix(matrix $m): matrix {
398
        if ($this->checkShape($this, $m) && $this->checkDtype($this, $m)) {
399
            $ar = self::factory($this->row, $this->col, $this->dtype);
400
            for ($i = 0; $i < $this->ndim; ++$i) {
401
                $ar->data[$i] = $this->data[$i] + $m->data[$i];
402
            }
403
            return $ar;
404
        }
405
    }
406
407
    protected function sumVector(vector $v): matrix {
408
        if ($this->checkDimensions($v, $this) && $this->checkDtype($this, $v)) {
409
            $ar = self::factory($this->row, $this->col, $this->dtype);
410
            for ($i = 0; $i < $this->row; ++$i) {
411
                for ($j = 0; $j < $this->col; ++$j) {
412
                    $ar->data[$i * $this->col + $j] = $v->data[$j] + $this->data[$i * $this->col + $j];
413
                }
414
            }
415
            return $ar;
416
        }
417
    }
418
419
    /**
420
     * subtract another matrix, vector or a scalar to this matrix
421
     * @param int|float|matrix|vector $d matrix|$scalar to subtract this matrix
422
     * @return \Np\matrix
423
     */
424
    public function subtract(int|float|matrix|vector $d): matrix {
425
        if ($d instanceof self) {
426
            return $this->subtractMatrix($d);
427
        } elseif ($d instanceof vector) {
428
            return $this->subtractVector($d);
429
        } else {
430
            return $this->subtractScalar($d);
431
        }
432
    }
433
434
    protected function subtractScalar(int|float $s): matrix {
435
        $ar = self::factory($this->row, $this->col, $this->dtype);
436
        for ($i = 0; $i < $this->ndim; ++$i) {
437
            $ar->data[$i] = $this->data[$i] - $s;
438
        }
439
        return $ar;
440
    }
441
442
    /**
443
     * 
444
     * @param matrix $m
445
     * @return matrix
446
     */
447
    protected function subtractMatrix(matrix $m): matrix {
448
        if ($this->checkShape($this, $m) && $this->checkDtype($this, $m)) {
449
            $ar = self::factory($this->row, $this->col, $this->dtype);
450
            for ($i = 0; $i < $this->ndim; ++$i) {
451
                $ar->data[$i] = $this->data[$i] - $m->data[$i];
452
            }
453
            return $ar;
454
        }
455
    }
456
457
    /**
458
     * 
459
     * @param vector $v
460
     * @return matrix
461
     */
462
    protected function subtractVector(vector $v): matrix {
463
        if ($this->checkDimensions($v, $this) && $this->checkDtype($this,$v)) {
464
            $ar = self::factory($this->row, $this->col, $this->dtype);
465
            for ($i = 0; $i < $this->row; ++$i) {
466
                for ($j = 0; $j < $this->col; ++$j) {
467
                    $ar->data[$i * $this->col + $j] = $this->data[$i * $this->col + $j] - $v->data[$j];
468
                }
469
            }
470
            return $ar;
471
        }
472
    }
473
474
    /**
475
     * 
476
     * @param vector $v
477
     * @return matrix
478
     */
479
    public function subtractColumnVector(vector $v): matrix {
480
        if ($this->checkDimensions($v, $this) && $this->checkDtype($this, $v)) {
481
            $ar = self::factory($this->row, $this->col, $this->dtype);
482
            for ($j = 0; $j < $this->col; ++$j) {
483
                for ($i = 0; $i < $this->row; ++$i) {
484
                    $ar->data[$i * $this->col + $j] = $this->data[$i * $this->col + $j] - $v->data[$i];
485
                }
486
            }
487
            return $ar;
488
        }
489
    }
490
491
    /**
492
     * Return the division of two elements, element-wise.
493
     * @param int|float|matrix $d
494
     * @return matrix
495
     */
496
    public function divide(int|float|matrix|vector $d): matrix {
497
        if ($d instanceof self) {
498
            return $this->divideMatrix($d);
499
        } elseif ($d instanceof vector) {
500
            return $this->divideVector($d);
501
        } else {
502
            return $this->divideScalar($d);
503
        }
504
    }
505
506
    protected function divideMatrix(matrix $m): matrix {
507
        if ($this->checkShape($this, $m) && $this->checkDtype($this, $m)) {
508
            $ar = self::factory($this->row, $this->col, $this->dtype);
509
            for ($i = 0; $i < $this->ndim; ++$i) {
510
                $ar->data[$i] = $this->data[$i] / $m->data[$i];
511
            }
512
            return $ar;
513
        }
514
    }
515
516
    protected function divideVector(vector $v): matrix {
517
        if ($this->checkDimensions($v, $this) && $this->checkDtype($this, $v)) {
518
            $ar = self::factory($this->row, $this->col, $this->dtype);
519
            for ($i = 0; $i < $this->row; ++$i) {
520
                for ($j = 0; $j < $this->col; ++$j) {
521
                    $ar->data[$i * $this->col + $j] = $this->data[$i * $this->col + $j] / $v->data[$j];
522
                }
523
            }
524
            return $ar;
525
        }
526
    }
527
528
    protected function divideScalar(int|float $s): matrix {
529
        $ar = self::factory($this->row, $this->col, $this->dtype);
530
        for ($i = 0; $i < $this->ndim; ++$i) {
531
            $ar->data[$i] = $this->data[$i] / $s;
532
        }
533
        return $ar;
534
    }
535
536
    /**
537
     * 
538
     * Raise this matrix to the power of the element-wise entry in another matrix.
539
     * 
540
     * @param int|float|matrix $m
541
     * @return matrix
542
     */
543
    public function pow(int|float|matrix|vector $d): matrix {
544
        if ($d instanceof self) {
545
            return $this->powMatrix($d);
546
        } else if ($d instanceof vector) {
547
            return $this->powVector($d);
548
        } else {
549
            return $this->powScalar($d);
550
        }
551
    }
552
553
    protected function powMatrix(matrix $m): matrix {
554
        if ($this->checkShape($this, $m) && $this->checkDtype($this, $m)) {
555
            $ar = self::factory($this->row, $this->col, $this->dtype);
556
            for ($i = 0; $i < $this->ndim; ++$i) {
557
                $ar->data[$i] = $this->data[$i] ** $m->data[$i];
558
            }
559
            return $ar;
560
        }
561
    }
562
563
    protected function powVector(vector $v): matrix {
564
        if ($this->checkDimensions($v, $this) && $this->checkDtype($this, $v)) {
565
            $ar = self::factory($this->row, $this->col, $this->dtype);
566
            for ($i = 0; $i < $this->row; ++$i) {
567
                for ($j = 0; $j < $this->col; ++$j) {
568
                    $ar->data[$i * $this->col + $j] = $this->data[$i * $this->col + $j] ** $v->data[$j];
569
                }
570
            }
571
            return $ar;
572
        }
573
    }
574
575
    protected function powScalar(int|float $s): matrix {
576
        $ar = $this->copy();
577
        for ($i = 0; $i < $this->ndim; ++$i) {
578
            $ar->data[$i] **= $s;
579
        }
580
        return $ar;
581
    }
582
583
    /**
584
     * Calculate the modulus i.e remainder of division between this matrix and another matrix.
585
     * @param int|float|matrix|vector $d
586
     * @return matrix
587
     */
588
    public function mod(int|float|matrix|vector $d): matrix {
589
        if ($d instanceof self) {
590
            $this->modMatrix($d);
591
        } else if ($d instanceof vector) {
592
            $this->modVector($d);
593
        } else {
594
            $this->modScalar($d);
595
        }
596
    }
597
598
    protected function modMatrix(matrix $m): matrix {
599
        if ($this->checkShape($this, $m) && $this->checkDtype($this, $m)) {
600
            $ar = self::factory($this->row, $this->col, $this->dtype);
601
            for ($i = 0; $i < $this->ndim; ++$i) {
602
                $ar->data[$i] = $this->data[$i] % $m->data[$i];
603
            }
604
            return $ar;
605
        } 
606
    }
607
608
    protected function modVector(vector $v): matrix {
609
        if ($this->checkDimensions($v, $this) && $this->checkDtype($this, $v)) {
610
            $ar = self::factory($this->row, $this->col, $this->dtype);
611
            for ($i = 0; $i < $this->row; ++$i) {
612
                for ($j = 0; $j < $this->col; ++$j) {
613
                    $ar->data[$i * $this->col + $j] = $this->data[$i * $this->col + $j] % $v->data[$j];
614
                }
615
            }
616
            return $ar;
617
        }
618
    }
619
620
    protected function modScalar(int|float $s): matrix {
621
        $ar = $this->copy();
622
        for ($i = 0; $i < $this->ndim; ++$i) {
623
            $ar->data[$i] %= $s;
624
        }
625
        return $ar;
626
    }
627
628
    /**
629
     * Return the element-wise reciprocal of the matrix.
630
     * 
631
     * @return matrix
632
     */
633
    public function reciprocal(): matrix {
634
        return self::ones($this->row, $this->col, $this->dtype)->divideMatrix($this);
635
    }
636
637
    /**
638
     * 
639
     * @param int|float $d
640
     * @return bool
641
     */
642
    public static function is_zero($d): bool {
643
        if (abs($d) < self::EPSILON) {
644
            return true;
645
        }
646
        return false;
647
    }
648
649
    /**
650
     *  is row zero
651
     * @param int $row
652
     * @return bool
653
     */
654
    public function is_rowZero(int $row): bool {
655
        for ($i = 0; $i < $this->col; ++$i) {
656
            if ($this->data[$row * $this->col + $i] != 0) {
657
                return false;
658
            }
659
        }
660
        return true;
661
    }
662
663
    /**
664
     * 
665
     * @return bool
666
     */
667
    public function has_ZeroRow(): bool {
668
        for ($i = 0; $i < $this->row; ++$i) {
669
            if ($this->is_rowZero($i)) {
670
                return true;
671
            }
672
        }
673
        return false;
674
    }
675
676
    /**
677
     * Transpose the matrix i.e row become cols and cols become rows.
678
     * @return \Np\matrix
679
     */
680
    public function transpose(): matrix {
681
        $ar = self::factory($this->col, $this->row, $this->dtype);
682
        for ($i = 0; $i < $ar->row; ++$i) {
683
            for ($j = 0; $j < $ar->col; ++$j) {
684
                $ar->data[$i * $ar->col + $j] = $this->data[$j * $ar->col + $i];
685
            }
686
        }
687
        return $ar;
688
    }
689
690
    /**
691
     * swap specific values in matrix
692
     * @param int $i1
693
     * @param int $i2
694
     */
695
    public function swapValue(int $i1, int $i2) {
696
        $tmp = $this->data[$i1];
697
        $this->data[$i1] = $this->data[$i2];
698
        $this->data[$i2] = $tmp;
699
    }
700
701
    /**
702
     * swap specific rows in matrix
703
     * @param int $r1
704
     * @param int $r2
705
     */
706
    public function swapRows(int $r1, int $r2) {
707
        for ($i = 0; $i < $this->col; ++$i) {
708
            $tmp = $this->data[$r1 * $this->col + $i];
709
            $this->data[$r1 * $this->col + $i] = $this->data[$r2 * $this->col + $i];
710
            $this->data[$r2 * $this->col + $i] = $tmp;
711
        }
712
    }
713
714
    /**
715
     * swap specific cols in matrix
716
     * @param int $c1
717
     * @param int $c2
718
     */
719
    public function swapCols(int $c1, int $c2) {
720
        for ($i = 0; $i < $this->row; ++$i) {
721
            $tmp = $this->data[$i * $this->row + $c1];
722
            $this->data[$i * $this->row + $c1] = $this->data[$i * $this->row + $c2];
723
            $this->data[$i * $this->row + $c2] = $tmp;
724
        }
725
    }
726
    
727
    /**
728
     * 
729
     * @param int|float $scalar
730
     * @return matrix
731
     */
732
    public function scale(int|float $scalar): matrix {
733
        if ($scalar == 0) {
734
            return self::zeros($this->row, $this->col, $this->dtype);
735
        }
736
737
        $ar = $this->copy();
738
        for ($i = 0; $i < $this->ndim; ++$i) {
739
            $ar->data[$i] *= $scalar;
740
        }
741
742
        return $ar;
743
    }
744
    
745
    /**
746
     * scale all the elements of a row 
747
     * @param int $row
748
     * @param int|float $c
749
     */
750
    public function scaleRow(int $row, int|float $c) {
751
        for ($i = 0; $i < $this->col; ++$i) {
752
            $this->data[$row * $this->col + $i] *= $c;
753
        }
754
    }
755
    
756
    /**
757
     * scale all the elements of 
758
     * @param int $col
759
     * @param int|float $c
760
     */
761
    public function scaleCol(int $col, int|float $c) {
762
        for ($i = 0; $i < $this->row; ++$i) {
763
            $this->data[$i * $this->col + $col] *= $c;
764
        }
765
    }
766
    
767
    /**
768
     * Scale digonally 
769
     * @param int|float $c
770
     * @param bool $lDig
771
     */
772
    public function scaleDigonalCol(int|float $c, bool $lDig = true) {
773
        if($lDig){
774
            for ($i = 0; $i < $this->row ; ++$i) {
775
                $this->data[$i * $this->col + $i] *= $c;
776
            }
777
        }
778
        else{
779
            for ($i = $this->row; $i > 0 ; --$i) {
780
                $this->data[$i * $this->col - $i] *= $c;
781
            }
782
        }
783
    }
784
785
    /**
786
     * 
787
     * @param int $r1
788
     * @param int $r2
789
     * @param float $c
790
     */
791
    public function addScaleRow(int $r1, int $r2, float $c) {
792
        for ($i = 0; $i < $this->col; ++$i) {
793
            $this->data[$r2 * $this->col + $i] += $this->data[$r1 * $this->col + $i] * $c;
794
        }
795
    }
796
797
    /**
798
     * Attach given matrix to the left of this matrix.
799
     * 
800
     * @param \Np\matrix $m
801
     * @return \Np\matrix
802
     */
803
    public function joinLeft(matrix $m): matrix {
804
        if ($this->row != $m->row && !$this->checkDtype($this, $m)) {
805
            self::_err('Error::Invalid size! or DataType!');
806
        }
807
        $col = $this->col + $m->col;
808
        $ar = self::factory($this->row, $col, $this->dtype);
809
        for ($i = 0; $i < $this->row; ++$i) {
810
            for ($j = 0; $j < $this->col; ++$j) {
811
                $ar->data[$i * $col + $j] = $this->data[$i * $this->col + $j];
812
            }
813
            for ($j = 0; $j < $m->col; ++$j) {
814
                $ar->data[$i * $col + ($this->col + $j)] = $m->data[$i * $m->col + $j];
815
            }
816
        }
817
        return $ar;
818
    }
819
820
    /**
821
     * Join matrix m to the Right of this matrix.
822
     * @param \Np\matrix $m
823
     * @return matrix
824
     */
825
    public function joinRight(matrix $m): matrix {
826
        if ($this->row != $m->row && !$this->checkDtype($this,$m)) {
827
            self::_err('Error::Invalid size! or DataType!');
828
        }
829
        $col = $this->col + $m->col;
830
        $ar = self::factory($this->row, $col, $this->dtype);
831
        for ($i = 0; $i < $m->row; ++$i) {
832
            for ($j = 0; $j < $m->col; ++$j) {
833
                $ar->data[$i * $col + $j] = $m->data[$i * $m->col + $j];
834
            }
835
            for ($j = 0; $j < $this->col; ++$j) {
836
                $ar->data[$i * $col + ($this->col + $j)] = $this->data[$i * $this->col + $j];
837
            }
838
        }
839
        return $ar;
840
    }
841
842
    /**
843
     * Join matrix m Above this matrix.
844
     * @param \Np\matrix $m
845
     * @return matrix
846
     */
847
    public function joinAbove(matrix $m): matrix {
848
        if ($this->col !== $m->col && !$this->checkDtype($this, $m)) {
849
            self::_err('Error::Invalid size! or DataType!');
850
        }
851
        $row = $this->row + $m->row;
852
        $ar = self::factory($row, $this->col, $this->dtype);
853
        for ($i = 0; $i < $m->row; ++$i) {
854
            for ($j = 0; $j < $m->col; ++$j) {
855
                $ar->data[$i * $m->col + $j] = $m->data[$i * $m->col + $j];
856
            }
857
            for ($j = 0; $j < $this->col; ++$j) {
858
                $ar->data[($i + $this->row) * $this->col + $j] = $this->data[$i * $this->col + $j];
859
            }
860
        }
861
        return $ar;
862
    }
863
864
    /**
865
     * Join matrix m below this matrix.
866
     * @param \Np\matrix $m
867
     * @return matrix
868
     */
869
    public function joinBelow(matrix $m): matrix {
870
        if ($this->col !== $m->col && !$this->checkDtype($this, $m)) {
871
            self::_err('Error::Invalid size! or DataType!');
872
        }
873
        $row = $this->row + $m->row;
874
        $ar = self::factory($row, $this->col, $this->dtype);
875
        for ($i = 0; $i < $this->row; ++$i) {
876
            for ($j = 0; $j < $this->col; ++$j) {
877
                $ar->data[$i * $this->col + $j] = $this->data[$i * $this->col + $j];
878
            }
879
            for ($j = 0; $j < $m->col; ++$j) {
880
                $ar->data[($i + $m->row) * $m->col + $j] = $m->data[$i * $m->col + $j];
881
            }
882
        }
883
        return $ar;
884
    }
885
886
    /**
887
     * Calculate the row echelon form of the matrix. 
888
     * Return the reduced matrix.
889
     *
890
     * @return matrix|null
891
     */
892
    public function ref(): matrix|null {
893
        return ref::factory($this);
894
    }
895
896
    /**
897
     * Return the lower triangular matrix of the Cholesky decomposition.
898
     *
899
     * @return matrix|null
900
     */
901
    public function cholesky(): matrix|null {
902
        return cholesky::factory($this);
903
    }
904
905
    /**
906
     * FIXME--------------
907
     * RREF
908
     * The reduced row echelon form (RREF) of a matrix.
909
     * @return \Np\matrix
910
     */
911
    public function rref(): matrix {
912
        return rref::factory($this);
913
    }
914
915
    /**
916
     * 
917
     * @param int $cols
918
     * @return \Np\matrix
919
     */
920
    public function diminish_left(int $cols): matrix {
921
        $ar = self::factory($this->row, $cols, $this->dtype);
922
        for ($i = 0; $i < $ar->row; ++$i) {
923
            for ($j = 0; $j < $ar->col; ++$j) {
924
                $ar->data[$i * $ar->col + $j] = $this->data[$i * $this->col + $j];
925
            }
926
        }
927
        return $ar;
928
    }
929
930
    /**
931
     * 
932
     * @param int $cols
933
     * @return \Np\matrix
934
     */
935
    public function diminish_right(int $cols): matrix {
936
        $ar = self::factory($this->row, $cols, $this->dtype);
937
        for ($i = 0; $i < $ar->row; ++$i) {
938
            for ($j = 0; $j < $ar->col; ++$j) {
939
                $ar->data[$i * $ar->col + $j] = $this->data[$i * $this->col - $cols + $j];
940
            }
941
        }
942
        return $ar;
943
    }
944
945
    /**
946
     * Return the index of the maximum element in every row of the matrix.
947
     * @return \Np\vector int
948
     */
949
    public function argMax(): vector {
950
        $v = vector::factory($this->row, vector::INT);
951
        for ($i = 0; $i < $this->row; ++$i) {
952
            $v->data[$i] = blas::max($this->rowAsVector($i));
953
        }
954
        return $v;
955
    }
956
957
    /**
958
     * Return the index of the minimum element in every row of the matrix.
959
     * @return \Np\vector int
960
     */
961
    public function argMin(): vector {
962
        $v = vector::factory($this->row, vector::INT);
963
        for ($i = 0; $i < $this->row; ++$i) {
964
            $v->data[$i] = blas::min($this->rowAsVector($i));
965
        }
966
967
        return $v;
968
    }
969
970
    /**
971
     * Set given data in matrix
972
     * @param int|float|array $data
973
     * @param bool $dignoal
974
     * @return void
975
     */
976
    public function setData(int|float|array $data, bool $dignoal = false): void {
977
        if ($dignoal == false) {
978
            if (is_array($data) && is_array($data[0])) {
979
                $f = $this->flattenArray($data);
980
                foreach ($f as $k => $v) {
981
                    $this->data[$k] = $v;
982
                }
983
            } elseif (is_numeric($data)) {
984
                for ($i = 0; $i < $this->ndim; ++$i) {
985
                    $this->data[$i] = $data;
986
                }
987
            }
988
        } elseif (is_numeric($data) || is_array($data) && !is_array($data[0])) {
989
            for ($i = 0; $i < $this->row; ++$i) {
990
                $this->data[$i * $this->col * $i] = $data;
991
            }
992
        }
993
    }
994
    
995
    /**
996
     * get the matrix data type
997
     * @return int
998
     */
999
    public function getDtype():int {
1000
        return $this->dtype;
1001
    }
1002
1003
    /**
1004
     * get the shape of matrix
1005
     * @return object
1006
     */
1007
    public function getShape(): object {
1008
        return (object) ['m' => $this->row, 'n' => $this->col];
1009
    }
1010
1011
    /**
1012
     * get the number of elements in the matrix.
1013
     * @return int
1014
     */
1015
    public function getSize(): int {
1016
        return $this->ndim;
1017
    }
1018
1019
    /**
1020
     * is matrix squred
1021
     * @return bool
1022
     */
1023
    public function isSquare(): bool {
1024
        if ($this->row === $this->col) {
1025
            return true;
1026
        }
1027
        return false;
1028
    }
1029
1030
    /**
1031
     * Return a row as vector from the matrix.
1032
     * @param int $index
1033
     * @return \Np\vector
1034
     */
1035
    public function rowAsVector(int $index): vector {
1036
        $vr = vector::factory($this->col, $this->dtype);
1037
        for ($j = 0; $j < $this->col; ++$j) {
1038
            $vr->data[$j] = $this->data[$index * $this->col + $j];
1039
        }
1040
        return $vr;
1041
    }
1042
1043
    /**
1044
     * Return a col as vector from the matrix.
1045
     * @param int $index
1046
     * @return \Np\vector
1047
     */
1048
    public function colAsVector(int $index): vector {
1049
        $vr = vector::factory($this->row, $this->dtype);
1050
        for ($i = 0; $i < $this->row; ++$i) {
1051
            $vr->data[$i] = $this->data[$i * $this->row + $index];
1052
        }
1053
        return $vr;
1054
    }
1055
1056
    /**
1057
     * Return the diagonal elements of a square matrix as a vector.
1058
     * @return \Np\vector
1059
     */
1060
    public function diagonalAsVector(): vector {
1061
        if (!$this->isSquare()) {
1062
            self::_err('Can not trace of a none square matrix');
1063
        }
1064
        $vr = vector::factory($this->row, $this->dtype);
1065
        for ($i = 0; $i < $this->row; ++$i) {
1066
            $vr->data[$i] = $this->getDiagonalVal($i);
1067
        }
1068
        return $vr;
1069
    }
1070
1071
    /**
1072
     * Flatten i.e unravel the matrix into a vector.
1073
     *
1074
     * @return \Np\vector
1075
     */
1076
    public function asVector(): vector {
1077
        $vr = vector::factory($this->ndim, $this->dtype);
1078
        for ($i = 0; $i < $this->ndim; ++$i) {
1079
            $vr->data[$i] = $this->data[$i];
1080
        }
1081
        return $vr;
1082
    }
1083
1084
    /**
1085
     * Return the elements of the matrix in a 2-d array.
1086
     * @return array
1087
     */
1088
    public function asArray(): array {
1089
        $ar = array_fill(0, $this->row, array_fill(0, $this->col, null));
1090
        for ($i = 0; $i < $this->row; ++$i) {
1091
            for ($j = 0; $j < $this->col; ++$j) {
1092
                $ar[$i][$j] = $this->data[$i * $this->col + $j];
1093
            }
1094
        }
1095
        return $ar;
1096
    }
1097
1098
    /**
1099
     * get a diagonal value from matrix
1100
     * @param int $i
1101
     * @return float
1102
     */
1103
    public function getDiagonalVal(int $i) {
1104
        if ($this->isSquare()) {
1105
            return $this->data[$i * $this->row + $i];
1106
        }
1107
    }
1108
1109
    
1110
1111
    /**
1112
     * Compute the singular value decomposition of a matrix and 
1113
     * return an object of the singular values and unitary matrices
1114
     *
1115
     * @return object (u,s,v)
1116
     */
1117
    public function svd(): svd {
1118
        return svd::factory($this);
1119
    }
1120
1121
    /**
1122
     * Compute the eigen decomposition of a general matrix.
1123
     * return the eigenvalues and eigenvectors as object
1124
     * 
1125
     * @param bool $symmetric
1126
     * @return eigen
1127
     */
1128
    public function eign(bool $symmetric = false): eigen {
1129
        return eigen::factory($this, $symmetric);
1130
    }
1131
1132
    /**
1133
     *  
1134
     * Compute the LU factorization of matrix.
1135
     * return lower, upper, and permutation matrices as object.
1136
     * 
1137
     * @return lu
1138
     */
1139
    public function lu(): lu {
1140
        return lu::factory($this);
1141
    }
1142
1143
    /**
1144
     * Return the L1 norm of the matrix.
1145
     * @return float
1146
     */
1147
    public function normL1(): float {
1148
        return lapack::lange('l', $this);
1149
    }
1150
1151
    /**
1152
     * Return the L2 norm of the matrix.
1153
     * @return float
1154
     */
1155
    public function normL2(): float {
1156
        return lapack::lange('f', $this);
1157
    }
1158
1159
    /**
1160
     * Return the L1 norm of the matrix.
1161
     * @return float
1162
     */
1163
    public function normINF(): float {
1164
        return lapack::lange('i', $this);
1165
    }
1166
1167
    /**
1168
     * Return the Frobenius norm of the matrix.
1169
     * @return float
1170
     */
1171
    public function normFrob(): float {
1172
        return $this->normL2();
1173
    }
1174
1175
    /**
1176
     * Compute the means of each row and return them in a vector.
1177
     *
1178
     * @return vector
1179
     */
1180
    public function mean(): vector {
1181
        return $this->sumRows()->divide($this->col);
1182
    }
1183
1184
    /**
1185
     * Compute the row variance of the matrix.
1186
     * 
1187
     * @param vector|null $mean
1188
     * @return vector
1189
     */
1190
    public function variance(vector|null $mean = null): vector {
1191
        if (isset($mean)) {
1192
            if (!$mean instanceof vector) {
1193
                self::_invalidArgument('mean must be a vector!');
1194
            }
1195
            if ($this->row !== $mean->col) {
1196
                self::_err('Err:: given mean vector dimensionality mismatched!');
1197
            }
1198
        } else {
1199
            $mean = $this->mean();
1200
        }
1201
        return $this->subtractColumnVector($mean)->square()
1202
                        ->sumRows()->divide($this->row);
1203
    }
1204
1205
    /**
1206
     *  Return the median vector of this matrix.
1207
     * @return vector
1208
     */
1209
    public function median(): vector {
1210
        $mid = intdiv($this->col, 2);
1211
        $odd = $this->col % 2 === 1;
1212
        $vr = vector::factory($this->row, $this->dtype);
1213
        for ($i = 0; $i < $this->row; ++$i) {
1214
            $a = $this->rowAsVector($i)->sort();
1215
            if ($odd) {
1216
                $median = $a->data[$mid];
1217
            } else {
1218
                $median = ($a->data[$mid - 1] + $a->data[$mid]) / 2.0;
1219
            }
1220
            $vr->data[$i] = $median;
1221
        }
1222
        unset($a);
1223
        return $vr;
1224
    }
1225
1226
    /**
1227
     * Compute the covariance matrix.
1228
     * 
1229
     * @param vector|null $mean
1230
     * @return matrix
1231
     */
1232
    public function covariance(vector|null $mean = null): matrix {
1233
        if (isset($mean)) {
1234
            if ($mean->col !== $this->row) {
1235
                self::_err('Err:: given mean vector dimensionality mismatched!');
1236
            }
1237
        } else {
1238
            $mean = $this->mean();
1239
        }
1240
1241
        $b = $this->subtractColumnVector($mean);
1242
1243
        return $b->dot($b->transpose())
1244
                        ->divideScalar($this->row);
1245
    }
1246
1247
    /**
1248
     * Square of matrix
1249
     * @return matrix
1250
     */
1251
    public function square(): matrix {
1252
        return $this->multiplyMatrix($this);
1253
    }
1254
1255
    /**
1256
     * 
1257
     * @param int|float|matrix|vector $d
1258
     * @return matrix
1259
     */
1260
    public function equal(int|float|matrix|vector $d): matrix {
1261
        if ($d instanceof self) {
1262
            return $this->equalMatrix($d);
1263
        } elseif ($d instanceof vector) {
1264
            return $this->equalVector($d);
1265
        } else {
1266
            return $this->equalScalar($d);
1267
        }
1268
    }
1269
1270
    protected function equalMatrix(matrix $m): matrix {
1271
        if ($this->checkShape($this, $m) && $this->checkDtype($this, $m)) {
1272
            $ar = self::factory($this->row, $this->col, $this->dtype);
1273
            for ($i = 0; $i < $this->ndim; ++$i) {
1274
                $ar->data[$i] = $this->data[$i] == $m->data[$i] ? 1 : 0;
1275
            }
1276
            return $ar;
1277
        }
1278
    }
1279
1280
    protected function equalVector(vector $v): matrix {
1281
        if ($this->checkDimensions($v, $this) && $this->checkDtype($this, $v)) {
1282
            $ar = self::factory($this->row, $this->col, $this->dtype);
1283
            for ($i = 0; $i < $this->row; ++$i) {
1284
                for ($j = 0; $j < $this->col; ++$j) {
1285
                    $ar->data[$i * $this->col + $j] = $this->data[$i * $this->col + $j] == $v->data[$j] ? 1 : 0;
1286
                }
1287
            }
1288
            return $ar;
1289
        }
1290
    }
1291
1292
    protected function equalScalar(int|float $s): matrix {
1293
        $ar = self::factory($this->row, $this->col, $this->dtype);
1294
        for ($i = 0; $i < $this->ndim; ++$i) {
1295
            $ar->data[$i] = $this->data[$i] == $s ? 1 : 0;
1296
        }
1297
        return $ar;
1298
    }
1299
1300
    /**
1301
     * 
1302
     * @param int|float|matrix|vector $d
1303
     * @return matrix
1304
     */
1305
    public function greater(int|float|matrix|vector $d): matrix {
1306
        if ($d instanceof self) {
1307
            return $this->greaterMatrix($d);
1308
        } elseif ($d instanceof vector) {
1309
            return $this->greaterVector($d);
1310
        } else {
1311
            return $this->greaterScalar($d);
1312
        }
1313
    }
1314
1315
    protected function greaterMatrix(matrix $m): matrix {
1316
        if ($this->checkShape($this, $m) && $this->checkDtype($this,$m)) {
1317
            $ar = self::factory($this->row, $this->col, $this->dtype);
1318
            for ($i = 0; $i < $this->ndim; ++$i) {
1319
                $ar->data[$i] = $this->data[$i] > $m->data[$i] ? 1 : 0;
1320
            }
1321
            return $ar;
1322
        }
1323
    }
1324
1325
    protected function greaterVector(vector $v): matrix {
1326
        if ($this->checkDimensions($v, $this) && $this->checkDtype($this,$v)) {
1327
            $ar = self::factory($this->row, $this->col, $this->dtype);
1328
            for ($i = 0; $i < $this->row; ++$i) {
1329
                for ($j = 0; $j < $this->col; ++$j) {
1330
                    $ar->data[$i * $this->col + $j] = $this->data[$i * $this->col + $j] > $v->data[$j] ? 1 : 0;
1331
                }
1332
            }
1333
            return $ar;
1334
        }
1335
    }
1336
1337
    protected function greaterScalar(int|float $s): matrix {
1338
        $ar = self::factory($this->row, $this->col, $this->dtype);
1339
        for ($i = 0; $i < $this->ndim; ++$i) {
1340
            $ar->data[$i] = $this->data[$i] > $s ? 1 : 0;
1341
        }
1342
        return $ar;
1343
    }
1344
1345
    /**
1346
     * 
1347
     * @param int|float|matrix $m
1348
     * @return matrix
1349
     */
1350
    public function less(int|float|matrix $m): matrix {
1351
        $ar = self::factory($this->row, $this->col, $this->dtype);
1352
        if ($m instanceof self) {
1353
            if ($this->checkShape($this,$m)) {
1354
                for ($i = 0; $i < $this->ndim; ++$i) {
1355
                    $ar->data[$i] = $this->data[$i] < $m->data[$i] ? 1 : 0;
1356
                }
1357
                return $ar;
1358
            }
1359
        } else {
1360
            for ($i = 0; $i < $this->ndim; ++$i) {
1361
                $ar->data[$i] = $this->data[$i] < $m ? 1 : 0;
1362
            }
1363
            return $ar;
1364
        }
1365
    }
1366
1367
    /**
1368
     * Is the matrix symmetric i.e. is it equal to its own transpose?
1369
     *
1370
     * @return bool
1371
     */
1372
    public function isSymmetric(): bool {
1373
        if (!$this->isSquare()) {
1374
            return false;
1375
        }
1376
        $ar = $this->transpose();
1377
        for ($i = 0; $i < $ar->ndim; ++$i) {
1378
            if ($ar->data[$i] != $this->data[$i]) {
1379
                unset($ar);
1380
                return false;
1381
            }
1382
        }
1383
        unset($ar);
1384
        return true;
1385
    }
1386
1387
    /**
1388
     * print the matrix in consol
1389
     */
1390
    public function printMatrix() {
1391
        echo __CLASS__ . PHP_EOL;
1392
        for ($i = 0; $i < $this->row; ++$i) {
1393
            for ($j = 0; $j < $this->col; ++$j) {
1394
                printf('%lf  ', $this->data[$i * $this->col + $j]);
1395
            }
1396
            echo PHP_EOL;
1397
        }
1398
    }
1399
1400
    public function __toString() {
1401
        return (string) $this->printMatrix();
1402
    }
1403
1404
    private function flattenArray(array $ar) {
1405
        if (is_array($ar) && is_array($ar[0])) {
1406
            $a = [];
1407
            foreach ($ar as $y => $value) {
1408
                foreach ($value as $k => $v) {
1409
                    $a[] = $v;
1410
                }
1411
            }
1412
            return $a;
1413
        }
1414
    }
1415
    
1416
    /**
1417
     * 
1418
     * @param int $row
1419
     * @param int $col
1420
     * @param int $dtype
1421
     * @return $this
1422
     */
1423
    protected function __construct(public int $row, public int $col, int $dtype = self::Float) {
1424
        if ($this->row < 1 || $this->col < 1) {
1425
            self::_invalidArgument('* To create Numphp/Matrix row & col must be greater than 0!, Op Failed! * ');
1426
        }
1427
        parent::__construct($this->row * $this->col, $dtype);
1428
        return $this;
1429
    }
1430
}
1431