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: arrays combinators combinators.short-circuit kernel math
4 math.functions math.order math.private math.vectors ranges
5 sequences sequences.deep sequences.private slots.private ;
8 ! defined here because of issue #1943
10 : regular-matrix? ( object -- ? )
12 dup first-unsafe length
16 ! the MRO (class linearization) is performed in the order the predicates appear here
17 ! except that null-matrix is last (but it is relied upon by zero-matrix)
19 ! sequence > matrix > zero-matrix > square-matrix > zero-square-matrix > null-matrix
21 ! Factor bug that's hard to repro: using `bi and` in these predicates
22 ! instead of 1&& will cause spirious no-method and bounds-error errors in <square-cols>
23 ! and the tests/docs for no apparent reason
24 PREDICATE: matrix < sequence
25 { [ [ sequence? ] all? ] [ regular-matrix? ] } 1&& ;
27 PREDICATE: irregular-matrix < sequence
28 { [ [ sequence? ] all? ] [ regular-matrix? not ] } 1&& ;
31 ! can't define dim using this predicate for this reason,
32 ! unless we are going to write two versions of dim, one of which is generic
33 PREDICATE: square-matrix < matrix
34 dimension first2-unsafe = ;
36 PREDICATE: null-matrix < matrix
39 PREDICATE: zero-matrix < matrix
40 dup null-matrix? [ drop f ] [ flatten [ zero? ] all? ] if ;
42 PREDICATE: zero-square-matrix < square-matrix
43 { [ zero-matrix? ] [ square-matrix? ] } 1&& ;
45 ! Benign matrix constructors
46 : <matrix> ( m n element -- matrix )
47 '[ _ _ <array> ] replicate ; inline
49 : <matrix-by> ( m n quot: ( ... -- elt ) -- matrix )
50 '[ _ _ replicate ] replicate ; inline
52 : <matrix-by-indices> ( ... m n quot: ( ... m' n' -- ... elt ) -- ... matrix )
53 [ [ <iota> ] bi@ ] dip cartesian-map ; inline
55 : <zero-matrix> ( m n -- matrix )
58 : <zero-square-matrix> ( n -- matrix )
59 dup <zero-matrix> ; inline
62 : (nth-from-end) ( n seq -- n )
63 length 1 - swap - ; inline flushable
65 : nth-end ( n seq -- elt )
66 [ (nth-from-end) ] keep nth ; inline flushable
68 : nth-end-unsafe ( n seq -- elt )
69 [ (nth-from-end) ] keep nth-unsafe ; inline flushable
71 : array-nth-end-unsafe ( n seq -- elt )
72 [ (nth-from-end) ] keep swap 2 fixnum+fast slot ; inline flushable
74 : set-nth-end ( elt n seq -- )
75 [ (nth-from-end) ] keep set-nth ; inline
77 : set-nth-end-unsafe ( elt n seq -- )
78 [ (nth-from-end) ] keep set-nth-unsafe ; inline
81 ! main-diagonal matrix
82 : <diagonal-matrix> ( diagonal-seq -- matrix )
83 [ length <zero-square-matrix> ] keep over
84 '[ dup _ nth set-nth-unsafe ] each-index ; inline
86 ! could also be written slower as: <diagonal-matrix> [ reverse ] map
87 : <anti-diagonal-matrix> ( diagonal-seq -- matrix )
88 [ length <zero-square-matrix> ] keep over
89 '[ dup _ nth set-nth-end-unsafe ] each-index ; inline
91 : <identity-matrix> ( n -- matrix )
92 1 <repetition> <diagonal-matrix> ; inline
94 : <eye> ( m n k z -- matrix )
95 [ [ <iota> ] bi@ ] 2dip
97 cartesian-map ; inline
99 ! if m = n and k = 0 then <identity-matrix> is (possibly) more efficient
100 :: <simple-eye> ( m n k -- matrix )
102 [ n <identity-matrix> ]
103 [ m n k 1 <eye> ] if ; inline
105 : <coordinate-matrix> ( dim -- coordinates )
106 first2 [ <iota> ] bi@ cartesian-product ; inline
108 ALIAS: <cartesian-indices> <coordinate-matrix>
110 : <cartesian-square-indices> ( n -- matrix )
111 dup 2array <cartesian-indices> ; inline
113 ALIAS: transpose flip
116 : array-matrix? ( matrix -- ? )
118 [ [ array? ] all? ] bi and ; inline foldable flushable
120 : matrix-cols-iota ( matrix -- cols-iota )
121 first-unsafe length <iota> ; inline
123 : unshaped-cols-iota ( matrix -- cols-iota )
124 [ first-unsafe length ] keep
125 [ length min ] 1 (each-from) (each-integer) <iota> ; inline
127 : generic-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
128 [ <reversed> [ nth-end-unsafe ] with { } map-as ] curry { } map-as ; inline
130 : array-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
131 [ <reversed> [ array-nth-end-unsafe ] with map ] curry map ; inline
134 ! much faster than [ reverse ] map flip [ reverse ] map
135 : anti-transpose ( matrix -- newmatrix )
137 [ dup regular-matrix?
138 [ matrix-cols-iota ] [ unshaped-cols-iota ] if
142 array-anti-transpose-unsafe
144 generic-anti-transpose-unsafe
148 ALIAS: anti-flip anti-transpose
150 : row ( n matrix -- row )
153 : rows ( seq matrix -- rows )
154 '[ _ row ] map ; inline
156 : col ( n matrix -- col )
157 swap '[ _ swap nth ] map ; inline
159 : cols ( seq matrix -- cols )
160 '[ _ col ] map ; inline
162 :: >square-matrix ( m -- subset )
163 m dimension first2 :> ( x y ) {
165 { [ x y < ] [ x <iota> m cols transpose ] }
166 { [ x y > ] [ y <iota> m rows ] }
169 GENERIC: <square-rows> ( desc -- matrix )
170 M: integer <square-rows>
171 <iota> <square-rows> ;
172 M: sequence <square-rows>
173 [ length ] keep >array '[ _ clone ] { } replicate-as ;
175 M: square-matrix <square-rows> ;
176 M: matrix <square-rows> >square-matrix ; ! coercing to square is more useful than no-method
178 GENERIC: <square-cols> ( desc -- matrix )
179 M: integer <square-cols>
180 <iota> <square-cols> ;
181 M: sequence <square-cols>
184 M: square-matrix <square-cols> ;
185 M: matrix <square-cols>
188 <PRIVATE ! implementation details of <lower-matrix> and <upper-matrix>
189 : dimension-range ( matrix -- dim range )
190 dimension [ <coordinate-matrix> ] [ first [1..b] ] bi ;
192 : upper-matrix-indices ( matrix -- matrix' )
193 dimension-range <reversed> [ tail-slice* >array ] 2map concat ;
195 : lower-matrix-indices ( matrix -- matrix' )
196 dimension-range [ head-slice >array ] 2map concat ;
200 DEFER: matrix-set-nths
201 : <lower-matrix> ( object m n -- matrix )
202 <zero-matrix> [ lower-matrix-indices ] [ matrix-set-nths ] [ ] tri ;
204 : <upper-matrix> ( object m n -- matrix )
205 <zero-matrix> [ upper-matrix-indices ] [ matrix-set-nths ] [ ] tri ;
207 ! element- and sequence-wise operations, getters and setters
209 [ ] [ [ append ] 2map ] map-reduce ;
211 : matrix-map ( matrix quot: ( ... elt -- ... elt' ) -- matrix' )
212 '[ _ map ] map ; inline
214 : matrix-map-index ( matrix quot: ( ... elt i j -- ... elt' ) -- matrix' )
215 '[ [ swap @ ] curry map-index ] map-index ; inline
217 : column-map ( matrix quot: ( ... col -- ... col' ) -- matrix' )
218 [ transpose ] dip map transpose ; inline
220 : matrix-nth ( pair matrix -- elt )
221 [ first2 swap ] dip nth nth ; inline
223 : matrix-nths ( pairs matrix -- elts )
224 '[ _ matrix-nth ] map ; inline
226 : matrix-set-nth ( obj pair matrix -- )
227 [ first2 swap ] dip nth set-nth ; inline
229 : matrix-set-nths ( obj pairs matrix -- )
230 '[ _ matrix-set-nth ] with each ; inline
232 ! -------------------------------------------
233 ! simple math of matrices follows
234 : mneg ( m -- m' ) [ vneg ] map ;
235 : mabs ( m -- m' ) [ vabs ] map ;
237 : n+m ( n m -- m ) [ n+v ] with map ;
238 : m+n ( m n -- m ) [ v+n ] curry map ;
239 : n-m ( n m -- m ) [ n-v ] with map ;
240 : m-n ( m n -- m ) [ v-n ] curry map ;
241 : n*m ( n m -- m ) [ n*v ] with map ;
242 : m*n ( m n -- m ) [ v*n ] curry map ;
243 : n/m ( n m -- m ) [ n/v ] with map ;
244 : m/n ( m n -- m ) [ v/n ] curry map ;
246 : m+ ( m1 m2 -- m ) [ v+ ] 2map ;
247 : m- ( m1 m2 -- m ) [ v- ] 2map ;
248 : m* ( m1 m2 -- m ) [ v* ] 2map ;
249 : m/ ( m1 m2 -- m ) [ v/ ] 2map ;
251 : vdotm ( v m -- p ) flip [ vdot ] with map ;
252 : mdotv ( m v -- p ) [ vdot ] curry map ;
253 : mdot ( m m -- m ) flip [ swap mdotv ] curry map ;
255 : m~ ( m1 m2 epsilon -- ? ) [ v~ ] curry 2all? ;
257 : mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ;
258 : mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ;
260 : matrix-l-infinity-norm ( m -- n )
261 dup zero-matrix? [ drop 0 ] [
262 [ [ abs ] map-sum ] map supremum
263 ] if ; inline foldable
265 : matrix-l1-norm ( m -- n )
266 dup zero-matrix? [ drop 0 ] [
267 flip matrix-l-infinity-norm
268 ] if ; inline foldable
270 : matrix-l2-norm ( m -- n )
271 dup zero-matrix? [ drop 0 ] [
272 [ [ sq ] map-sum ] map-sum sqrt
273 ] if ; inline foldable
275 ! XXX: M: zero-matrix l1-norm drop 0 ; inline
276 ! XXX: M: matrix l1-norm matrix-l1-norm ; inline
278 ! XXX: M: zero-matrix l2-norm drop 0 ; inline
279 ! XXX: M: matrix l2-norm matrix-l2-norm ; inline
281 ! XXX: M: zero-matrix l-infinity-norm drop 0 ; inline
282 ! XXX: M: matrix l-infinity-norm matrix-l-infinity-norm ; inline
284 ALIAS: frobenius-norm matrix-l2-norm
285 ALIAS: hilbert-schmidt-norm matrix-l2-norm
287 :: matrix-p-q-norm ( m p q -- n )
288 m dup zero-matrix? [ drop 0 ] [
289 [ [ sq ] map-sum q p / ^ ] map-sum q recip ^
290 ] if ; inline foldable
292 : matrix-p-norm-entrywise ( m p -- n )
293 [ flatten1 V{ } like ] dip p-norm-default ; inline
295 ! XXX: M: zero-matrix p-norm-default 2drop 0 ; inline
296 ! XXX: M: matrix p-norm-default matrix-p-norm-entrywise ; inline
298 : matrix-p-norm ( m p -- n )
299 over zero-matrix? [ 2drop 0 ] [
301 { [ dup 1 number= ] [ drop matrix-l1-norm ] }
302 { [ dup 2 number= ] [ drop matrix-l2-norm ] }
303 { [ dup fp-infinity? ] [ drop matrix-l-infinity-norm ] }
304 [ matrix-p-norm-entrywise ]
306 ] if ; inline foldable
308 ! XXX: M: zero-matrix p-norm 2drop 0 ; inline
309 ! XXX: M: matrix p-norm matrix-p-norm ; inline
311 : matrix-normalize ( m -- m' )
314 ] unless ; inline foldable
316 ! well-defined for square matrices; but works on nonsquare too
317 : main-diagonal ( matrix -- seq )
318 >square-matrix [ swap nth-unsafe ] map-index ; inline
320 ! top right to bottom left; reverse the result if you expected it to start in the lower left
321 : anti-diagonal ( matrix -- seq )
322 >square-matrix [ swap nth-end-unsafe ] map-index ; inline
325 : (rows-iota) ( matrix -- rows-iota )
326 dimension first-unsafe <iota> ;
327 : (cols-iota) ( matrix -- cols-iota )
328 dimension second-unsafe <iota> ;
330 : simple-rows-except ( matrix desc quot -- others )
331 curry [ dup (rows-iota) ] dip
332 pick reject-as swap rows ; inline
334 : simple-cols-except ( matrix desc quot -- others )
335 curry [ dup (cols-iota) ] dip
336 pick reject-as swap cols transpose ; inline ! need to un-transpose the result of cols
338 CONSTANT: scalar-except-quot [ = ]
339 CONSTANT: sequence-except-quot [ member? ]
342 GENERIC: rows-except ( matrix desc -- others )
343 M: integer rows-except scalar-except-quot simple-rows-except ;
344 M: sequence rows-except sequence-except-quot simple-rows-except ;
346 GENERIC: cols-except ( matrix desc -- others )
347 M: integer cols-except scalar-except-quot simple-cols-except ;
348 M: sequence cols-except sequence-except-quot simple-cols-except ;
350 ! well-defined for any regular matrix
351 : matrix-except ( matrix exclude-pair -- submatrix )
352 first2 [ rows-except ] dip cols-except ;
354 ALIAS: submatrix-excluding matrix-except
356 :: matrix-except-all ( matrix -- submatrices )
357 matrix dimension [ <iota> ] map first2-unsafe cartesian-product
358 [ [ matrix swap matrix-except ] map ] map ;
360 ALIAS: all-submatrices matrix-except-all
362 : dimension ( matrix -- dimension )
364 [ [ length ] [ first-unsafe length ] bi 2array ] if-empty ;