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 ;
10 : make-matrix ( m n quot -- matrix )
11 '[ _ _ replicate ] replicate ; inline
13 : <matrix> ( m n element -- matrix )
14 '[ _ _ <array> ] replicate ; inline
16 : zero-matrix ( m n -- matrix )
19 : diagonal-matrix ( diagonal-seq -- matrix )
20 dup length dup zero-matrix
21 [ '[ dup _ nth set-nth ] each-index ] keep ; inline
23 : identity-matrix ( n -- matrix )
24 1 <repetition> diagonal-matrix ; inline
26 : eye ( m n k -- matrix )
27 [ [ iota ] bi@ ] dip neg '[ _ + = 1 0 ? ] cartesian-map ;
29 : hilbert-matrix ( m n -- matrix )
30 [ iota ] bi@ [ + 1 + recip ] cartesian-map ;
32 : toeplitz-matrix ( n -- matrix )
33 iota dup [ - abs 1 + ] cartesian-map ;
35 : hankel-matrix ( n -- matrix )
36 [ iota dup ] keep '[ + abs 1 + dup _ > [ drop 0 ] when ] cartesian-map ;
38 : box-matrix ( r -- matrix )
39 2 * 1 + dup '[ _ 1 <array> ] replicate ;
41 : vandermonde-matrix ( u n -- matrix )
42 iota [ v^n ] with map reverse flip ;
44 :: rotation-matrix3 ( axis theta -- matrix )
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
53 :: rotation-matrix4 ( axis theta -- matrix )
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 ;
62 :: translation-matrix4 ( offset -- matrix )
63 offset first3 :> ( x y z )
71 : >scale-factors ( number/sequence -- x y z )
72 dup number? [ dup dup ] [ first3 ] if ;
74 :: scale-matrix3 ( factors -- matrix )
75 factors >scale-factors :> ( x y z )
82 :: scale-matrix4 ( factors -- matrix )
83 factors >scale-factors :> ( x y z )
91 : ortho-matrix4 ( dim -- matrix )
92 [ recip ] map scale-matrix4 ;
94 :: frustum-matrix4 ( xy-dim near far -- matrix )
95 xy-dim first2 :> ( x y )
98 near far + near far - /f :> zf
99 2 near far * * near far - /f :> wf
108 :: skew-matrix4 ( theta -- matrix )
119 : mneg ( m -- m ) [ vneg ] map ;
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 ;
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 ;
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 ;
139 : m~ ( m m epsilon -- ? ) [ v~ ] curry 2all? ;
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 ;
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
149 :: normal ( vec1 vec2 vec3 -- vec4 )
150 vec2 vec1 v- vec3 vec1 v- cross normalize ; inline
153 [ [ v. ] [ norm-sq ] bi / ] keep n*v ;
158 : angle-between ( v u -- a )
159 [ normalize ] bi@ h. acos ;
161 : (gram-schmidt) ( v seq -- newseq )
162 [ dupd proj v- ] each ;
164 : gram-schmidt ( seq -- orthogonal )
165 V{ } clone [ over (gram-schmidt) over push ] reduce ;
167 : norm-gram-schmidt ( seq -- orthonormal )
168 gram-schmidt [ normalize ] map ;
170 ERROR: negative-power-matrix m n ;
173 make-bits over first length identity-matrix
174 [ [ dupd m. ] when [ dup m. ] dip ] reduce nip ;
177 dup 0 >= [ (m^n) ] [ negative-power-matrix ] if ;
180 [ ] [ [ append ] 2map ] map-reduce ;
182 : kron ( m1 m2 -- m )
183 '[ [ _ n*m ] map ] map stitch stitch ;
188 : row ( n matrix -- col )
191 : rows ( seq matrix -- cols )
192 '[ _ row ] map ; inline
194 : col ( n matrix -- col )
195 swap '[ _ swap nth ] map ; inline
197 : cols ( seq matrix -- cols )
198 '[ _ col ] map ; inline
200 : set-index ( object pair matrix -- )
201 [ first2 swap ] dip nth set-nth ; inline
203 : set-indices ( object sequence matrix -- )
204 '[ _ set-index ] with each ; inline
206 : matrix-map ( matrix quot -- )
207 '[ _ map ] map ; inline
209 : column-map ( matrix quot -- seq )
210 [ [ first length iota ] keep ] dip '[ _ col @ ] map ; inline
212 : cartesian-square-indices ( n -- matrix )
213 iota dup cartesian-product ; inline
215 : cartesian-matrix-map ( matrix quot -- matrix' )
216 [ [ first length cartesian-square-indices ] keep ] dip
217 '[ _ @ ] matrix-map ; inline
219 : cartesian-matrix-column-map ( matrix quot -- matrix' )
220 [ cols first2 ] prepose cartesian-matrix-map ; inline
222 : cov-matrix-ddof ( matrix ddof -- cov )
223 '[ _ cov-ddof ] cartesian-matrix-column-map ; inline
225 : cov-matrix ( matrix -- cov ) 0 cov-matrix-ddof ; inline
227 : sample-cov-matrix ( matrix -- cov ) 1 cov-matrix-ddof ; inline
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 ;
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 ;
239 : make-matrix-with-indices ( m n quot -- matrix )
240 [ [ iota ] bi@ ] dip '[ @ ] cartesian-map ; inline
242 : null-matrix? ( matrix -- ? ) empty? ; inline
244 : well-formed-matrix? ( matrix -- ? )
246 [ ] [ first length ] bi
250 : dim ( matrix -- pair/f )
252 [ [ length ] [ first length ] bi 2array ] if-empty ;
254 : square-matrix? ( matrix -- ? )
255 { [ well-formed-matrix? ] [ dim all-eq? ] } 1&& ;
257 : matrix-coordinates ( dim -- coordinates )
258 first2 [ iota ] bi@ cartesian-product ; inline
260 : dimension-range ( matrix -- dim range )
261 dim [ matrix-coordinates ] [ first [1,b] ] bi ;
263 : upper-matrix-indices ( matrix -- matrix' )
264 dimension-range <reversed> [ tail-slice* >array ] 2map concat ;
266 : lower-matrix-indices ( matrix -- matrix' )
267 dimension-range [ head-slice >array ] 2map concat ;
269 : make-lower-matrix ( object m n -- matrix )
270 zero-matrix [ lower-matrix-indices ] [ set-indices ] [ ] tri ;
272 : make-upper-matrix ( object m n -- matrix )
273 zero-matrix [ upper-matrix-indices ] [ set-indices ] [ ] tri ;