001 (ns org.clojars.punit-naik.clj-ml.utils.matrix
002 (:require [org.clojars.punit-naik.clj-ml.utils.generic :as gu]
003 [org.clojars.punit-naik.clj-ml.utils.linear-algebra :as lau]))
004
005 (defn equal-dimensions?
006 "Checks if the nested matrices of a matrix have euqal dimensions or not"
007 [m]
008 (if (every? coll? m)
009 (not (nil?
010 (reduce
011 (fn [acc r]
012 (if (nil? acc)
013 acc
014 (when (= (count acc) (count r)) r))) m)))
015 true))
016
017 (defn matrix?
018 "Returns true if a data structure is a valid N-D matrix, else false"
019 [m]
020 (when (coll? m)
021 (loop [mat m
022 nested? false
023 result false]
024 (if (or (and nested? (false? result)) (empty? mat))
025 result
026 (let [f (first mat)]
027 (if-not (coll? f)
028 (if nested? result (pos-int? (count mat)))
029 (recur (first mat) true
030 (and (equal-dimensions? mat)
031 (every? true? (map equal-dimensions? mat))))))))))
032
033 (defn get-val
034 "Get's a specific value from the martix `m` based on the path provided in `index-path`"
035 [m index-path]
036 (loop [ip index-path
037 result m]
038 (if (empty? ip)
039 result
040 (recur (rest ip) (nth result (first ip))))))
041
042 (defn index-matrix-rows
043 "This function indexes `matrix`'s rows and returns a map where key is the row number and value is the row itself"
044 [matrix]
045 (dissoc
046 (reduce
047 (fn [{:keys [size] :as acc} v]
048 (assoc acc (first size) v :size (rest size)))
049 {:size (range (count matrix))} matrix) :size))
050
051 (defmacro random-fn
052 "Executes the function `f` repeatedly `n` times"
053 [n f]
054 `(repeatedly ~n (fn [] ~f)))
055
056 (defn create-matrix
057 "Creates an N-D matrix of `dimenstions` specified in order with random float values in it.
058 Optional seed (Integer Value) for the random matrix generator can be specified as well."
059 [{:keys [dimensions seed] :or {seed 10}}]
060 (loop [d dimensions
061 result (lazy-seq [])]
062 (if (empty? d)
063 result
064 (recur (butlast d)
065 (if (seq result)
066 (random-fn (last d) result)
067 (concat result (random-fn (last d) (rand seed))))))))
068
069 (defn create-identity-matrix
070 "Same as `org.clojars.punit-naik.clj-ml.utils.matrix/create-matrix`
071 But this one creates a 2-D identity matrix"
072 [m]
073 (map (fn [i] (map (fn [j] (if (= i j) 1 0)) (range m))) (range m)))
074
075 (defn dimension
076 "Returns the dimenston of a 2-D matrix in a vector two elements"
077 [m]
078 [(count m) (count (first m))])
079
080 (defn square?
081 "Checks if a matrix is quare or not"
082 [matrix]
083 (let [[m n] (dimension matrix)]
084 (= m n)))
085
086 (defn identity-matrix?
087 "Checks if the matrix `m` is a 2-D identity matrix or not"
088 [m]
089 (let [[p] (dimension m)]
090 (and (matrix? m)
091 (square? m)
092 (every? true?
093 (flatten
094 (map (fn [i]
095 (map (fn [j]
096 (let [ij (get-val m [i j])]
097 (if (= i j)
098 (= ij (if (double? ij) 1.0 1))
099 (zero? ij))))
100 (range p)))
101 (range p)))))))
102
103 (defn transpose
104 "Returns the transpose of a 2-D matrix"
105 [m]
106 (apply map list m))
107
108 (defn perform-arithmetic-op
109 "Performs arithmetic operation on a matrix with a scalar or another matrix"
110 [mat operand operation]
111 (if-not (matrix? operand)
112 (map (fn [row] (map (fn [col-elem] (operation col-elem operand)) row)) mat)
113 (if (= (dimension mat) (dimension operand))
114 (map
115 (fn [row-1 row-2]
116 (map
117 operation
118 row-1 row-2)) mat operand)
119 (throw
120 (AssertionError. "Dimensions of matrices are not the same")))))
121
122 (defn matrix-multiply
123 "Multiplies two matrices of MxN and NxP dimensions
124 Calculates the dot product"
125 [a b]
126 (if (= (second (dimension a)) (first (dimension b)))
127 (let [b-transpose (transpose b)]
128 (map
129 (fn [row-a]
130 (map (fn [col-b] (reduce + (map * row-a col-b))) b-transpose)) a))
131 (throw
132 (AssertionError.
133 (str "The number of columns of the first matrix "
134 "are not equal to the number of rows of the second matrix")))))
135
136 (defn reciprocal
137 "Calculates the reciprocal of each and every element of a 2-D matrix"
138 [m]
139 (map (fn [row] (map (fn [col-elem] (double (/ 1 col-elem))) row)) m))
140
141 (defn exponential
142 "Calculates the exponential of each and every element of a 2-D matrix"
143 [m]
144 (map (fn [row] (map (fn [col-elem] (double (Math/exp col-elem))) row)) m))
145
146 (defn absolute
147 "Calculates the absolute value of each and every element of a 2-D matrix"
148 [m]
149 (map (fn [row] (map (fn [col-elem] (double (Math/abs col-elem))) row)) m))
150
151 (defn sigmoid
152 "Returns the sigmoid/logistic values of a 2-D matrix"
153 ([m]
154 (-> (perform-arithmetic-op m -1 *)
155 exponential
156 (perform-arithmetic-op 1 +)
157 reciprocal))
158 ([m _]
159 (-> (perform-arithmetic-op m -1 *)
160 (perform-arithmetic-op 1 +)
161 (perform-arithmetic-op m *))))
162
163 (defn swap-rows
164 "Swaps the rows at the `i`th and `j`th indexes"
165 [m i j]
166 (let [ith-row (get-val m [i])
167 jth-row (get-val m [j])]
168 (-> m
169 (gu/replace-nth i jth-row)
170 (gu/replace-nth j ith-row))))
171
172 (defn mean
173 "Calculates the mean of a 2-D matrix"
174 [m]
175 (let [averaged-rows (map (fn [row] (double (/ (reduce + row) (count row)))) m)]
176 (gu/mean-coll averaged-rows)))
177
178 (defn covariance
179 "Calculates the covariance matrix of a 2-D matrix"
180 [m]
181 (let [n (count m)
182 t (map (fn [i v] {:row i :value v}) (range n) (transpose m))]
183 (loop [tm t
184 result (lazy-seq [])]
185 (if (empty? tm)
186 (let [intermediate-result (reduce merge result)]
187 (reduce (fn [acc v]
188 (concat acc [(map #(get intermediate-result % 0.0) v)]))
189 (lazy-seq [])
190 (map (fn [i]
191 (map (fn [j] (vector i j))
192 (range (count t))))
193 (range (count t)))))
194 (let [{:keys [row value]} (first tm)
195 i-mean (gu/mean-coll value)]
196 (recur (rest tm)
197 (concat result
198 (map
199 (fn [x]
200 (let [j (:value (get-val t [x]))
201 j-mean (gu/mean-coll j)]
202 (assoc {} (sort [row x])
203 (double (/ (reduce +
204 (map #(* (- %1 i-mean) (- %2 j-mean))
205 value j))
206 (dec n))))))
207 (range row (count t))))))))))
208
209 (defn upper-triangular-matrix?
210 [m]
211 (->> m
212 (map-indexed
213 (fn [i row]
214 (if (zero? i)
215 true
216 (every? zero? (take i row)))))
217 (every? true?)))
218
219 (defn lower-triangular-matrix?
220 [m]
221 (->> m
222 (map-indexed
223 (fn [i row]
224 (if (= i (dec (count m)))
225 true
226 (every? zero? (take (- (dec (count m)) i) row)))))
227 (every? true?)))
228
229 (defn triangular-matrix?
230 [m]
231 (or (upper-triangular-matrix? m)
232 (lower-triangular-matrix? m)))
233
234 (defn row-adjust
235 "using `row-1` to adjust row elements of `row-2` so that their first `n` values are equal to zeros"
236 [row-1 row-2 n]
237 (if-not (every? zero? (take n row-2))
238 (loop [m (range (gu/first-n-zeros row-1) n)
239 r-2 row-2
240 prev-r-2 row-2]
241 (if (< (gu/first-n-zeros r-2) (gu/first-n-zeros prev-r-2))
242 prev-r-2
243 (if (or (>= (gu/first-n-zeros r-2) n)
244 (empty? m))
245 r-2
246 (let [mth-row-1 (double (nth row-1 (first m)))
247 mth-r-2 (double (nth r-2 (first m)))
248 mth-row-1-multiplier (Math/abs (/ mth-r-2 (if (zero? mth-row-1)
249 1 mth-row-1)))]
250 (recur (rest m)
251 (if-not (and (zero? mth-r-2) (zero? mth-row-1))
252 (as-> [row-1] $
253 (perform-arithmetic-op $ (Math/abs mth-row-1-multiplier) *)
254 (perform-arithmetic-op [r-2] $
255 (cond
256 (and (zero? mth-row-1)
257 (zero? mth-r-2)) +
258 (zero? mth-row-1) +
259 (zero? mth-r-2) +
260 :else (if (not= (/ mth-row-1 (Math/abs mth-row-1))
261 (/ mth-r-2 (Math/abs mth-r-2))) + -)))
262 (first $)
263 (map #(if (and (pos? %) (< % 0.0000001)) 0.0 %) $))
264 r-2)
265 r-2)))))
266 row-2))
267
268 (defn recursive-row-adjust
269 [matrix row-index-to-be-processed]
270 (if-not (and (= (dec (count matrix)) row-index-to-be-processed)
271 (every? zero? (nth matrix row-index-to-be-processed)))
272 (loop [row-idxs (range row-index-to-be-processed)
273 result nil]
274 (if (or (and (seq result)
275 (or (>= (gu/first-n-zeros result)
276 row-index-to-be-processed)
277 (>= (gu/first-n-zeros (get-val matrix [row-index-to-be-processed]))
278 row-index-to-be-processed)))
279 (empty? row-idxs))
280 result
281 (recur (rest row-idxs)
282 (if (or (= (first row-idxs) row-index-to-be-processed)
283 (and (zero? (get-val matrix [(first row-idxs) (dec row-index-to-be-processed)]))
284 (zero? (get-val matrix [(first row-idxs)
285 (gu/first-n-zeros
286 (get-val matrix [row-index-to-be-processed]))]))
287 (not (nil? result))))
288 result
289 (row-adjust (get-val matrix [(first row-idxs)])
290 (or result
291 (get-val matrix [row-index-to-be-processed]))
292 row-index-to-be-processed)))))
293 (nth matrix row-index-to-be-processed)))
294
295 (defn upper-triangular-matrix
296 "Converts any square matrix into an upper-triangular matrix
297 where all the matrix elements below the diagonal elements are zero"
298 [matrix]
299 (if (triangular-matrix? matrix)
300 {:upper-triangular matrix :num-swaps 0}
301 (let [[num-rows _] (dimension matrix)]
302 (loop [m matrix
303 nr (range 1 num-rows)
304 swap-count 0]
305 (if (empty? nr)
306 {:upper-triangular m :num-swaps swap-count}
307 (let [adjusted-row (recursive-row-adjust m (first nr))
308 first-n-zeros-adjusted-row (gu/first-n-zeros adjusted-row)
309 m-adjusted (gu/replace-nth m (first nr) adjusted-row)
310 should-be-swapped? (and (> first-n-zeros-adjusted-row (first nr))
311 (not= first-n-zeros-adjusted-row num-rows))
312 m-swapped (cond-> m-adjusted
313 should-be-swapped? (swap-rows (first nr) first-n-zeros-adjusted-row))]
314 (recur m-swapped
315 (cond-> nr
316 (not should-be-swapped?) rest)
317 (cond-> swap-count
318 should-be-swapped? inc))))))))
319
320 (defn determinant
321 "Caluclates the determinant of a square matrix by first calculating it's uper triangular matrix
322 and then multiplying it's diagonal elements together while taking into account the number of row swaps made in the process"
323 [matrix]
324 (let [{:keys [upper-triangular num-swaps]} (upper-triangular-matrix matrix)]
325 (-> (reduce
326 (fn [{:keys [i] :as acc} v]
327 (update (update acc :result #(* % (nth v (first i)))) :i rest))
328 {:i (range (count upper-triangular)) :result 1} upper-triangular)
329 :result (* (if (pos? num-swaps) -1.0 1.0)) gu/round-decimal)))
330
331 (defn cross-product
332 "Finds the cross product of two (indexed) rows of a matrix"
333 [row-1-indexed row-2-indexed]
334 (reduce #(merge-with + %1 %2)
335 (map (fn [x]
336 (reduce merge
337 (map
338 (fn [y] {(sort (flatten [(first x) (first y)])) (* (second x) (second y))})
339 row-2-indexed)))
340 row-1-indexed)))
341
342 (defn concat-matrix-rows
343 "Concatenates matrix rows with the first `(num-cols - 1)` values of the same row
344 This helps in finding the characteristic equation of the matrix"
345 [matrix num-cols]
346 (map #(concat % (take (dec num-cols) %)) matrix))
347
348 (defn characteristic-equation-parts
349 "Finds the positive or the negative parts of the characteristic equation
350 Both the positive and negative parts will be added to form the final eqution"
351 ([concatenated-matrix-minus-lambda-i num-rows]
352 (characteristic-equation-parts concatenated-matrix-minus-lambda-i num-rows true))
353 ([concatenated-matrix-minus-lambda-i num-rows positive-part?]
354 (let [end-col-idx (dec (second (dimension concatenated-matrix-minus-lambda-i)))
355 i-range (if positive-part?
356 (range num-rows)
357 (range end-col-idx (- end-col-idx num-rows) -1))]
358 (loop [i i-range
359 result []]
360 (if (empty? i)
361 result
362 (let [product (loop [paths (map vector
363 (range num-rows)
364 (range (first i)
365 ((if positive-part? + -) (first i) num-rows)
366 (if positive-part? 1 -1)))
367 result {}]
368 (if (empty? paths)
369 result
370 (recur (rest paths)
371 (let [ij (get-val concatenated-matrix-minus-lambda-i (first paths))
372 ij-mod (if (coll? ij) (index-matrix-rows ij) ij)]
373 (if (and (coll? result)
374 (empty? result))
375 ij-mod
376 (if (coll? result)
377 (if (coll? ij-mod)
378 (cross-product result ij-mod)
379 (into {} (map (fn [[k v]] [k (* v ij-mod)]) result)))
380 (if (coll? ij-mod)
381 (into {} (map (fn [[k v]] [k (* v result)]) ij-mod))
382 (* ij-mod result))))))))]
383 (recur (rest i)
384 (if (and (coll? result)
385 (empty? result))
386 (if (map? product)
387 (vals product) product)
388 (if (coll? result)
389 (if (coll? product)
390 (let [bigger-coll (if (>= (count product) (count result)) (vals product) result)
391 smaller-coll (if (< (count product) (count result)) (vals product) result)
392 smaller-coll-padded (concat (take (- (count bigger-coll) (count smaller-coll)) (repeat 0)) smaller-coll)]
393 (map + bigger-coll smaller-coll-padded))
394 (update (vec result) (dec (count result)) + product))
395 (if (coll? product)
396 (update (vec (vals product)) (dec (count product)) + result)
397 (+ result product)))))))))))
398
399 (defn matrix-minus-lambda-i
400 "Does A-λI when λ is known as well as unkonwn"
401 ([matrix] (matrix-minus-lambda-i matrix nil))
402 ([matrix lambda]
403 (map #(update (vec %1) %2 (fn [e]
404 (if (nil? lambda)
405 [-1 e] (- e lambda))))
406 matrix (range (first (dimension matrix))))))
407
408 (defn eigen-values
409 "Gets the eigen values of a matrix from it's characteristic equation"
410 [matrix]
411 (let [solve-eq (fn [eq]
412 (->> (remove zero? eq)
413 lau/solve-equation
414 (concat (filter zero? (last (partition-by identity eq))))
415 (map (fn [ev]
416 (if (re-find #"\." (str ev))
417 (gu/round-decimal ev) ev)))
418 sort))
419 [m n] (dimension matrix)
420 matrix-minus-unkown-lambda-i (matrix-minus-lambda-i matrix)]
421 (if (triangular-matrix? matrix)
422 (map (fn [i row]
423 (nth row i)) (range (count matrix)) matrix)
424 (let [concatenated-matrix-minus-lambda-i (concat-matrix-rows matrix-minus-unkown-lambda-i n)
425 first-part-product (characteristic-equation-parts concatenated-matrix-minus-lambda-i m)
426 second-part-product (characteristic-equation-parts concatenated-matrix-minus-lambda-i m false)
427 bigger-coll (if (>= (count first-part-product) (count second-part-product)) first-part-product second-part-product)
428 smaller-coll (if (< (count first-part-product) (count second-part-product)) first-part-product second-part-product)
429 smaller-coll-padded (concat (make-array Integer/TYPE (- (count bigger-coll)
430 (count smaller-coll))) smaller-coll)
431 smaller-coll-padded-negative (map #(* % -1) smaller-coll-padded)
432 eq (map + bigger-coll smaller-coll-padded-negative)]
433 (solve-eq eq)))))
434
435 (defn row-adjust-rref
436 "Adjusts the element of `row-1` at index `i` and makes it zero using the elements of `row-2`"
437 [row-1 row-2 i]
438 (let [row-1-i (nth row-1 i)]
439 (if-not (zero? row-1-i)
440 (let [row-1-i-abs (Math/abs row-1-i)
441 row-2-i (nth row-2 i)
442 row-2-i-abs (Math/abs row-2-i)
443 row-2-multiplier (/ row-1-i-abs row-2-i-abs)
444 minus-or-plus (if (= (double (/ row-1-i row-1-i-abs))
445 (double (/ row-2-i row-2-i-abs))) - +)]
446 (map (fn [row-1-elem row-2-elem]
447 (gu/round-decimal
448 (minus-or-plus row-1-elem row-2-elem)))
449 row-1 (map #(double (* % row-2-multiplier))
450 row-2)))
451 row-1)))
452
453 (defn zero-above-below-i-j
454 "Makes all the elements above and below `matrix[i, j]` zero using row transformations"
455 [matrix [i j] num-rows]
456 (let [matrix-i (nth matrix i)]
457 (loop [rows-to-be-fixed (concat (range 0 i) (range (inc i) num-rows))
458 result matrix]
459 (if (empty? rows-to-be-fixed)
460 result
461 (let [p (first rows-to-be-fixed)
462 fixed-row (row-adjust-rref (nth result p) matrix-i j)
463 all-zeros? (every? zero? fixed-row)]
464 (recur ((if all-zeros? butlast rest) rows-to-be-fixed)
465 (if all-zeros?
466 (concat (remove nil? (gu/replace-nth result p nil)) (list fixed-row))
467 (gu/replace-nth result p fixed-row))))))))
468
469 (defn reduced-row-echelon-form
470 [matrix]
471 (let [[m n] (dimension matrix)
472 indexes (range m)
473 sorted-matrix (sort-by #(vector (count (filter zero? %))
474 (gu/first-n-zeros %)) matrix)]
475 (loop [idxs indexes
476 result nil]
477 (if (empty? idxs)
478 result
479 (recur (rest idxs)
480 (let [i (first idxs)
481 row-i (nth (or result sorted-matrix) i)
482 all-zeros? (every? zero? row-i)
483 j (gu/first-n-zeros row-i)
484 ij-th (nth row-i (cond-> j
485 (= j n) dec))
486 new-result (cond-> (or result sorted-matrix)
487 (not all-zeros?) (zero-above-below-i-j [i j] m))]
488 (if (and (not all-zeros?)
489 (not= ij-th 1))
490 (gu/replace-nth new-result i (map (fn [e] (double (/ e (cond-> ij-th
491 (zero? ij-th) inc)))) row-i))
492 new-result)))))))
493
494 (defn adjust-rref-indices
495 "Adusts the indices of the RREF of a matrix by looking at the first <x> zeros of the rows"
496 [n rref]
497 (as-> (reduce (fn [acc row]
498 (let [index (gu/first-n-zeros row)
499 all-zeros? (every? zero? row)]
500 (cond-> acc
501 (not all-zeros?) (update :indices #(disj % index))
502 (not all-zeros?) (assoc index row))))
503 {:indices (into (hash-set) (range n))} rref) $
504 (reduce (fn [acc row]
505 (let [all-zeros? (every? zero? row)
506 index (first (acc :indices))]
507 (cond-> acc
508 all-zeros? (update :indices #(disj % index))
509 all-zeros? (assoc index row)))) $ rref)
510 (dissoc $ :indices)
511 (into (sorted-map) $)
512 (vals $)))
513
514 (defn eigen-vector-for-lamba
515 "Finds the eigenvector for a matrix with a particular eigenvalue"
516 ([matrix lambda] (eigen-vector-for-lamba matrix lambda false))
517 ([matrix lambda already-calculated?]
518 (let [[_ n] (dimension matrix)
519 rref-reversed (->> (matrix-minus-lambda-i matrix lambda)
520 reduced-row-echelon-form
521 (adjust-rref-indices n)
522 (map-indexed (fn [i row] [i row])) reverse)
523 zero-row-count (count (filter #(every? zero? (second %)) rref-reversed))]
524 (loop [rref-r rref-reversed
525 default-val-counter 1.0
526 result {}]
527 (if (empty? rref-r)
528 (cond->> (into (sorted-map) result)
529 (<= zero-row-count 1) vals
530 (> zero-row-count 1) (map (fn [[_ v]] ((if already-calculated? + -) v 1))))
531 (let [first-rref-reversed (first rref-r)
532 first-rref-reversed-index (first first-rref-reversed)
533 first-rref-reversed-row (second first-rref-reversed)
534 all-zeros? (every? zero? first-rref-reversed-row)]
535 (recur (rest rref-r) (cond-> default-val-counter
536 all-zeros? inc)
537 (cond-> result
538 all-zeros? (assoc first-rref-reversed-index (cond->> default-val-counter
539 already-calculated? (- 1.0)))
540 (not all-zeros?)
541 (assoc first-rref-reversed-index
542 (* -1.0 (reduce + (map
543 (fn [row-v i] (* row-v (get result i 0)))
544 first-rref-reversed-row (range n)))))))))))))
545
546 (defn eigen-vectors
547 "Finds the Eigenvectors of a matrix by using it's Eigenvalues"
548 [matrix eigen-values]
549 (loop [evals eigen-values
550 prev-eval nil
551 evecs []]
552 (if (empty? evals)
553 evecs
554 (recur (rest evals)
555 (first evals)
556 (conj evecs (eigen-vector-for-lamba matrix (first evals) (= prev-eval (first evals))))))))
557
558 (defn invertible?
559 "Checks if a matrix is invertible or not"
560 [matrix]
561 (and (square? matrix)
562 (not (zero? (determinant matrix)))))
563
564 (defn concat-identity-matrix
565 "Concatenates identity matrix of same size to the original matrix to generate [A|I]"
566 [matrix]
567 (map concat matrix (create-identity-matrix (first (dimension matrix)))))
568
569 (defn inverse
570 [matrix]
571 (when (invertible? matrix)
572 (let [matrix-cat-i (concat-identity-matrix matrix)
573 matrix-cat-i-rref (reduced-row-echelon-form matrix-cat-i)
574 [m] (dimension matrix)
575 left-part (map #(take m %) matrix-cat-i-rref)
576 right-part (map #(drop m %) matrix-cat-i-rref)]
577 (when (identity-matrix? left-part)
578 right-part))))
579
580 (defn rank
581 "Finds the rank of a mtrix by counting the number of non zero rows"
582 [m]
583 (->> m
584 (keep #(when-not (every? zero? %) true))
585 count))
586
587 (defn solve-linear-equation-system
588 "Given a set of linear equations where each equation is represented as:
589 ax + by + cz + ... = d
590 In a matrix for like:
591 [[a1 b1 c1 d1]
592 [a2 b2 c2 d2]
593 .
594 .
595 [an bn cn dn]]
596 Tries to find the solution for every unknown, x, y, z in this case"
597 [eq-matrix]
598 (let [rref-m (reduced-row-echelon-form eq-matrix)
599 a (map butlast rref-m)]
600 (when (= (rank rref-m)
601 (rank a))
602 (reduce
603 (fn [solution [i c-row]]
604 (let [a-row (butlast c-row)
605 b (last c-row)]
606 (assoc solution i
607 (float
608 (if (= i (dec (count eq-matrix)))
609 (/ b (last a-row))
610 (-> (->> (perform-arithmetic-op [(take-last i a-row)] [(take-last i solution)] *)
611 first
612 (apply +))
613 (* -1)
614 (+ b)
615 (/ (nth a-row i))))))))
616 (->> (first a)
617 count
618 (make-array Float/TYPE)
619 (into []))
620 (->> rref-m
621 (map-indexed vector)
622 reverse)))))