basis/math/matrices/matrices.factor

   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 ;
   6 IN: math.matrices
   7
   8 ! defined here because of issue #1943
   9 DEFER: regular-matrix?
  10 : regular-matrix? ( object -- ? )
  11     [ t ] [
  12         dup first-unsafe length
  13         '[ length _ = ] all?
  14     ] if-empty ;
  15
  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)
  18 ! in other words:
  19 ! sequence > matrix > zero-matrix > square-matrix > zero-square-matrix > null-matrix
  20
  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&& ;
  26
  27 PREDICATE: irregular-matrix < sequence
  28     { [ [ sequence? ] all? ] [ regular-matrix? not ] } 1&& ;
  29
  30 DEFER: dimension
  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 = ;
  35
  36 PREDICATE: null-matrix < matrix
  37     flatten empty? ;
  38
  39 PREDICATE: zero-matrix < matrix
  40     dup null-matrix? [ drop f ] [ flatten [ zero? ] all? ] if ;
  41
  42 PREDICATE: zero-square-matrix < square-matrix
  43     { [ zero-matrix? ] [ square-matrix? ] } 1&& ;
  44
  45 ! Benign matrix constructors
  46 : <matrix> ( m n element -- matrix )
  47     '[ _ _ <array> ] replicate ; inline
  48
  49 : <matrix-by> ( m n quot: ( ... -- elt ) -- matrix )
  50     '[ _ _ replicate ] replicate ; inline
  51
  52 : <matrix-by-indices> ( ... m n quot: ( ... m' n' -- ... elt ) -- ... matrix )
  53     [ [ <iota> ] bi@ ] dip cartesian-map ; inline
  54
  55 : <zero-matrix> ( m n -- matrix )
  56     0 <matrix> ; inline
  57
  58 : <zero-square-matrix> ( n -- matrix )
  59     dup <zero-matrix> ; inline
  60
  61 <PRIVATE
  62 : (nth-from-end) ( n seq -- n )
  63     length 1 - swap - ; inline flushable
  64
  65 : nth-end ( n seq -- elt )
  66     [ (nth-from-end) ] keep nth ; inline flushable
  67
  68 : nth-end-unsafe ( n seq -- elt )
  69     [ (nth-from-end) ] keep nth-unsafe ; inline flushable
  70
  71 : array-nth-end-unsafe ( n seq -- elt )
  72     [ (nth-from-end) ] keep swap 2 fixnum+fast slot ; inline flushable
  73
  74 : set-nth-end ( elt n seq -- )
  75     [ (nth-from-end) ] keep set-nth ; inline
  76
  77 : set-nth-end-unsafe ( elt n seq -- )
  78     [ (nth-from-end) ] keep set-nth-unsafe ; inline
  79 PRIVATE>
  80
  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
  85
  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
  90
  91 : <identity-matrix> ( n -- matrix )
  92     1 <repetition> <diagonal-matrix> ; inline
  93
  94 : <eye> ( m n k z -- matrix )
  95     [ [ <iota> ] bi@ ] 2dip
  96     '[ _ neg + = _ 0 ? ]
  97     cartesian-map ; inline
  98
  99 ! if m = n and k = 0 then <identity-matrix> is (possibly) more efficient
 100 :: <simple-eye> ( m n k -- matrix )
 101     m n = k 0 = and
 102     [ n <identity-matrix> ]
 103     [ m n k 1 <eye> ] if ; inline
 104
 105 : <coordinate-matrix> ( dim -- coordinates )
 106   first2 [ <iota> ] bi@ cartesian-product ; inline
 107
 108 ALIAS: <cartesian-indices> <coordinate-matrix>
 109
 110 : <cartesian-square-indices> ( n -- matrix )
 111     dup 2array <cartesian-indices> ; inline
 112
 113 ALIAS: transpose flip
 114
 115 <PRIVATE
 116 : array-matrix? ( matrix -- ? )
 117     [ array? ]
 118     [ [ array? ] all? ] bi and ; inline foldable flushable
 119
 120 : matrix-cols-iota ( matrix -- cols-iota )
 121   first-unsafe length <iota> ; inline
 122
 123 : unshaped-cols-iota ( matrix -- cols-iota )
 124   [ first-unsafe length ] keep
 125   [ length min ] 1 sequence-iterator-from each-integer-from <iota> ; inline
 126
 127 : generic-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
 128     [ <reversed> [ nth-end-unsafe ] with { } map-as ] curry { } map-as ; inline
 129
 130 : array-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
 131     [ <reversed> [ array-nth-end-unsafe ] with map ] curry map ; inline
 132 PRIVATE>
 133
 134 ! much faster than [ reverse ] map flip [ reverse ] map
 135 : anti-transpose ( matrix -- newmatrix )
 136     dup empty? [ ] [
 137         [ dup regular-matrix?
 138             [ matrix-cols-iota ] [ unshaped-cols-iota ] if
 139         ] keep
 140
 141         dup array-matrix? [
 142             array-anti-transpose-unsafe
 143         ] [
 144             generic-anti-transpose-unsafe
 145         ] if
 146     ] if ;
 147
 148 ALIAS: anti-flip anti-transpose
 149
 150 : row ( n matrix -- row )
 151     nth ; inline
 152
 153 : rows ( seq matrix -- rows )
 154     '[ _ row ] map ; inline
 155
 156 : col ( n matrix -- col )
 157     swap '[ _ swap nth ] map ; inline
 158
 159 : cols ( seq matrix -- cols )
 160     '[ _ col ] map ; inline
 161
 162 :: >square-matrix ( m -- subset )
 163     m dimension first2 :> ( x y ) {
 164         { [ x y = ] [ m ] }
 165         { [ x y < ] [ x <iota> m cols transpose ] }
 166         { [ x y > ] [ y <iota> m rows ] }
 167     } cond ;
 168
 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 ;
 174
 175 M: square-matrix <square-rows> ;
 176 M: matrix <square-rows> >square-matrix ; ! coercing to square is more useful than no-method
 177
 178 GENERIC: <square-cols> ( desc -- matrix )
 179 M: integer <square-cols>
 180     <iota> <square-cols> ;
 181 M: sequence <square-cols>
 182     <square-rows> flip ;
 183
 184 M: square-matrix <square-cols> ;
 185 M: matrix <square-cols>
 186     >square-matrix ;
 187
 188 <PRIVATE ! implementation details of <lower-matrix> and <upper-matrix>
 189 : dimension-range ( matrix -- dim range )
 190     dimension [ <coordinate-matrix> ] [ first [1..b] ] bi ;
 191
 192 : upper-matrix-indices ( matrix -- matrix' )
 193     dimension-range <reversed> [ tail-slice* >array ] 2map concat ;
 194
 195 : lower-matrix-indices ( matrix -- matrix' )
 196     dimension-range [ head-slice >array ] 2map concat ;
 197 PRIVATE>
 198
 199 ! triangulars
 200 DEFER: matrix-set-nths
 201 : <lower-matrix> ( object m n -- matrix )
 202     <zero-matrix> [ lower-matrix-indices ] [ matrix-set-nths ] [ ] tri ;
 203
 204 : <upper-matrix> ( object m n -- matrix )
 205     <zero-matrix> [ upper-matrix-indices ] [ matrix-set-nths ] [ ] tri ;
 206
 207 ! element- and sequence-wise operations, getters and setters
 208 : stitch ( m -- m' )
 209     [ ] [ [ append ] 2map ] map-reduce ;
 210
 211 : matrix-map ( matrix quot: ( ... elt -- ... elt' ) -- matrix' )
 212     '[ _ map ] map ; inline
 213
 214 : matrix-map-index ( matrix quot: ( ... elt i j -- ... elt' ) -- matrix' )
 215     '[ [ swap @ ] curry map-index ] map-index ; inline
 216
 217 : column-map ( matrix quot: ( ... col -- ... col' ) -- matrix' )
 218     [ transpose ] dip map transpose ; inline
 219
 220 : matrix-nth ( pair matrix -- elt )
 221     [ first2 swap ] dip nth nth ; inline
 222
 223 : matrix-nths ( pairs matrix -- elts )
 224     '[ _ matrix-nth ] map ; inline
 225
 226 : matrix-set-nth ( obj pair matrix -- )
 227     [ first2 swap ] dip nth set-nth ; inline
 228
 229 : matrix-set-nths ( obj pairs matrix -- )
 230     '[ _ matrix-set-nth ] with each ; inline
 231
 232 ! -------------------------------------------
 233 ! simple math of matrices follows
 234 : mneg ( m -- m' ) [ vneg ] map ;
 235 : mabs ( m -- m' ) [ vabs ] map ;
 236
 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 ;
 245
 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 ;
 250
 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 ;
 254
 255 : m~  ( m1 m2 epsilon -- ? ) [ v~ ] curry 2all? ;
 256
 257 : mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ;
 258 : mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ;
 259
 260 : matrix-l-infinity-norm ( m -- n )
 261     dup zero-matrix? [ drop 0 ] [
 262         [ [ abs ] map-sum ] map supremum
 263     ] if ; inline foldable
 264
 265 : matrix-l1-norm ( m -- n )
 266     dup zero-matrix? [ drop 0 ] [
 267         flip matrix-l-infinity-norm
 268     ] if ; inline foldable
 269
 270 : matrix-l2-norm ( m -- n )
 271     dup zero-matrix? [ drop 0 ] [
 272         [ [ sq ] map-sum ] map-sum sqrt
 273     ] if ; inline foldable
 274
 275 ! XXX: M: zero-matrix l1-norm drop 0 ; inline
 276 ! XXX: M: matrix l1-norm matrix-l1-norm ; inline
 277
 278 ! XXX: M: zero-matrix l2-norm drop 0 ; inline
 279 ! XXX: M: matrix l2-norm matrix-l2-norm ; inline
 280
 281 ! XXX: M: zero-matrix l-infinity-norm drop 0 ; inline
 282 ! XXX: M: matrix l-infinity-norm matrix-l-infinity-norm ; inline
 283
 284 ALIAS: frobenius-norm matrix-l2-norm
 285 ALIAS: hilbert-schmidt-norm matrix-l2-norm
 286
 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
 291
 292 : matrix-p-norm-entrywise ( m p -- n )
 293     [ flatten1 V{ } like ] dip p-norm-default ; inline
 294
 295 ! XXX: M: zero-matrix p-norm-default 2drop 0 ; inline
 296 ! XXX: M: matrix p-norm-default matrix-p-norm-entrywise ; inline
 297
 298 : matrix-p-norm ( m p -- n )
 299     over zero-matrix? [ 2drop 0 ] [
 300         {
 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 ]
 305         } cond
 306     ] if ; inline foldable
 307
 308 ! XXX: M: zero-matrix p-norm 2drop 0 ; inline
 309 ! XXX: M: matrix p-norm matrix-p-norm ; inline
 310
 311 : matrix-normalize ( m -- m' )
 312     dup zero-matrix? [
 313         dup mabs mmax m/n
 314     ] unless ; inline foldable
 315
 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
 319
 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
 323
 324 <PRIVATE
 325 : (rows-iota) ( matrix -- rows-iota )
 326     dimension first-unsafe <iota> ;
 327 : (cols-iota) ( matrix -- cols-iota )
 328     dimension second-unsafe <iota> ;
 329
 330 : simple-rows-except ( matrix desc quot -- others )
 331     curry [ dup (rows-iota) ] dip
 332     pick reject-as swap rows ; inline
 333
 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
 337
 338 CONSTANT: scalar-except-quot [ = ]
 339 CONSTANT: sequence-except-quot [ member? ]
 340 PRIVATE>
 341
 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 ;
 345
 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 ;
 349
 350 ! well-defined for any regular matrix
 351 : matrix-except ( matrix exclude-pair -- submatrix )
 352     first2 [ rows-except ] dip cols-except ;
 353
 354 ALIAS: submatrix-excluding matrix-except
 355
 356 :: matrix-except-all ( matrix -- submatrices )
 357     matrix dimension [ <iota> ] map first2-unsafe cartesian-product
 358     [ [ matrix swap matrix-except ] map ] map ;
 359
 360 ALIAS: all-submatrices matrix-except-all
 361
 362 : dimension ( matrix -- dimension )
 363     [ { 0 0 } ]
 364     [ [ length ] [ first-unsafe length ] bi 2array ] if-empty ;