1 ! Copyright (C) 2005, 2010, 2018, 2020 Slava Pestov, Joe Groff, and Cat Stevens.
2 ! See http://factorcode.org/license.txt for BSD license.
3 USING: accessors arrays classes.singleton columns combinators
4 combinators.short-circuit combinators.smart formatting fry
5 grouping kernel locals math math.bits math.functions math.order
6 math.private math.ranges math.statistics math.vectors
7 math.vectors.private sequences sequences.deep sequences.private
8 slots.private summary ;
11 ! defined here because of issue #1943
12 DEFER: regular-matrix?
13 : regular-matrix? ( object -- ? )
15 dup first-unsafe length
19 ! the MRO (class linearization) is performed in the order the predicates appear here
20 ! except that null-matrix is last (but it is relied upon by zero-matrix)
22 ! sequence > matrix > zero-matrix > square-matrix > zero-square-matrix > null-matrix
24 ! Factor bug that's hard to repro: using `bi and` in these predicates
25 ! instead of 1&& will cause spirious no-method and bounds-error errors in <square-cols>
26 ! and the tests/docs for no apparent reason
27 PREDICATE: matrix < sequence
28 { [ [ sequence? ] all? ] [ regular-matrix? ] } 1&& ;
30 PREDICATE: irregular-matrix < sequence
31 { [ [ sequence? ] all? ] [ regular-matrix? not ] } 1&& ;
34 ! can't define dim using this predicate for this reason,
35 ! unless we are going to write two versions of dim, one of which is generic
36 PREDICATE: square-matrix < matrix
37 dimension first2-unsafe = ;
39 PREDICATE: null-matrix < matrix
42 PREDICATE: zero-matrix < matrix
43 dup null-matrix? [ drop f ] [ flatten [ zero? ] all? ] if ;
45 PREDICATE: zero-square-matrix < square-matrix
46 { [ zero-matrix? ] [ square-matrix? ] } 1&& ;
48 ! Benign matrix constructors
49 : <matrix> ( m n element -- matrix )
50 '[ _ _ <array> ] replicate ; inline
52 : <matrix-by> ( m n quot: ( ... -- elt ) -- matrix )
53 '[ _ _ replicate ] replicate ; inline
55 : <matrix-by-indices> ( ... m n quot: ( ... m' n' -- ... elt ) -- ... matrix )
56 [ [ <iota> ] bi@ ] dip cartesian-map ; inline
58 : <zero-matrix> ( m n -- matrix )
61 : <zero-square-matrix> ( n -- matrix )
62 dup <zero-matrix> ; inline
65 : (nth-from-end) ( n seq -- n )
66 length 1 - swap - ; inline flushable
68 : nth-end ( n seq -- elt )
69 [ (nth-from-end) ] keep nth ; inline flushable
71 : nth-end-unsafe ( n seq -- elt )
72 [ (nth-from-end) ] keep nth-unsafe ; inline flushable
74 : array-nth-end-unsafe ( n seq -- elt )
75 [ (nth-from-end) ] keep swap 2 fixnum+fast slot ; inline flushable
77 : set-nth-end ( elt n seq -- )
78 [ (nth-from-end) ] keep set-nth ; inline
80 : set-nth-end-unsafe ( elt n seq -- )
81 [ (nth-from-end) ] keep set-nth-unsafe ; inline
84 ! main-diagonal matrix
85 : <diagonal-matrix> ( diagonal-seq -- matrix )
86 [ length <zero-square-matrix> ] keep over
87 '[ dup _ nth set-nth-unsafe ] each-index ; inline
89 ! could also be written slower as: <diagonal-matrix> [ reverse ] map
90 : <anti-diagonal-matrix> ( diagonal-seq -- matrix )
91 [ length <zero-square-matrix> ] keep over
92 '[ dup _ nth set-nth-end-unsafe ] each-index ; inline
94 : <identity-matrix> ( n -- matrix )
95 1 <repetition> <diagonal-matrix> ; inline
97 : <eye> ( m n k z -- matrix )
98 [ [ <iota> ] bi@ ] 2dip
100 cartesian-map ; inline
102 ! if m = n and k = 0 then <identity-matrix> is (possibly) more efficient
103 :: <simple-eye> ( m n k -- matrix )
105 [ n <identity-matrix> ]
106 [ m n k 1 <eye> ] if ; inline
108 : <coordinate-matrix> ( dim -- coordinates )
109 first2 [ <iota> ] bi@ cartesian-product ; inline
111 ALIAS: <cartesian-indices> <coordinate-matrix>
113 : <cartesian-square-indices> ( n -- matrix )
114 dup 2array <cartesian-indices> ; inline
116 ALIAS: transpose flip
119 : array-matrix? ( matrix -- ? )
121 [ [ array? ] all? ] bi and ; inline foldable flushable
123 : matrix-cols-iota ( matrix -- cols-iota )
124 first-unsafe length <iota> ; inline
126 : unshaped-cols-iota ( matrix -- cols-iota )
127 [ first-unsafe length 1 ] keep
128 [ length min ] (each) (each-integer) <iota> ; inline
130 : generic-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
131 [ <reversed> [ nth-end-unsafe ] with { } map-as ] curry { } map-as ; inline
133 : array-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
134 [ <reversed> [ array-nth-end-unsafe ] with map ] curry map ; inline
137 ! much faster than [ reverse ] map flip [ reverse ] map
138 : anti-transpose ( matrix -- newmatrix )
140 [ dup regular-matrix?
141 [ matrix-cols-iota ] [ unshaped-cols-iota ] if
145 array-anti-transpose-unsafe
147 generic-anti-transpose-unsafe
151 ALIAS: anti-flip anti-transpose
153 : row ( n matrix -- row )
156 : rows ( seq matrix -- rows )
157 '[ _ row ] map ; inline
159 : col ( n matrix -- col )
160 swap '[ _ swap nth ] map ; inline
162 : cols ( seq matrix -- cols )
163 '[ _ col ] map ; inline
165 :: >square-matrix ( m -- subset )
166 m dimension first2 :> ( x y ) {
168 { [ x y < ] [ x <iota> m cols transpose ] }
169 { [ x y > ] [ y <iota> m rows ] }
172 GENERIC: <square-rows> ( desc -- matrix )
173 M: integer <square-rows>
174 <iota> <square-rows> ;
175 M: sequence <square-rows>
176 [ length ] keep >array '[ _ clone ] { } replicate-as ;
178 M: square-matrix <square-rows> ;
179 M: matrix <square-rows> >square-matrix ; ! coercing to square is more useful than no-method
181 GENERIC: <square-cols> ( desc -- matrix )
182 M: integer <square-cols>
183 <iota> <square-cols> ;
184 M: sequence <square-cols>
187 M: square-matrix <square-cols> ;
188 M: matrix <square-cols>
191 <PRIVATE ! implementation details of <lower-matrix> and <upper-matrix>
192 : dimension-range ( matrix -- dim range )
193 dimension [ <coordinate-matrix> ] [ first [1,b] ] bi ;
195 : upper-matrix-indices ( matrix -- matrix' )
196 dimension-range <reversed> [ tail-slice* >array ] 2map concat ;
198 : lower-matrix-indices ( matrix -- matrix' )
199 dimension-range [ head-slice >array ] 2map concat ;
203 DEFER: matrix-set-nths
204 : <lower-matrix> ( object m n -- matrix )
205 <zero-matrix> [ lower-matrix-indices ] [ matrix-set-nths ] [ ] tri ;
207 : <upper-matrix> ( object m n -- matrix )
208 <zero-matrix> [ upper-matrix-indices ] [ matrix-set-nths ] [ ] tri ;
210 ! element- and sequence-wise operations, getters and setters
212 [ ] [ [ append ] 2map ] map-reduce ;
214 : matrix-map ( matrix quot: ( ... elt -- ... elt' ) -- matrix' )
215 '[ _ map ] map ; inline
217 : matrix-map-index ( matrix quot: ( ... elt i j -- ... elt' ) -- matrix' )
218 '[ [ swap @ ] curry map-index ] map-index ; inline
220 : column-map ( matrix quot: ( ... col -- ... col' ) -- matrix' )
221 [ transpose ] dip map transpose ; inline
223 : matrix-nth ( pair matrix -- elt )
224 [ first2 swap ] dip nth nth ; inline
226 : matrix-nths ( pairs matrix -- elts )
227 '[ _ matrix-nth ] map ; inline
229 : matrix-set-nth ( obj pair matrix -- )
230 [ first2 swap ] dip nth set-nth ; inline
232 : matrix-set-nths ( obj pairs matrix -- )
233 '[ _ matrix-set-nth ] with each ; inline
235 ! -------------------------------------------
236 ! simple math of matrices follows
237 : mneg ( m -- m' ) [ vneg ] map ;
238 : mabs ( m -- m' ) [ vabs ] map ;
240 : n+m ( n m -- m ) [ n+v ] with map ;
241 : m+n ( m n -- m ) [ v+n ] curry map ;
242 : n-m ( n m -- m ) [ n-v ] with map ;
243 : m-n ( m n -- m ) [ v-n ] curry map ;
244 : n*m ( n m -- m ) [ n*v ] with map ;
245 : m*n ( m n -- m ) [ v*n ] curry map ;
246 : n/m ( n m -- m ) [ n/v ] with map ;
247 : m/n ( m n -- m ) [ v/n ] curry map ;
249 : m+ ( m1 m2 -- m ) [ v+ ] 2map ;
250 : m- ( m1 m2 -- m ) [ v- ] 2map ;
251 : m* ( m1 m2 -- m ) [ v* ] 2map ;
252 : m/ ( m1 m2 -- m ) [ v/ ] 2map ;
254 : vdotm ( v m -- p ) flip [ vdot ] with map ;
255 : mdotv ( m v -- p ) [ vdot ] curry map ;
256 : mdot ( m m -- m ) flip [ swap mdotv ] curry map ;
258 : m~ ( m1 m2 epsilon -- ? ) [ v~ ] curry 2all? ;
260 : mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ;
261 : mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ;
263 : m-infinity-norm ( m -- n )
264 dup zero-matrix? [ drop 0. ] [
265 [ [ abs ] map-sum ] map supremum
269 dup zero-matrix? [ drop 0. ] [
273 : frobenius-norm ( m -- n )
274 dup zero-matrix? [ drop 0. ] [
275 [ [ sq ] map-sum ] map-sum sqrt
278 ALIAS: hilbert-schmidt-norm frobenius-norm
280 :: matrix-p-q-norm ( m p q -- n )
281 m dup zero-matrix? [ drop 0. ] [
282 [ [ sq ] map-sum q p / ^ ] map-sum q recip ^
285 : matrix-p-norm-entrywise ( m p -- n )
286 dup zero-matrix? [ 2drop 0. ] [
287 [ flatten1 V{ } like ] dip p-norm-default
290 : matrix-p-norm ( m p -- n )
291 dup zero-matrix? [ 2drop 0. ] [
293 { [ dup 1 = ] [ drop m-1norm ] }
294 { [ dup 2 = ] [ drop frobenius-norm ] }
295 { [ dup fp-infinity? ] [ drop m-infinity-norm ] }
296 [ matrix-p-norm-entrywise ]
300 : normalize-matrix ( m -- m' )
301 dup zero-matrix? [ ] [
305 ! well-defined for square matrices; but works on nonsquare too
306 : main-diagonal ( matrix -- seq )
307 >square-matrix [ swap nth-unsafe ] map-index ; inline
309 ! top right to bottom left; reverse the result if you expected it to start in the lower left
310 : anti-diagonal ( matrix -- seq )
311 >square-matrix [ swap nth-end-unsafe ] map-index ; inline
314 : (rows-iota) ( matrix -- rows-iota )
315 dimension first-unsafe <iota> ;
316 : (cols-iota) ( matrix -- cols-iota )
317 dimension second-unsafe <iota> ;
319 : simple-rows-except ( matrix desc quot -- others )
320 curry [ dup (rows-iota) ] dip
321 pick reject-as swap rows ; inline
323 : simple-cols-except ( matrix desc quot -- others )
324 curry [ dup (cols-iota) ] dip
325 pick reject-as swap cols transpose ; inline ! need to un-transpose the result of cols
327 CONSTANT: scalar-except-quot [ = ]
328 CONSTANT: sequence-except-quot [ member? ]
331 GENERIC: rows-except ( matrix desc -- others )
332 M: integer rows-except scalar-except-quot simple-rows-except ;
333 M: sequence rows-except sequence-except-quot simple-rows-except ;
335 GENERIC: cols-except ( matrix desc -- others )
336 M: integer cols-except scalar-except-quot simple-cols-except ;
337 M: sequence cols-except sequence-except-quot simple-cols-except ;
339 ! well-defined for any regular matrix
340 : matrix-except ( matrix exclude-pair -- submatrix )
341 first2 [ rows-except ] dip cols-except ;
343 ALIAS: submatrix-excluding matrix-except
345 :: matrix-except-all ( matrix -- submatrices )
346 matrix dimension [ <iota> ] map first2-unsafe cartesian-product
347 [ [ matrix swap matrix-except ] map ] map ;
349 ALIAS: all-submatrices matrix-except-all
351 : dimension ( matrix -- dimension )
353 [ [ length ] [ first-unsafe length ] bi 2array ] if-empty ;