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 : matrix-l-infinity-norm ( m -- n )
264 dup zero-matrix? [ drop 0 ] [
265 [ [ abs ] map-sum ] map supremum
266 ] if ; inline foldable
268 : matrix-l1-norm ( m -- n )
269 dup zero-matrix? [ drop 0 ] [
270 flip matrix-l-infinity-norm
271 ] if ; inline foldable
273 : matrix-l2-norm ( m -- n )
274 dup zero-matrix? [ drop 0 ] [
275 [ [ sq ] map-sum ] map-sum sqrt
276 ] if ; inline foldable
278 ! XXX: M: zero-matrix l1-norm drop 0 ; inline
279 ! XXX: M: matrix l1-norm matrix-l1-norm ; inline
281 ! XXX: M: zero-matrix l2-norm drop 0 ; inline
282 ! XXX: M: matrix l2-norm matrix-l2-norm ; inline
284 ! XXX: M: zero-matrix l-infinity-norm drop 0 ; inline
285 ! XXX: M: matrix l-infinity-norm matrix-l-infinity-norm ; inline
287 ALIAS: frobenius-norm matrix-l2-norm
288 ALIAS: hilbert-schmidt-norm matrix-l2-norm
290 :: matrix-p-q-norm ( m p q -- n )
291 m dup zero-matrix? [ drop 0 ] [
292 [ [ sq ] map-sum q p / ^ ] map-sum q recip ^
293 ] if ; inline foldable
295 : matrix-p-norm-entrywise ( m p -- n )
296 [ flatten1 V{ } like ] dip p-norm-default ; inline
298 ! XXX: M: zero-matrix p-norm-default 2drop 0 ; inline
299 ! XXX: M: matrix p-norm-default matrix-p-norm-entrywise ; inline
301 : matrix-p-norm ( m p -- n )
302 over zero-matrix? [ 2drop 0 ] [
304 { [ dup 1 number= ] [ drop matrix-l1-norm ] }
305 { [ dup 2 number= ] [ drop matrix-l2-norm ] }
306 { [ dup fp-infinity? ] [ drop matrix-l-infinity-norm ] }
307 [ matrix-p-norm-entrywise ]
309 ] if ; inline foldable
311 ! XXX: M: zero-matrix p-norm 2drop 0 ; inline
312 ! XXX: M: matrix p-norm matrix-p-norm ; inline
314 : matrix-normalize ( m -- m' )
317 ] unless ; inline foldable
319 ! well-defined for square matrices; but works on nonsquare too
320 : main-diagonal ( matrix -- seq )
321 >square-matrix [ swap nth-unsafe ] map-index ; inline
323 ! top right to bottom left; reverse the result if you expected it to start in the lower left
324 : anti-diagonal ( matrix -- seq )
325 >square-matrix [ swap nth-end-unsafe ] map-index ; inline
328 : (rows-iota) ( matrix -- rows-iota )
329 dimension first-unsafe <iota> ;
330 : (cols-iota) ( matrix -- cols-iota )
331 dimension second-unsafe <iota> ;
333 : simple-rows-except ( matrix desc quot -- others )
334 curry [ dup (rows-iota) ] dip
335 pick reject-as swap rows ; inline
337 : simple-cols-except ( matrix desc quot -- others )
338 curry [ dup (cols-iota) ] dip
339 pick reject-as swap cols transpose ; inline ! need to un-transpose the result of cols
341 CONSTANT: scalar-except-quot [ = ]
342 CONSTANT: sequence-except-quot [ member? ]
345 GENERIC: rows-except ( matrix desc -- others )
346 M: integer rows-except scalar-except-quot simple-rows-except ;
347 M: sequence rows-except sequence-except-quot simple-rows-except ;
349 GENERIC: cols-except ( matrix desc -- others )
350 M: integer cols-except scalar-except-quot simple-cols-except ;
351 M: sequence cols-except sequence-except-quot simple-cols-except ;
353 ! well-defined for any regular matrix
354 : matrix-except ( matrix exclude-pair -- submatrix )
355 first2 [ rows-except ] dip cols-except ;
357 ALIAS: submatrix-excluding matrix-except
359 :: matrix-except-all ( matrix -- submatrices )
360 matrix dimension [ <iota> ] map first2-unsafe cartesian-product
361 [ [ matrix swap matrix-except ] map ] map ;
363 ALIAS: all-submatrices matrix-except-all
365 : dimension ( matrix -- dimension )
367 [ [ length ] [ first-unsafe length ] bi 2array ] if-empty ;