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