basis/math/matrices/matrices.factor

   1 ! Copyright (C) 2005, 2010 Slava Pestov, Joe Groff.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: accessors arrays columns kernel locals math math.bits
   4 math.functions math.order math.vectors sequences
   5 sequences.private fry math.statistics grouping
   6 combinators.short-circuit math.ranges combinators.smart ;
   7 IN: math.matrices
   8
   9 ! Matrices
  10 : make-matrix ( m n quot -- matrix )
  11     '[ _ _ replicate ] replicate ; inline
  12
  13 : <matrix> ( m n element -- matrix )
  14     '[ _ _ <array> ] replicate ; inline
  15
  16 : zero-matrix ( m n -- matrix )
  17     0 <matrix> ; inline
  18
  19 : diagonal-matrix ( diagonal-seq -- matrix )
  20     dup length dup zero-matrix
  21     [ '[ dup _ nth set-nth ] each-index ] keep ; inline
  22
  23 : identity-matrix ( n -- matrix )
  24     1 <repetition> diagonal-matrix ; inline
  25
  26 : eye ( m n k -- matrix )
  27     [ [ iota ] bi@ ] dip neg '[ _ + = 1 0 ? ] cartesian-map ;
  28
  29 : hilbert-matrix ( m n -- matrix )
  30     [ iota ] bi@ [ + 1 + recip ] cartesian-map ;
  31
  32 : toeplitz-matrix ( n -- matrix )
  33     iota dup [ - abs 1 + ] cartesian-map ;
  34
  35 : hankel-matrix ( n -- matrix )
  36     [ iota dup ] keep '[ + abs 1 + dup _ > [ drop 0 ] when ] cartesian-map ;
  37
  38 : box-matrix ( r -- matrix )
  39     2 * 1 + dup '[ _ 1 <array> ] replicate ;
  40
  41 : vandermonde-matrix ( u n -- matrix )
  42     iota [ v^n ] with map reverse flip ;
  43
  44 :: rotation-matrix3 ( axis theta -- matrix )
  45     theta cos :> c
  46     theta sin :> s
  47     axis first3 :> ( x y z )
  48     x sq 1.0 x sq - c * +     x y * 1.0 c - * z s * -   x z * 1.0 c - * y s * + 3array
  49     x y * 1.0 c - * z s * +   y sq 1.0 y sq - c * +     y z * 1.0 c - * x s * - 3array
  50     x z * 1.0 c - * y s * -   y z * 1.0 c - * x s * +   z sq 1.0 z sq - c * +   3array
  51     3array ;
  52
  53 :: rotation-matrix4 ( axis theta -- matrix )
  54     theta cos :> c
  55     theta sin :> s
  56     axis first3 :> ( x y z )
  57     x sq 1.0 x sq - c * +     x y * 1.0 c - * z s * -   x z * 1.0 c - * y s * +   0 4array
  58     x y * 1.0 c - * z s * +   y sq 1.0 y sq - c * +     y z * 1.0 c - * x s * -   0 4array
  59     x z * 1.0 c - * y s * -   y z * 1.0 c - * x s * +   z sq 1.0 z sq - c * +     0 4array
  60     { 0.0 0.0 0.0 1.0 } 4array ;
  61
  62 :: translation-matrix4 ( offset -- matrix )
  63     offset first3 :> ( x y z )
  64     {
  65         { 1.0 0.0 0.0 x   }
  66         { 0.0 1.0 0.0 y   }
  67         { 0.0 0.0 1.0 z   }
  68         { 0.0 0.0 0.0 1.0 }
  69     } ;
  70
  71 : >scale-factors ( number/sequence -- x y z )
  72     dup number? [ dup dup ] [ first3 ] if ;
  73
  74 :: scale-matrix3 ( factors -- matrix )
  75     factors >scale-factors :> ( x y z )
  76     {
  77         { x   0.0 0.0 }
  78         { 0.0 y   0.0 }
  79         { 0.0 0.0 z   }
  80     } ;
  81
  82 :: scale-matrix4 ( factors -- matrix )
  83     factors >scale-factors :> ( x y z )
  84     {
  85         { x   0.0 0.0 0.0 }
  86         { 0.0 y   0.0 0.0 }
  87         { 0.0 0.0 z   0.0 }
  88         { 0.0 0.0 0.0 1.0 }
  89     } ;
  90
  91 : ortho-matrix4 ( dim -- matrix )
  92     [ recip ] map scale-matrix4 ;
  93
  94 :: frustum-matrix4 ( xy-dim near far -- matrix )
  95     xy-dim first2 :> ( x y )
  96     near x /f :> xf
  97     near y /f :> yf
  98     near far + near far - /f :> zf
  99     2 near far * * near far - /f :> wf
 100
 101     {
 102         { xf  0.0  0.0 0.0 }
 103         { 0.0 yf   0.0 0.0 }
 104         { 0.0 0.0  zf  wf  }
 105         { 0.0 0.0 -1.0 0.0 }
 106     } ;
 107
 108 :: skew-matrix4 ( theta -- matrix )
 109     theta tan :> zf
 110
 111     {
 112         { 1.0 0.0 0.0 0.0 }
 113         { 0.0 1.0 0.0 0.0 }
 114         { 0.0 zf  1.0 0.0 }
 115         { 0.0 0.0 0.0 1.0 }
 116     } ;
 117
 118 ! Matrix operations
 119 : mneg ( m -- m ) [ vneg ] map ;
 120
 121 : n+m  ( n m -- m ) [ n+v ] with map ;
 122 : m+n  ( m n -- m ) [ v+n ] curry map ;
 123 : n-m  ( n m -- m ) [ n-v ] with map ;
 124 : m-n  ( m n -- m ) [ v-n ] curry map ;
 125 : n*m ( n m -- m ) [ n*v ] with map ;
 126 : m*n ( m n -- m ) [ v*n ] curry map ;
 127 : n/m ( n m -- m ) [ n/v ] with map ;
 128 : m/n ( m n -- m ) [ v/n ] curry map ;
 129
 130 : m+   ( m m -- m ) [ v+ ] 2map ;
 131 : m-   ( m m -- m ) [ v- ] 2map ;
 132 : m*   ( m m -- m ) [ v* ] 2map ;
 133 : m/   ( m m -- m ) [ v/ ] 2map ;
 134
 135 : v.m ( v m -- v ) flip [ v. ] with map ;
 136 : m.v ( m v -- v ) [ v. ] curry map ;
 137 : m.  ( m m -- m ) flip [ swap m.v ] curry map ;
 138
 139 : m~  ( m m epsilon -- ? ) [ v~ ] curry 2all? ;
 140
 141 : mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ;
 142 : mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ;
 143 : mnorm ( m -- n ) dup mmax abs m/n ;
 144
 145 : cross ( vec1 vec2 -- vec3 )
 146     [ [ { 1 2 0 } vshuffle ] [ { 2 0 1 } vshuffle ] bi* v* ]
 147     [ [ { 2 0 1 } vshuffle ] [ { 1 2 0 } vshuffle ] bi* v* ] 2bi v- ; inline
 148
 149 :: normal ( vec1 vec2 vec3 -- vec4 )
 150     vec2 vec1 v- vec3 vec1 v- cross normalize ; inline
 151
 152 : proj ( v u -- w )
 153     [ [ v. ] [ norm-sq ] bi / ] keep n*v ;
 154
 155 : perp ( v u -- w )
 156     dupd proj v- ;
 157
 158 : angle-between ( v u -- a )
 159     [ normalize ] bi@ h. acos ;
 160
 161 : (gram-schmidt) ( v seq -- newseq )
 162     [ dupd proj v- ] each ;
 163
 164 : gram-schmidt ( seq -- orthogonal )
 165     V{ } clone [ over (gram-schmidt) suffix! ] reduce ;
 166
 167 : norm-gram-schmidt ( seq -- orthonormal )
 168     gram-schmidt [ normalize ] map ;
 169
 170 ERROR: negative-power-matrix m n ;
 171
 172 : (m^n) ( m n -- n )
 173     make-bits over first length identity-matrix
 174     [ [ dupd m. ] when [ dup m. ] dip ] reduce nip ;
 175
 176 : m^n ( m n -- n )
 177     dup 0 >= [ (m^n) ] [ negative-power-matrix ] if ;
 178
 179 : stitch ( m -- m' )
 180     [ ] [ [ append ] 2map ] map-reduce ;
 181
 182 : kron ( m1 m2 -- m )
 183     '[ [ _ n*m  ] map ] map stitch stitch ;
 184
 185 : outer ( u v -- m )
 186     [ n*v ] curry map ;
 187
 188 : row ( n matrix -- col )
 189     nth ; inline
 190
 191 : rows ( seq matrix -- cols )
 192     '[ _ row ] map ; inline
 193
 194 : col ( n matrix -- col )
 195     swap '[ _ swap nth ] map ; inline
 196
 197 : cols ( seq matrix -- cols )
 198     '[ _ col ] map ; inline
 199
 200 : set-index ( object pair matrix -- )
 201     [ first2 swap ] dip nth set-nth ; inline
 202
 203 : set-indices ( object sequence matrix -- )
 204     '[ _ set-index ] with each ; inline
 205
 206 : matrix-map ( matrix quot -- )
 207     '[ _ map ] map ; inline
 208
 209 : column-map ( matrix quot -- seq )
 210     [ [ first length iota ] keep ] dip '[ _ col @ ] map ; inline
 211
 212 : cartesian-square-indices ( n -- matrix )
 213     iota dup cartesian-product ; inline
 214
 215 : cartesian-matrix-map ( matrix quot -- matrix' )
 216     [ [ first length cartesian-square-indices ] keep ] dip
 217     '[ _ @ ] matrix-map ; inline
 218
 219 : cartesian-matrix-column-map ( matrix quot -- matrix' )
 220     [ cols first2 ] prepose cartesian-matrix-map ; inline
 221
 222 : cov-matrix-ddof ( matrix ddof -- cov )
 223     '[ _ cov-ddof ] cartesian-matrix-column-map ; inline
 224
 225 : cov-matrix ( matrix -- cov ) 0 cov-matrix-ddof ; inline
 226
 227 : sample-cov-matrix ( matrix -- cov ) 1 cov-matrix-ddof ; inline
 228
 229 GENERIC: square-rows ( object -- matrix )
 230 M: integer square-rows iota square-rows ;
 231 M: sequence square-rows
 232     [ length ] keep >array '[ _ clone ] { } replicate-as ;
 233
 234 GENERIC: square-cols ( object -- matrix )
 235 M: integer square-cols iota square-cols ;
 236 M: sequence square-cols
 237     [ length ] keep [ <array> ] with { } map-as ;
 238
 239 : make-matrix-with-indices ( m n quot -- matrix )
 240     [ [ iota ] bi@ ] dip '[ @ ] cartesian-map ; inline
 241
 242 : null-matrix? ( matrix -- ? ) empty? ; inline
 243
 244 : well-formed-matrix? ( matrix -- ? )
 245     [ t ] [
 246         [ ] [ first length ] bi
 247         '[ length _ = ] all?
 248     ] if-empty ;
 249
 250 : dim ( matrix -- pair/f )
 251     [ 2 0 <array> ]
 252     [ [ length ] [ first length ] bi 2array ] if-empty ;
 253
 254 : square-matrix? ( matrix -- ? )
 255     { [ well-formed-matrix? ] [ dim all-eq? ] } 1&& ;
 256
 257 : matrix-coordinates ( dim -- coordinates )
 258     first2 [ iota ] bi@ cartesian-product ; inline
 259
 260 : dimension-range ( matrix -- dim range )
 261     dim [ matrix-coordinates ] [ first [1,b] ] bi ;
 262
 263 : upper-matrix-indices ( matrix -- matrix' )
 264     dimension-range <reversed> [ tail-slice* >array ] 2map concat ;
 265
 266 : lower-matrix-indices ( matrix -- matrix' )
 267     dimension-range [ head-slice >array ] 2map concat ;
 268
 269 : make-lower-matrix ( object m n -- matrix )
 270     zero-matrix [ lower-matrix-indices ] [ set-indices ] [ ] tri ;
 271
 272 : make-upper-matrix ( object m n -- matrix )
 273     zero-matrix [ upper-matrix-indices ] [ set-indices ] [ ] tri ;