steve
/
factor
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

math.matrices: rewrite, modernize and overhaul 

math.matrices.elimination: move to extra
math.matrices.extras: expand with esoteric, less-used and unfinished code from basis

- math.matrices and .extras receive more words, tests, and docs
- matrix has become a predicate class
- 94% of matrices words have complete docs
- 77% of matrices.extras words have complete docs
- much more consistent naming for constructors etc
- added missing words / features such as main-diagonal and anti-transpose
- optimizations
- lots of documentation
fix-linux
Cat Stevens 2018-02-04 00:00:21 -05:00 committed by John Benediktsson
parent c71b92eba9
commit 4350bcbfcd
27 changed files with 3462 additions and 650 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
basis/compression/run-length/run-length.factor
View File

@ -13,7 +13,7 @@ IN: compression.run-length
:: run-length-uncompress-bitmap4 ( byte-array m n -- byte-array' ) :: run-length-uncompress-bitmap4 ( byte-array m n -- byte-array' )
byte-array <sequence-parser> :> sp byte-array <sequence-parser> :> sp
m 1 + n zero-matrix :> matrix m 1 + n <zero-matrix> :> matrix
n 4 mod n + :> stride n 4 mod n + :> stride
0 :> i! 0 :> i!
0 :> j! 0 :> j!
@ -45,7 +45,7 @@ IN: compression.run-length
:: run-length-uncompress-bitmap8 ( byte-array m n -- byte-array' ) :: run-length-uncompress-bitmap8 ( byte-array m n -- byte-array' )
byte-array <sequence-parser> :> sp byte-array <sequence-parser> :> sp
m 1 + n zero-matrix :> matrix m 1 + n <zero-matrix> :> matrix
n 4 mod n + :> stride n 4 mod n + :> stride
0 :> i! 0 :> i!
0 :> j! 0 :> j!

3
basis/math/matrices/authors.txt
View File

@ -1 +1,4 @@
Slava Pestov Slava Pestov
Joe Groff
Doug Coleman
Cat Stevens

1170
basis/math/matrices/matrices-docs.factor
View File

File diff suppressed because it is too large Load Diff

932
basis/math/matrices/matrices-tests.factor
View File

File diff suppressed because it is too large Load Diff

468
basis/math/matrices/matrices.factor
View File

@ -1,122 +1,238 @@
! Copyright (C) 2005, 2010 Slava Pestov, Joe Groff. ! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens.
! See http://factorcode.org/license.txt for BSD license. ! See http://factorcode.org/license.txt for BSD license.
USING: accessors arrays columns kernel locals math math.bits USING: accessors arrays classes.singleton columns combinators
math.functions math.order math.vectors sequences combinators.short-circuit combinators.smart formatting fry
sequences.private fry math.statistics grouping grouping kernel locals math math.bits math.functions math.order
combinators.short-circuit math.ranges combinators.smart ; math.private math.ranges math.statistics math.vectors
math.vectors.private sequences sequences.deep sequences.private
slots.private summary ;
IN: math.matrices IN: math.matrices
! Matrices ! defined here because of issue #1943
: make-matrix ( m n quot -- matrix ) DEFER: regular-matrix?
'[ _ _ replicate ] replicate ; inline : regular-matrix? ( object -- ? )
[ t ] [
dup first-unsafe length
'[ length _ = ] all?
] if-empty ;
! the MRO (class linearization) is performed in the order the predicates appear here
! except that null-matrix is last (but it is relied upon by zero-matrix)
! in other words:
! sequence > matrix > zero-matrix > square-matrix > zero-square-matrix > null-matrix
! Factor bug that's hard to repro: using `bi and` in these predicates
! instead of 1&& will cause spirious no-method and bounds-error errors in <square-cols>
! and the tests/docs for no apparent reason
PREDICATE: matrix < sequence
{ [ [ sequence? ] all? ] [ regular-matrix? ] } 1&& ;
PREDICATE: irregular-matrix < sequence
{ [ [ sequence? ] all? ] [ regular-matrix? not ] } 1&& ;
DEFER: dimension
! can't define dim using this predicate for this reason,
! unless we are going to write two versions of dim, one of which is generic
PREDICATE: square-matrix < matrix
dimension first2-unsafe = ;
PREDICATE: null-matrix < matrix
flatten empty? ;
PREDICATE: zero-matrix < matrix
dup null-matrix? [ drop f ] [ flatten [ zero? ] all? ] if ;
PREDICATE: zero-square-matrix < square-matrix
{ [ zero-matrix? ] [ square-matrix? ] } 1&& ;
! Benign matrix constructors
: <matrix> ( m n element -- matrix ) : <matrix> ( m n element -- matrix )
'[ _ _ <array> ] replicate ; inline '[ _ _ <array> ] replicate ; inline
: zero-matrix ( m n -- matrix ) : <matrix-by> ( m n quot: ( ... -- elt ) -- matrix )
'[ _ _ replicate ] replicate ; inline
: <matrix-by-indices> ( ... m n quot: ( ... m' n' -- ... elt ) -- ... matrix )
[ [ <iota> ] bi@ ] dip cartesian-map ; inline
: <zero-matrix> ( m n -- matrix )
0 <matrix> ; inline 0 <matrix> ; inline
: diagonal-matrix ( diagonal-seq -- matrix ) : <zero-square-matrix> ( n -- matrix )
dup length dup zero-matrix dup <zero-matrix> ; inline
[ '[ dup _ nth set-nth ] each-index ] keep ; inline
: identity-matrix ( n -- matrix ) <PRIVATE
1 <repetition> diagonal-matrix ; inline : (nth-from-end) ( n seq -- n )
length 1 - swap - ; inline flushable
: eye ( m n k -- matrix ) : nth-end ( n seq -- elt )
[ [ <iota> ] bi@ ] dip neg '[ _ + = 1 0 ? ] cartesian-map ; [ (nth-from-end) ] keep nth ; inline flushable
: hilbert-matrix ( m n -- matrix ) : nth-end-unsafe ( n seq -- elt )
[ <iota> ] bi@ [ + 1 + recip ] cartesian-map ; [ (nth-from-end) ] keep nth-unsafe ; inline flushable
: toeplitz-matrix ( n -- matrix ) : array-nth-end-unsafe ( n seq -- elt )
<iota> dup [ - abs 1 + ] cartesian-map ; [ (nth-from-end) ] keep swap 2 fixnum+fast slot ; inline flushable
: hankel-matrix ( n -- matrix ) : set-nth-end ( elt n seq -- )
[ <iota> dup ] keep '[ + abs 1 + dup _ > [ drop 0 ] when ] cartesian-map ; [ (nth-from-end) ] keep set-nth ; inline
: box-matrix ( r -- matrix ) : set-nth-end-unsafe ( elt n seq -- )
2 * 1 + dup '[ _ 1 <array> ] replicate ; [ (nth-from-end) ] keep set-nth-unsafe ; inline
PRIVATE>
: vandermonde-matrix ( u n -- matrix ) ! main-diagonal matrix
<iota> [ v^n ] with map reverse flip ; : <diagonal-matrix> ( diagonal-seq -- matrix )
[ length <zero-square-matrix> ] keep over
'[ dup _ nth set-nth-unsafe ] each-index ; inline
:: rotation-matrix3 ( axis theta -- matrix ) ! could also be written slower as: <diagonal-matrix> [ reverse ] map
theta cos :> c : <anti-diagonal-matrix> ( diagonal-seq -- matrix )
theta sin :> s [ length <zero-square-matrix> ] keep over
axis first3 :> ( x y z ) '[ dup _ nth set-nth-end-unsafe ] each-index ; inline
x sq 1.0 x sq - c * + x y * 1.0 c - * z s * - x z * 1.0 c - * y s * + 3array
x y * 1.0 c - * z s * + y sq 1.0 y sq - c * + y z * 1.0 c - * x s * - 3array
x z * 1.0 c - * y s * - y z * 1.0 c - * x s * + z sq 1.0 z sq - c * + 3array
3array ;
:: rotation-matrix4 ( axis theta -- matrix ) : <identity-matrix> ( n -- matrix )
theta cos :> c 1 <repetition> <diagonal-matrix> ; inline
theta sin :> s
axis first3 :> ( x y z )
x sq 1.0 x sq - c * + x y * 1.0 c - * z s * - x z * 1.0 c - * y s * + 0 4array
x y * 1.0 c - * z s * + y sq 1.0 y sq - c * + y z * 1.0 c - * x s * - 0 4array
x z * 1.0 c - * y s * - y z * 1.0 c - * x s * + z sq 1.0 z sq - c * + 0 4array
{ 0.0 0.0 0.0 1.0 } 4array ;
:: translation-matrix4 ( offset -- matrix ) : <eye> ( m n k z -- matrix )
offset first3 :> ( x y z ) [ [ <iota> ] bi@ ] 2dip
{ '[ _ neg + = _ 0 ? ]
{ 1.0 0.0 0.0 x } cartesian-map ; inline
{ 0.0 1.0 0.0 y }
{ 0.0 0.0 1.0 z }
{ 0.0 0.0 0.0 1.0 }
} ;
: >scale-factors ( number/sequence -- x y z ) ! if m = n and k = 0 then <identity-matrix> is (possibly) more efficient
dup number? [ dup dup ] [ first3 ] if ; :: <simple-eye> ( m n k -- matrix )
m n = k 0 = and
[ n <identity-matrix> ]
[ m n k 1 <eye> ] if ; inline
:: scale-matrix3 ( factors -- matrix ) : <coordinate-matrix> ( dim -- coordinates )
factors >scale-factors :> ( x y z ) first2 [ <iota> ] bi@ cartesian-product ; inline
{
{ x 0.0 0.0 }
{ 0.0 y 0.0 }
{ 0.0 0.0 z }
} ;
:: scale-matrix4 ( factors -- matrix ) ALIAS: <cartesian-indices> <coordinate-matrix>
factors >scale-factors :> ( x y z )
{
{ x 0.0 0.0 0.0 }
{ 0.0 y 0.0 0.0 }
{ 0.0 0.0 z 0.0 }
{ 0.0 0.0 0.0 1.0 }
} ;
: ortho-matrix4 ( dim -- matrix ) : <cartesian-square-indices> ( n -- matrix )
[ recip ] map scale-matrix4 ; dup 2array <cartesian-indices> ; inline
:: frustum-matrix4 ( xy-dim near far -- matrix ) ALIAS: transpose flip
xy-dim first2 :> ( x y )
near x /f :> xf
near y /f :> yf
near far + near far - /f :> zf
2 near far * * near far - /f :> wf
{ <PRIVATE
{ xf 0.0 0.0 0.0 } : array-matrix? ( matrix -- ? )
{ 0.0 yf 0.0 0.0 } [ array? ]
{ 0.0 0.0 zf wf } [ [ array? ] all? ] bi and ; inline foldable flushable
{ 0.0 0.0 -1.0 0.0 }
} ;
:: skew-matrix4 ( theta -- matrix ) : matrix-cols-iota ( matrix -- cols-iota )
theta tan :> zf first-unsafe length <iota> ; inline
{ : unshaped-cols-iota ( matrix -- cols-iota )
{ 1.0 0.0 0.0 0.0 } [ first-unsafe length 1 ] keep
{ 0.0 1.0 0.0 0.0 } [ length min ] (each) (each-integer) <iota> ; inline
{ 0.0 zf 1.0 0.0 }
{ 0.0 0.0 0.0 1.0 }
} ;
! Matrix operations : generic-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
: mneg ( m -- m ) [ vneg ] map ; [ <reversed> [ nth-end-unsafe ] with { } map-as ] curry { } map-as ; inline
: array-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
[ <reversed> [ array-nth-end-unsafe ] with map ] curry map ; inline
PRIVATE>
! much faster than [ reverse ] map flip [ reverse ] map
: anti-transpose ( matrix -- newmatrix )
dup empty? [ ] [
[ dup regular-matrix?
[ matrix-cols-iota ] [ unshaped-cols-iota ] if
] keep
dup array-matrix? [
array-anti-transpose-unsafe
] [
generic-anti-transpose-unsafe
] if
] if ;
ALIAS: anti-flip anti-transpose
: row ( n matrix -- row )
nth ; inline
: rows ( seq matrix -- rows )
'[ _ row ] map ; inline
: col ( n matrix -- col )
swap '[ _ swap nth ] map ; inline
: cols ( seq matrix -- cols )
'[ _ col ] map ; inline
:: >square-matrix ( m -- subset )
m dimension first2 :> ( x y ) {
{ [ x y = ] [ m ] }
{ [ x y < ] [ x <iota> m cols transpose ] }
{ [ x y > ] [ y <iota> m rows ] }
} cond ;
GENERIC: <square-rows> ( desc -- matrix )
M: integer <square-rows>
<iota> <square-rows> ;
M: sequence <square-rows>
[ length ] keep >array '[ _ clone ] { } replicate-as ;
M: square-matrix <square-rows> ;
M: matrix <square-rows> >square-matrix ; ! coercing to square is more useful than no-method
GENERIC: <square-cols> ( desc -- matrix )
M: integer <square-cols>
<iota> <square-cols> ;
M: sequence <square-cols>
<square-rows> flip ;
M: square-matrix <square-cols> ;
M: matrix <square-cols>
>square-matrix ;
<PRIVATE ! implementation details of <lower-matrix> and <upper-matrix>
: dimension-range ( matrix -- dim range )
dimension [ <coordinate-matrix> ] [ first [1,b] ] bi ;
: upper-matrix-indices ( matrix -- matrix' )
dimension-range <reversed> [ tail-slice* >array ] 2map concat ;
: lower-matrix-indices ( matrix -- matrix' )
dimension-range [ head-slice >array ] 2map concat ;
PRIVATE>
! triangulars
DEFER: matrix-set-nths
: <lower-matrix> ( object m n -- matrix )
<zero-matrix> [ lower-matrix-indices ] [ matrix-set-nths ] [ ] tri ;
: <upper-matrix> ( object m n -- matrix )
<zero-matrix> [ upper-matrix-indices ] [ matrix-set-nths ] [ ] tri ;
! element- and sequence-wise operations, getters and setters
: stitch ( m -- m' )
[ ] [ [ append ] 2map ] map-reduce ;
: matrix-map ( matrix quot: ( ... elt -- ... elt' ) -- matrix' )
'[ _ map ] map ; inline
: column-map ( matrix quot: ( ... col -- ... col' ) -- matrix' )
[ transpose ] dip map transpose ; inline
: matrix-nth ( pair matrix -- elt )
[ first2 swap ] dip nth nth ; inline
: matrix-nths ( pairs matrix -- elts )
'[ _ matrix-nth ] map ; inline
: matrix-set-nth ( obj pair matrix -- )
[ first2 swap ] dip nth set-nth ; inline
: matrix-set-nths ( obj pairs matrix -- )
'[ _ matrix-set-nth ] with each ; inline
! -------------------------------------------
! simple math of matrices follows
: mneg ( m -- m' ) [ vneg ] map ;
: mabs ( m -- m' ) [ vabs ] map ;
: n+m ( n m -- m ) [ n+v ] with map ; : n+m ( n m -- m ) [ n+v ] with map ;
: m+n ( m n -- m ) [ v+n ] curry map ; : m+n ( m n -- m ) [ v+n ] curry map ;
@ -127,150 +243,70 @@ IN: math.matrices
: n/m ( n m -- m ) [ n/v ] with map ; : n/m ( n m -- m ) [ n/v ] with map ;
: m/n ( m n -- m ) [ v/n ] curry map ; : m/n ( m n -- m ) [ v/n ] curry map ;
: m+ ( m m -- m ) [ v+ ] 2map ; : m+ ( m1 m2 -- m ) [ v+ ] 2map ;
: m- ( m m -- m ) [ v- ] 2map ; : m- ( m1 m2 -- m ) [ v- ] 2map ;
: m* ( m m -- m ) [ v* ] 2map ; : m* ( m1 m2 -- m ) [ v* ] 2map ;
: m/ ( m m -- m ) [ v/ ] 2map ; : m/ ( m1 m2 -- m ) [ v/ ] 2map ;
: v.m ( v m -- v ) flip [ v. ] with map ; : v.m ( v m -- p ) flip [ v. ] with map ;
: m.v ( m v -- v ) [ v. ] curry map ; : m.v ( m v -- p ) [ v. ] curry map ;
: m. ( m m -- m ) flip [ swap m.v ] curry map ; : m. ( m m -- m ) flip [ swap m.v ] curry map ;
: m~ ( m m epsilon -- ? ) [ v~ ] curry 2all? ; : m~ ( m1 m2 epsilon -- ? ) [ v~ ] curry 2all? ;
: mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ; : mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ;
: mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ; : mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ;
: mnorm ( m -- n ) dup mmax abs m/n ; : mnorm ( m -- m' ) dup mmax abs m/n ;
: m-infinity-norm ( m -- n ) [ [ abs ] map-sum ] map supremum ; : m-infinity-norm ( m -- n ) [ [ abs ] map-sum ] map supremum ;
: m-1norm ( m -- n ) flip m-infinity-norm ; : m-1norm ( m -- n ) flip m-infinity-norm ;
: frobenius-norm ( m -- n ) [ [ sq ] map-sum ] map-sum sqrt ; : frobenius-norm ( m -- n ) [ [ sq ] map-sum ] map-sum sqrt ;
: cross ( vec1 vec2 -- vec3 ) ! well-defined for square matrices; but works on nonsquare too
[ [ { 1 2 0 } vshuffle ] [ { 2 0 1 } vshuffle ] bi* v* ] : main-diagonal ( matrix -- seq )
[ [ { 2 0 1 } vshuffle ] [ { 1 2 0 } vshuffle ] bi* v* ] 2bi v- ; inline >square-matrix [ swap nth-unsafe ] map-index ; inline
:: normal ( vec1 vec2 vec3 -- vec4 ) ! top right to bottom left; reverse the result if you expected it to start in the lower left
vec2 vec1 v- vec3 vec1 v- cross normalize ; inline : anti-diagonal ( matrix -- seq )
>square-matrix [ swap nth-end-unsafe ] map-index ; inline
: proj ( v u -- w ) <PRIVATE
[ [ v. ] [ norm-sq ] bi / ] keep n*v ; : (rows-iota) ( matrix -- rows-iota )
dimension first-unsafe <iota> ;
: (cols-iota) ( matrix -- cols-iota )
dimension second-unsafe <iota> ;
: perp ( v u -- w ) : simple-rows-except ( matrix desc quot -- others )
dupd proj v- ; curry [ dup (rows-iota) ] dip
pick reject-as swap rows ; inline
: angle-between ( v u -- a ) : simple-cols-except ( matrix desc quot -- others )
[ normalize ] bi@ h. acos ; curry [ dup (cols-iota) ] dip
pick reject-as swap cols transpose ; inline ! need to un-transpose the result of cols
: (gram-schmidt) ( v seq -- newseq ) CONSTANT: scalar-except-quot [ = ]
[ dupd proj v- ] each ; CONSTANT: sequence-except-quot [ member? ]
PRIVATE>
: gram-schmidt ( seq -- orthogonal ) GENERIC: rows-except ( matrix desc -- others )
V{ } clone [ over (gram-schmidt) suffix! ] reduce ; M: integer rows-except scalar-except-quot simple-rows-except ;
M: sequence rows-except sequence-except-quot simple-rows-except ;
: norm-gram-schmidt ( seq -- orthonormal ) GENERIC: cols-except ( matrix desc -- others )
gram-schmidt [ normalize ] map ; M: integer cols-except scalar-except-quot simple-cols-except ;
M: sequence cols-except sequence-except-quot simple-cols-except ;
ERROR: negative-power-matrix m n ; ! well-defined for any regular matrix
: matrix-except ( matrix exclude-pair -- submatrix )
first2 [ rows-except ] dip cols-except ;
: (m^n) ( m n -- n ) ALIAS: submatrix-excluding matrix-except
make-bits over first length identity-matrix
[ [ dupd m. ] when [ dup m. ] dip ] reduce nip ;
: m^n ( m n -- n ) :: matrix-except-all ( matrix -- submatrices )
dup 0 >= [ (m^n) ] [ negative-power-matrix ] if ; matrix dimension [ <iota> ] map first2-unsafe cartesian-product
[ [ matrix swap matrix-except ] map ] map ;
: stitch ( m -- m' ) ALIAS: all-submatrices matrix-except-all
[ ] [ [ append ] 2map ] map-reduce ;
: kron ( m1 m2 -- m ) : dimension ( matrix -- dimension )
'[ [ _ n*m ] map ] map stitch stitch ; [ { 0 0 } ]
[ [ length ] [ first-unsafe length ] bi 2array ] if-empty ;
: outer ( u v -- m )
[ n*v ] curry map ;
: row ( n matrix -- col )
nth ; inline
: rows ( seq matrix -- cols )
'[ _ row ] map ; inline
: col ( n matrix -- col )
swap '[ _ swap nth ] map ; inline
: cols ( seq matrix -- cols )
'[ _ col ] map ; inline
: set-index ( object pair matrix -- )
[ first2 swap ] dip nth set-nth ; inline
: set-indices ( object sequence matrix -- )
'[ _ set-index ] with each ; inline
: matrix-map ( matrix quot -- )
'[ _ map ] map ; inline
: column-map ( matrix quot -- seq )
[ [ first length <iota> ] keep ] dip '[ _ col @ ] map ; inline
: cartesian-square-indices ( n -- matrix )
<iota> dup cartesian-product ; inline
: cartesian-matrix-map ( matrix quot -- matrix' )
[ [ first length cartesian-square-indices ] keep ] dip
'[ _ @ ] matrix-map ; inline
: cartesian-matrix-column-map ( matrix quot -- matrix' )
[ cols first2 ] prepose cartesian-matrix-map ; inline
: cov-matrix-ddof ( matrix ddof -- cov )
'[ _ cov-ddof ] cartesian-matrix-column-map ; inline
: population-cov-matrix ( matrix -- cov ) 0 cov-matrix-ddof ; inline
: sample-cov-matrix ( matrix -- cov ) 1 cov-matrix-ddof ; inline
GENERIC: square-rows ( object -- matrix )
M: integer square-rows <iota> square-rows ;
M: sequence square-rows
[ length ] keep >array '[ _ clone ] { } replicate-as ;
GENERIC: square-cols ( object -- matrix )
M: integer square-cols <iota> square-cols ;
M: sequence square-cols
[ length ] keep [ <array> ] with { } map-as ;
: make-matrix-with-indices ( m n quot -- matrix )
[ [ <iota> ] bi@ ] dip cartesian-map ; inline
: null-matrix? ( matrix -- ? ) empty? ; inline
: well-formed-matrix? ( matrix -- ? )
[ t ] [
[ ] [ first length ] bi
'[ length _ = ] all?
] if-empty ;
: dim ( matrix -- pair/f )
[ 2 0 <array> ]
[ [ length ] [ first length ] bi 2array ] if-empty ;
: square-matrix? ( matrix -- ? )
{ [ well-formed-matrix? ] [ dim all-eq? ] } 1&& ;
: matrix-coordinates ( dim -- coordinates )
first2 [ <iota> ] bi@ cartesian-product ; inline
: dimension-range ( matrix -- dim range )
dim [ matrix-coordinates ] [ first [1,b] ] bi ;
: upper-matrix-indices ( matrix -- matrix' )
dimension-range <reversed> [ tail-slice* >array ] 2map concat ;
: lower-matrix-indices ( matrix -- matrix' )
dimension-range [ head-slice >array ] 2map concat ;
: make-lower-matrix ( object m n -- matrix )
zero-matrix [ lower-matrix-indices ] [ set-indices ] [ ] tri ;
: make-upper-matrix ( object m n -- matrix )
zero-matrix [ upper-matrix-indices ] [ set-indices ] [ ] tri ;

4
basis/math/vectors/vectors-docs.factor
View File

@ -315,11 +315,11 @@ HELP: vclamp
} ; } ;
HELP: v. HELP: v.
{ $values { "u" { $sequence real } } { "v" { $sequence real } } { "x" "a real number" } } { $values { "u" { $sequence real } } { "v" { $sequence real } } { "x" real } }
{ $description "Computes the dot product of two vectors." } ; { $description "Computes the dot product of two vectors." } ;
HELP: h. HELP: h.
{ $values { "u" { $sequence real } } { "v" { $sequence real } } { "x" "a real number" } } { $values { "u" { $sequence real } } { "v" { $sequence real } } { "x" real } }
{ $description "Computes the Hermitian inner product of two vectors." } ; { $description "Computes the Hermitian inner product of two vectors." } ;
HELP: vs+ HELP: vs+

8
basis/math/vectors/vectors-tests.factor
View File

@ -2,6 +2,7 @@ USING: math.vectors tools.test kernel specialized-arrays compiler
kernel.private alien.c-types math.functions ; kernel.private alien.c-types math.functions ;
SPECIALIZED-ARRAY: int SPECIALIZED-ARRAY: int
{ { 10 20 30 } } [ 10 { 1 2 3 } n*v ] unit-test
{ { 1 2 3 } } [ 1/2 { 2 4 6 } n*v ] unit-test { { 1 2 3 } } [ 1/2 { 2 4 6 } n*v ] unit-test
{ { 1 2 3 } } [ { 2 4 6 } 1/2 v*n ] unit-test { { 1 2 3 } } [ { 2 4 6 } 1/2 v*n ] unit-test
{ { 1 2 3 } } [ { 2 4 6 } 2 v/n ] unit-test { { 1 2 3 } } [ { 2 4 6 } 2 v/n ] unit-test
@ -49,3 +50,10 @@ SPECIALIZED-ARRAY: int
{ { 0 30 100 } } [ { { 0 30 100 } } [
{ -10 30 120 } { 0 0 0 } { 100 100 100 } vclamp { -10 30 120 } { 0 0 0 } { 100 100 100 } vclamp
] unit-test ] unit-test
{ { 0 0 1 } } [ { 1 0 0 } { 0 1 0 } cross ] unit-test
{ { 1 0 0 } } [ { 0 1 0 } { 0 0 1 } cross ] unit-test
{ { 0 1 0 } } [ { 0 0 1 } { 1 0 0 } cross ] unit-test
{ { 0.0 -0.707 0.707 } } [ { 1.0 0.0 0.0 } { 0.0 0.707 0.707 } cross ] unit-test
{ { 0 -2 2 } } [ { -1 -1 -1 } { 1 -1 -1 } cross ] unit-test
{ { 1 0 0 } } [ { 1 1 0 } { 1 0 0 } proj ] unit-test

16
basis/math/vectors/vectors.factor
View File

@ -279,3 +279,19 @@ PRIVATE>
: vclamp ( v min max -- w ) : vclamp ( v min max -- w )
rot vmin vmax ; inline rot vmin vmax ; inline
: cross ( vec1 vec2 -- vec3 )
[ [ { 1 2 0 } vshuffle ] [ { 2 0 1 } vshuffle ] bi* v* ]
[ [ { 2 0 1 } vshuffle ] [ { 1 2 0 } vshuffle ] bi* v* ] 2bi v- ; inline
:: normal ( vec1 vec2 vec3 -- vec4 )
vec2 vec1 v- vec3 vec1 v- cross normalize ; inline
: proj ( v u -- w )
[ [ v. ] [ norm-sq ] bi / ] keep n*v ;
: perp ( v u -- w )
dupd proj v- ;
: angle-between ( v u -- a )
[ normalize ] bi@ h. acos ;

12
extra/benchmark/3d-matrix-scalar/3d-matrix-scalar.factor
View File

@ -1,15 +1,15 @@
USING: kernel locals math math.matrices math.order math.vectors USING: kernel locals math math.matrices math.matrices.extras
prettyprint sequences ; math.order math.vectors prettyprint sequences ;
IN: benchmark.3d-matrix-scalar IN: benchmark.3d-matrix-scalar
:: p-matrix ( dim fov near far -- matrix ) :: p-matrix ( dim fov near far -- matrix )
dim dup first2 min v/n fov v*n near v*n dim dup first2 min v/n fov v*n near v*n
near far frustum-matrix4 ; near far <frustum-matrix4> ;
:: mv-matrix ( pitch yaw location -- matrix ) :: mv-matrix ( pitch yaw location -- matrix )
{ 1.0 0.0 0.0 } pitch rotation-matrix4 { 1.0 0.0 0.0 } pitch <rotation-matrix4>
{ 0.0 1.0 0.0 } yaw rotation-matrix4 { 0.0 1.0 0.0 } yaw <rotation-matrix4>
location vneg translation-matrix4 m. m. ; location vneg <translation-matrix4> m. m. ;
:: 3d-matrix-scalar-benchmark ( -- ) :: 3d-matrix-scalar-benchmark ( -- )
f :> result! f :> result!

5
extra/benchmark/3d-matrix-vector/3d-matrix-vector.factor
View File

@ -1,5 +1,6 @@
USING: kernel locals math math.matrices.simd math.order math.vectors USING: kernel locals math math.matrices math.matrices.simd
math.vectors.simd prettyprint sequences typed ; math.order math.vectors math.vectors.simd prettyprint sequences
typed ;
QUALIFIED-WITH: alien.c-types c QUALIFIED-WITH: alien.c-types c
IN: benchmark.3d-matrix-vector IN: benchmark.3d-matrix-vector

4
extra/benchmark/matrix-exponential-scalar/matrix-exponential-scalar.factor
View File

@ -1,4 +1,4 @@
USING: locals math math.combinatorics math.matrices USING: locals math math.combinatorics math.matrices math.matrices.extras
prettyprint sequences typed ; prettyprint sequences typed ;
IN: benchmark.matrix-exponential-scalar IN: benchmark.matrix-exponential-scalar
@ -15,7 +15,7 @@ IN: benchmark.matrix-exponential-scalar
:: matrix-exponential-scalar-benchmark ( -- ) :: matrix-exponential-scalar-benchmark ( -- )
f :> result! f :> result!
4 identity-matrix :> i4 4 <identity-matrix> :> i4
10000 [ 10000 [
i4 20 e^m result! i4 20 e^m result!
] times ] times

7
extra/game/debug/tests/tests.factor
View File

@ -2,8 +2,9 @@
! See http://factorcode.org/license.txt for BSD license. ! See http://factorcode.org/license.txt for BSD license.
USING: accessors colors.constants game.debug game.loop USING: accessors colors.constants game.debug game.loop
game.worlds gpu gpu.framebuffers gpu.util.wasd kernel literals game.worlds gpu gpu.framebuffers gpu.util.wasd kernel literals
locals make math math.matrices math.parser math.trig sequences locals make math math.matrices math.matrices.extras math.parser
specialized-arrays ui.gadgets.worlds ui.pixel-formats ; math.trig sequences specialized-arrays ui.gadgets.worlds
ui.pixel-formats ;
FROM: alien.c-types => float ; FROM: alien.c-types => float ;
SPECIALIZED-ARRAY: float SPECIALIZED-ARRAY: float
IN: game.debug.tests IN: game.debug.tests
@ -21,7 +22,7 @@ IN: game.debug.tests
{ 0 0 0 } { 1 0 0 } COLOR: red debug-line { 0 0 0 } { 1 0 0 } COLOR: red debug-line
{ 0 0 0 } { 0 1 0 } COLOR: green debug-line { 0 0 0 } { 0 1 0 } COLOR: green debug-line
{ 0 0 0 } { 0 0 1 } COLOR: blue debug-line { 0 0 0 } { 0 0 1 } COLOR: blue debug-line
{ -1.2 0 0 } { 0 1 0 } 0 deg>rad rotation-matrix3 debug-axes { -1.2 0 0 } { 0 1 0 } 0 deg>rad <rotation-matrix3> debug-axes
{ 3 5 -2 } { 3 2 1 } COLOR: white debug-box { 3 5 -2 } { 3 2 1 } COLOR: white debug-box
{ 0 9 0 } 8 2 COLOR: blue debug-cylinder { 0 9 0 } 8 2 COLOR: blue debug-cylinder
] float-array{ } make ] float-array{ } make

6
extra/gml/geometry/geometry.factor
View File

@ -1,8 +1,8 @@
! Copyright (C) 2010 Slava Pestov. ! Copyright (C) 2010 Slava Pestov.
USING: arrays kernel math.matrices math.vectors.simd.cords USING: arrays gml.runtime kernel math.matrices
math.trig gml.runtime ; math.matrices.extras math.trig math.vectors.simd.cords ;
IN: gml.geometry IN: gml.geometry
GML: rot_vec ( v n alpha -- v ) GML: rot_vec ( v n alpha -- v )
! Inefficient! ! Inefficient!
deg>rad rotation-matrix4 swap >array m.v >double-4 ; deg>rad <rotation-matrix4> swap >array m.v >double-4 ;

8
extra/gpu/demos/raytrace/raytrace.factor
View File

@ -2,9 +2,9 @@
! See http://factorcode.org/license.txt for BSD license. ! See http://factorcode.org/license.txt for BSD license.
USING: accessors arrays combinators.tuple game.loop game.worlds USING: accessors arrays combinators.tuple game.loop game.worlds
generalizations gpu gpu.render gpu.shaders gpu.util gpu.util.wasd generalizations gpu gpu.render gpu.shaders gpu.util gpu.util.wasd
kernel literals math math.libm math.matrices math.order math.vectors kernel literals math math.libm math.matrices math.matrices.extras
method-chains sequences ui ui.gadgets ui.gadgets.worlds math.order math.vectors method-chains sequences ui ui.gadgets
ui.pixel-formats audio.engine audio.loader locals ; ui.gadgets.worlds ui.pixel-formats audio.engine audio.loader locals ;
IN: gpu.demos.raytrace IN: gpu.demos.raytrace
GLSL-SHADER-FILE: raytrace-vertex-shader vertex-shader "raytrace.v.glsl" GLSL-SHADER-FILE: raytrace-vertex-shader vertex-shader "raytrace.v.glsl"
@ -48,7 +48,7 @@ TUPLE: raytrace-world < wasd-world
dup dtheta>> [ + ] curry change-theta drop ; dup dtheta>> [ + ] curry change-theta drop ;
: sphere-center ( sphere -- center ) : sphere-center ( sphere -- center )
[ [ axis>> ] [ theta>> ] bi rotation-matrix4 ] [ [ axis>> ] [ theta>> ] bi <rotation-matrix4> ]
[ home>> ] bi m.v ; [ home>> ] bi m.v ;
M: sphere audio-position sphere-center ; inline M: sphere audio-position sphere-center ; inline

18
extra/gpu/util/wasd/wasd.factor
View File

@ -4,8 +4,8 @@ USING: accessors arrays combinators.smart game.input
game.input.scancodes game.loop game.worlds game.input.scancodes game.loop game.worlds
gpu.render gpu.state kernel literals gpu.render gpu.state kernel literals
locals math math.constants math.functions math.matrices locals math math.constants math.functions math.matrices
math.order math.vectors opengl.gl sequences math.matrices.extras math.order math.vectors opengl.gl
ui ui.gadgets.worlds specialized-arrays audio.engine ; sequences ui ui.gadgets.worlds specialized-arrays audio.engine ;
FROM: alien.c-types => float ; FROM: alien.c-types => float ;
SPECIALIZED-ARRAY: float SPECIALIZED-ARRAY: float
IN: gpu.util.wasd IN: gpu.util.wasd
@ -38,14 +38,14 @@ GENERIC: wasd-fly-vertically? ( world -- ? )
M: wasd-world wasd-fly-vertically? drop t ; M: wasd-world wasd-fly-vertically? drop t ;
: wasd-mv-matrix ( world -- matrix ) : wasd-mv-matrix ( world -- matrix )
[ { 1.0 0.0 0.0 } swap pitch>> rotation-matrix4 ] [ { 1.0 0.0 0.0 } swap pitch>> <rotation-matrix4> ]
[ { 0.0 1.0 0.0 } swap yaw>> rotation-matrix4 ] [ { 0.0 1.0 0.0 } swap yaw>> <rotation-matrix4> ]
[ location>> vneg translation-matrix4 ] tri m. m. ; [ location>> vneg <translation-matrix4> ] tri m. m. ;
: wasd-mv-inv-matrix ( world -- matrix ) : wasd-mv-inv-matrix ( world -- matrix )
[ location>> translation-matrix4 ] [ location>> <translation-matrix4> ]
[ { 0.0 -1.0 0.0 } swap yaw>> rotation-matrix4 ] [ { 0.0 -1.0 0.0 } swap yaw>> <rotation-matrix4> ]
[ { -1.0 0.0 0.0 } swap pitch>> rotation-matrix4 ] tri m. m. ; [ { -1.0 0.0 0.0 } swap pitch>> <rotation-matrix4> ] tri m. m. ;
: wasd-p-matrix ( world -- matrix ) : wasd-p-matrix ( world -- matrix )
p-matrix>> ; p-matrix>> ;
@ -63,7 +63,7 @@ CONSTANT: fov 0.7
world wasd-far-plane :> far-plane world wasd-far-plane :> far-plane
world wasd-fov-vector near-plane v*n world wasd-fov-vector near-plane v*n
near-plane far-plane frustum-matrix4 ; near-plane far-plane <frustum-matrix4> ;
:: wasd-pixel-ray ( world loc -- direction ) :: wasd-pixel-ray ( world loc -- direction )
loc world dim>> [ /f 0.5 - 2.0 * ] 2map loc world dim>> [ /f 0.5 - 2.0 * ] 2map

0
basis/math/matrices/elimination/authors.txt → extra/math/matrices/elimination/authors.txt
View File

5
basis/math/matrices/elimination/elimination-docs.factor → extra/math/matrices/elimination/elimination-docs.factor
View File

@ -8,8 +8,9 @@ HELP: inverse
{ $examples { $examples
"A matrix multiplied by its inverse is the identity matrix." "A matrix multiplied by its inverse is the identity matrix."
{ $example { $example
"USING: kernel math.matrices math.matrices.elimination prettyprint ;" "USING: kernel math.matrices prettyprint ;"
"{ { 3 4 } { 7 9 } } dup inverse m. 2 identity-matrix = ." "FROM: math.matrices.elimination => inverse ;"
"{ { 3 4 } { 7 9 } } dup inverse m. 2 <identity-matrix> = ."
"t" "t"
} }
} ; } ;

0
basis/math/matrices/elimination/elimination-tests.factor → extra/math/matrices/elimination/elimination-tests.factor
View File

4
basis/math/matrices/elimination/elimination.factor → extra/math/matrices/elimination/elimination.factor
View File

@ -93,7 +93,7 @@ SYMBOL: matrix
: nullspace ( matrix -- seq ) : nullspace ( matrix -- seq )
echelon reduced dup empty? [ echelon reduced dup empty? [
dup first length identity-matrix [ dup first length <identity-matrix> [
[ [
dup leading drop dup leading drop
[ basis-vector ] [ drop ] if* [ basis-vector ] [ drop ] if*
@ -109,5 +109,5 @@ SYMBOL: matrix
: inverse ( matrix -- matrix ) ! Assumes an invertible matrix : inverse ( matrix -- matrix ) ! Assumes an invertible matrix
dup length dup length
[ identity-matrix [ append ] 2map solution ] keep [ <identity-matrix> [ append ] 2map solution ] keep
[ tail ] curry map ; [ tail ] curry map ;

0
basis/math/matrices/elimination/summary.txt → extra/math/matrices/elimination/summary.txt
View File

4
extra/math/matrices/extras/authors.txt Normal file
View File

@ -0,0 +1,4 @@
Slava Pestov
Joe Groff
Doug Coleman
Cat Stevens

641
extra/math/matrices/extras/extras-docs.factor Normal file
View File

@ -0,0 +1,641 @@
USING: arrays generic.single help.markup help.syntax kernel math
math.matrices math.matrices.private math.matrices.extras
math.order math.ratios math.vectors opengl.gl random sequences
urls ;
IN: math.matrices.extras
ABOUT: "math.matrices.extras"
ARTICLE: "math.matrices.extras" "Extra matrix operations"
"These constructions have special mathematical properties:"
{ $subsections
<box-matrix>
<hankel-matrix>
<hilbert-matrix>
<toeplitz-matrix>
<vandermonde-matrix>
}
"Common transformation matrices:"
{ $subsections
<frustum-matrix4>
<ortho-matrix4>
<rotation-matrix3>
<rotation-matrix4>
<scale-matrix3>
<scale-matrix4>
<skew-matrix4>
<translation-matrix4>
<random-integer-matrix>
<random-unit-matrix>
}
{ $subsections
invertible-matrix?
linearly-independent-matrix?
}
"Common algorithms on matrices:"
{ $subsections
gram-schmidt
gram-schmidt-normalize
kronecker-product
outer-product
}
"Matrix algebra:"
{ $subsections
mmin
mmax
mnorm
rank
nullity
} { $subsections
determinant 1/det m*1/det
>minors >cofactors
multiplicative-inverse
}
"Covariance in matrices:"
{ $subsections
covariance-matrix
covariance-matrix-ddof
sample-covariance-matrix
}
"Errors thrown by this vocabulary:"
{ $subsections negative-power-matrix non-square-determinant undefined-inverse } ;
HELP: invertible-matrix?
{ $values { "matrix" matrix } { "?" boolean } }
{ $description "Tests whether the input matrix has a " { $link multiplicative-inverse } ". In order for a matrix to be invertible, it must be a " { $link square-matrix } ", " { $emphasis "or" } ", if it is non-square, it must not be of " { $link +deficient-rank+ } "." }
{ $examples { $example "USING: math.matrices.extras prettyprint ;" "" } } ;
HELP: linearly-independent-matrix?
{ $values { "matrix" matrix } { "?" boolean } }
{ $description "Tests whether the input matrix is linearly independent." }
{ $examples { $example "USING: math.matrices.extras prettyprint ;" "" } } ;
! SINGLETON RANK TYPES
HELP: rank-kind
{ $class-description "The class of matrix rank quantifiers." } ;
HELP: +full-rank+
{ $class-description "A " { $link rank-kind } " describing a matrix of full rank." } ;
HELP: +half-rank+
{ $class-description "A " { $link rank-kind } " describing a matrix of half rank." } ;
HELP: +zero-rank+
{ $class-description "A " { $link rank-kind } " describing a matrix of zero rank." } ;
HELP: +deficient-rank+
{ $class-description "A " { $link rank-kind } " describing a matrix of deficient rank." } ;
HELP: +uncalculated-rank+
{ $class-description "A " { $link rank-kind } " describing a matrix whose rank is not (yet) known." } ;
! ERRORS
HELP: negative-power-matrix
{ $values { "m" matrix } { "n" integer } }
{ $description "Throws a " { $link negative-power-matrix } " error." }
{ $error-description "Given the semantics of " { $link m^n } ", negative exponents are not within the domain of the power matrix function." } ;
HELP: non-square-determinant
{ $values { "m" integer } { "n" integer } }
{ $description "Throws a " { $link non-square-determinant } " error." }
{ $error-description { $link determinant } " was used with a non-square matrix whose dimensions are " { $snippet "m x n" } ". It is not generally possible to find the determinant of a non-square matrix." } ;
HELP: undefined-inverse
{ $values { "m" integer } { "n" integer } { "r" rank-kind } }
{ $description "Throws an " { $link undefined-inverse } " error." }
{ $error-description { $link multiplicative-inverse } " was used with a non-square matrix of rank " { $snippet "rank" } " whose dimensions are " { $snippet "m x n" } ". It is not generally possible to find the inverse of a " { $link +deficient-rank+ } " non-square " { $link matrix } "." } ;
HELP: <random-integer-matrix>
{ $values { "m" integer } { "n" integer } { "max" integer } { "matrix" matrix } }
{ $description "Creates a " { $snippet "m x n" } " " { $link matrix } " full of random, possibly signed " { $link integer } "s whose absolute values are less than or equal to " { $snippet "max" } ", as given by " { $link random-integers } "." }
{ $notelist
{ "The signedness of the numbers in the resulting matrix will be randomized. Use " { $link mabs } " with this word to generate a matrix of random positive integers." }
{ $equiv-word-note "integral" <random-unit-matrix> }
}
{ $errors { $link no-method } " if " { $snippet "max"} " is not an " { $link integer } "." }
{ $examples
{ $unchecked-example
"USING: math.matrices.extras prettyprint ;"
"2 4 15 <random-integer-matrix> ."
"{ { -9 -9 1 3 } { -14 -8 14 10 } }"
}
} ;
HELP: <random-unit-matrix>
{ $values { "m" integer } { "n" integer } { "max" number } { "matrix" matrix } }
{ $description "Creates a " { $snippet "m x n" } " " { $link matrix } " full of random, possibly signed " { $link float } "s as a fraction of " { $snippet "max" } "." }
{ $notelist
{ "The signedness of the numbers in the resulting matrix will be randomized. Use " { $link mabs } " with this word to generate a matrix of random positive numbers." }
{ $equiv-word-note "real" <random-integer-matrix> }
{ "This word is implemented by generating sub-integral floats through " { $link random-units } " and multiplying by random integers less than or equal to " { $snippet "max" } "." }
}
{ $examples
{ $unchecked-example
"USING: math.matrices.extras prettyprint ;"
"4 2 15 <random-unit-matrix> ."
"{
{ -3.713295909201797 3.815787135075961 }
{ -2.460506890603817 1.535222788710546 }
{ 3.692213981267878 -1.462963244399762 }
{ 13.8967592095433 -6.688509969360172 }
}"
}
} ;
HELP: <hankel-matrix>
{ $values { "n" integer } { "matrix" matrix } }
{ $description
"A Hankel matrix is a symmetric, " { $link square-matrix } " in which each ascending skew-diagonal from left to right is constant. See " { $url URL" https://en.wikipedia.org/wiki/Hankel_matrix" "hankel matrix" } "."
$nl
"The following is true of any Hankel matrix" { $snippet "A" } ": " { $snippet "A[i][j] = A[j][i] = a[i+j-2]" } "."
$nl
"The " { $link <toeplitz-matrix> } " is an upside-down Hankel matrix."
$nl
"The " { $link <hilbert-matrix> } " is a special case of the Hankel matrix."
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"4 <hankel-matrix> ."
"{ { 1 2 3 4 } { 2 3 4 0 } { 3 4 0 0 } { 4 0 0 0 } }"
}
} ;
HELP: <hilbert-matrix>
{ $values { "m" integer } { "n" integer } { "matrix" matrix } }
{ $description
"A Hilbert matrix is a " { $link square-matrix } " " { $snippet "A" } " in which entries are the unit fractions "
{ $snippet "A[i][j] = 1/(i+j-1)" }
". See " { $url URL" https://en.wikipedia.org/wiki/Hilbert_matrix" "hilbert matrix" } "."
$nl
"A Hilbert matrix is a special case of the " { $link <hankel-matrix> } "."
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"1 2 <hilbert-matrix> ."
"{ { 1 1/2 } }"
}
{ $example
"USING: math.matrices.extras prettyprint ;"
"3 6 <hilbert-matrix> ."
"{
{ 1 1/2 1/3 1/4 1/5 1/6 }
{ 1/2 1/3 1/4 1/5 1/6 1/7 }
{ 1/3 1/4 1/5 1/6 1/7 1/8 }
}"
}
} ;
HELP: <toeplitz-matrix>
{ $values { "n" integer } { "matrix" matrix } }
{ $description "A Toeplitz matrix is an upside-down " { $link <hankel-matrix> } ". Unlike the Hankel matrix, a Toeplitz matrix can be non-square. See " { $url URL" https://en.wikipedia.org/wiki/Hankel_matrix" "hankel matrix" } "."
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"4 <toeplitz-matrix> ."
"{ { 1 2 3 4 } { 2 1 2 3 } { 3 2 1 2 } { 4 3 2 1 } }"
}
} ;
HELP: <box-matrix>
{ $values { "r" integer } { "matrix" matrix } }
{ $description "Create a box matrix (a " { $link square-matrix } ") with the dimensions of " { $snippet "r x r" } ", filled with ones. The number of elements in the output scales linearly (" { $snippet "(r*2)+1" } ") with " { $snippet "r" } "." }
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"2 <box-matrix> ."
"{
{ 1 1 1 1 1 }
{ 1 1 1 1 1 }
{ 1 1 1 1 1 }
{ 1 1 1 1 1 }
{ 1 1 1 1 1 }
}"
}
{ $example
"USING: math.matrices.extras prettyprint ;"
"3 <box-matrix> ."
"{
{ 1 1 1 1 1 1 1 }
{ 1 1 1 1 1 1 1 }
{ 1 1 1 1 1 1 1 }
{ 1 1 1 1 1 1 1 }
{ 1 1 1 1 1 1 1 }
{ 1 1 1 1 1 1 1 }
{ 1 1 1 1 1 1 1 }
}"
}
} ;
HELP: <scale-matrix3>
{ $values { "factors" sequence } { "matrix" matrix } }
{ $description "Make a " { $snippet "3 x 3" } " scaling matrix, used to scale an object in 3 dimensions. See " { $url URL" https://en.wikipedia.org/wiki/Scaling_(geometry)#Matrix_representation" "scaling matrix on Wikipedia" } "." }
{ $notelist
{ $finite-input-note "three" "factors" }
{ $equiv-word-note "3-matrix" <scale-matrix4> }
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{ 22 33 -44 } <scale-matrix4> ."
"{
{ 22 0.0 0.0 0.0 }
{ 0.0 33 0.0 0.0 }
{ 0.0 0.0 -44 0.0 }
{ 0.0 0.0 0.0 1.0 }
}"
}
} ;
HELP: <scale-matrix4>
{ $values { "factors" sequence } { "matrix" matrix } }
{ $description "Make a " { $snippet "4 x 4" } " scaling matrix, used to scale an object in 3 or more dimensions. See " { $url URL" https://en.wikipedia.org/wiki/Scaling_(geometry)#Matrix_representation" "scaling matrix on Wikipedia" } "." }
{ $notelist
{ $finite-input-note "three" "factors" }
{ $equiv-word-note "4-matrix" <scale-matrix3> }
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{ 22 33 -44 } <scale-matrix4> ."
"{
{ 22 0.0 0.0 0.0 }
{ 0.0 33 0.0 0.0 }
{ 0.0 0.0 -44 0.0 }
{ 0.0 0.0 0.0 1.0 }
}"
}
} ;
HELP: <ortho-matrix4>
{ $values { "factors" sequence } { "matrix" matrix } }
{ $description "Create a " { $link <scale-matrix4> } ", with the scale factors inverted." }
{ $notelist
{ $finite-input-note "three" "factors" }
{ $equiv-word-note "inverse" <scale-matrix4> }
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{ -9.3 100 1/2 } <ortho-matrix4> ."
"{
{ -0.1075268817204301 0.0 0.0 0.0 }
{ 0.0 1/100 0.0 0.0 }
{ 0.0 0.0 2 0.0 }
{ 0.0 0.0 0.0 1.0 }
}"
}
} ;
HELP: <frustum-matrix4>
{ $values { "xy-dim" pair } { "near" number } { "far" number } { "matrix" matrix } }
{ $description "Make a " { $snippet "4 x 4" } " matrix suitable for representing an occlusion frustum. A viewing or occlusion frustum is the three-dimensional region of a three-dimensional object which is visible on the screen. See " { $url URL" https://en.wikipedia.org/wiki/Frustum" "frustum on Wikipedia" } "." }
{ $notes { $finite-input-note "two" "xy-dim" } }
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{ 5 4 } 5 6 <frustum-matrix4> ."
"{
{ 1.0 0.0 0.0 0.0 }
{ 0.0 1.25 0.0 0.0 }
{ 0.0 0.0 -11.0 -60.0 }
{ 0.0 0.0 -1.0 0.0 }
}"
}
} ;
{ <frustum-matrix4> glFrustum } related-words
HELP: cartesian-matrix-map
{ $values { "matrix" matrix } { "quot" { $quotation ( ... pair matrix -- ... matrix' ) } } { "matrix-seq" { $sequence matrix } } }
{ $description "Calls the quotation with the matrix and the coordinate pair of the current element on the stack, with the matrix on the top of the stack." }
{ $examples
{ $example
"USING: arrays math.matrices.extras prettyprint ;"
"{ { 21 22 } { 23 24 } } [ 2array ] cartesian-matrix-map ."
"{
{
{ { 0 0 } { { 21 22 } { 23 24 } } }
{ { 0 1 } { { 21 22 } { 23 24 } } }
}
{
{ { 1 0 } { { 21 22 } { 23 24 } } }
{ { 1 1 } { { 21 22 } { 23 24 } } }
}
}"
}
}
{ $notelist
{ $equiv-word-note "orthogonal" cartesian-column-map }
{ $equiv-word-note "two-dimensional" map-index }
$2d-only-note
} ;
HELP: cartesian-column-map
{ $values { "matrix" matrix } { "quot" { $quotation ( ... pair matrix -- ... matrix' ) } } { "matrix-seq" { $sequence matrix } } }
{ $notelist
{ $equiv-word-note "orthogonal" cartesian-matrix-map }
$2d-only-note
} ;
HELP: gram-schmidt
{ $values { "matrix" matrix } { "orthogonal" matrix } }
{ $description "Apply a Gram-Schmidt transform on the matrix." }
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{ { 1 2 } { 3 4 } { 5 6 } } gram-schmidt ."
"{ { 1 2 } { 4/5 -2/5 } { 0 0 } }"
}
} ;
HELP: gram-schmidt-normalize
{ $values { "matrix" matrix } { "orthonormal" matrix } }
{ $description "Apply a Gram-Schmidt transform on the matrix, and " { $link normalize } " each row of the result, resulting in an orthogonal and normalized matrix (orthonormal)." }
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{ { 1 2 } { 3 4 } { 5 6 } } gram-schmidt-normalize ."
"{
{ 0.4472135954999579 0.8944271909999159 }
{ 0.894427190999916 -0.447213595499958 }
{ NAN: 8000000000000 NAN: 8000000000000 }
}"
}
} ;
HELP: m^n
{ $values { "m" matrix } { "n" object } }
{ $description "Compute the " { $snippet "nth" } " power of the input matrix. If " { $snippet "n" } " is " { $snippet "-1" } ", the inverse of the matrix is calculated (but see " { $link multiplicative-inverse } " for pitfalls)." }
{ $errors
{ $link negative-power-matrix } " if " { $snippet "n" } " is a negative number other than " { $snippet "-1" } "."
$nl
{ $link undefined-inverse } " if " { $snippet "n" } " is " { $snippet "-1" } " and the " { $link multiplicative-inverse } " of " { $snippet "m" } " is undefined."
}
{ $notelist
{ $equiv-word-note "swapped" n^m }
$2d-only-note
{ $matrix-scalar-note max abs / }
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{ { 1 2 } { 3 4 } } 2 m^n ."
"{ { 7 10 } { 15 22 } }"
}
} ;
HELP: n^m
{ $values { "n" object } { "m" matrix } }
{ $description "Because it is nonsensical to raise a number to the power of a matrix, this word exists to save typing " { $snippet "swap m^n" } ". See " { $link m^n } " for more information." }
{ $errors
{ $link negative-power-matrix } " if " { $snippet "n" } " is a negative number other than " { $snippet "-1" } "."
$nl
{ $link undefined-inverse } " if " { $snippet "n" } " is " { $snippet "-1" } " and the " { $link multiplicative-inverse } " of " { $snippet "m" } " is undefined."
}
{ $notelist
{ $equiv-word-note "swapped" m^n }
$2d-only-note
{ $matrix-scalar-note max abs / }
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"2 { { 1 2 } { 3 4 } } n^m ."
"{ { 7 10 } { 15 22 } }"
}
} ;
HELP: kronecker-product
{ $values { "m1" matrix } { "m2" matrix } { "m" matrix } }
{ $description "Calculates the " { $url URL" http://enwp.org/Kronecker_product" "Kronecker product" } " of two matrices. This product can be described as a generalization of the vector-based " { $link outer-product } " to matrices. The Kronecker product gives the matrix of the tensor product with respect to a standard choice of basis." }
{ $notelist
{ $equiv-word-note "matrix" outer-product }
$2d-only-note
{ $matrix-scalar-note * }
}
{ $examples
{ $unchecked-example
"USING: math.matrices.extras prettyprint ;"
"{
{ 1 2 }
{ 3 4 }
} {
{ 0 5 }
{ 6 7 }
} kronecker-product ."
"{
{ 0 5 0 10 }
{ 6 7 12 14 }
{ 0 15 0 20 }
{ 18 21 24 28 }
}" }
} ;
HELP: outer-product
{ $values { "u" sequence } { "v" sequence } { "matrix" matrix } }
{ $description "Computes the " { $url URL" http:// enwp.org/Outer_product" "outer-product product" } " of " { $snippet "u" } " and " { $snippet "v" } "." }
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{ 5 6 7 } { 1 2 3 } outer-product ."
"{ { 5 10 15 } { 6 12 18 } { 7 14 21 } }" }
} ;
HELP: rank
{ $values { "matrix" matrix } { "rank" rank-kind } }
{ $contract "The " { $emphasis "rank" } " of a " { $link matrix } " is how its number of linearly independent columns compare to the maximal number of linearly independent columns for a matrix with the same dimension." }
{ $notes "See " { $url "https://en.wikipedia.org/wiki/Rank_(linear_algebra)" } " for more information." } ;
HELP: nullity
{ $values { "matrix" matrix } { "nullity" rank-kind } }
;
HELP: determinant
{ $values { "matrix" square-matrix } { "determinant" number } }
{ $contract "Compute the determinant of the input matrix. Generally, the determinant of a matrix is a scaling factor of the transformation described by the matrix." }
{ $notelist
$2d-only-note
{ $matrix-scalar-note max - * }
}
{ $errors { $link non-square-determinant } " if the input matrix is not a " { $link square-matrix } "." }
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{
{ 3 0 -1 }
{ -3 1 3 }
{ 2 -5 4 }
} determinant ."
"44"
}
{ $example
"USING: math.matrices.extras prettyprint ;"
"{
{ -8 -8 13 11 10 -5 -14 }
{ 3 -11 -8 3 -7 -3 4 }
{ 10 4 -5 3 0 -6 -12 }
{ -14 0 -3 -8 10 0 10 }
{ 3 -6 1 -10 -9 10 0 }
{ 5 -12 -14 6 5 -1 -7 }
{ -9 -14 -8 5 2 2 -2 }
} determinant ."
"-103488155"
}
} ;
HELP: 1/det
{ $values { "matrix" square-matrix } { "1/det" number } }
{ $description "Find the inverse (" { $link recip } ") of the " { $link determinant } " of the input matrix." }
{ $notelist
$2d-only-note
{ $matrix-scalar-note determinant recip }
}
{ $errors
{ $link non-square-determinant } " if the input matrix is not a " { $link square-matrix } "."
$nl
{ $link division-by-zero } " if the " { $link determinant } " of the input matrix is " { $snippet "0" } "."
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{
{ 0 10 -12 4 }
{ -9 6 -11 9 }
{ -5 -10 0 2 }
{ -7 -11 10 11 }
} 1/det ."
"-1/9086"
}
} ;
HELP: m*1/det
{ $values { "matrix" square-matrix } { "matrix'" square-matrix } }
{ $description "Multiply the input matrix by the inverse (" { $link recip } ") of its " { $link determinant } "." }
{ $notelist
{ "This word is used to implement " { $link recip } " for " { $link square-matrix } "." }
$2d-only-note
{ $matrix-scalar-note determinant recip }
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{
{ -14 0 -13 7 }
{ -4 11 7 -12 }
{ -3 2 9 -14 }
{ 3 -5 10 -2 }
} m*1/det ."
"{
{ 7/6855 0 13/13710 -7/13710 }
{ 2/6855 -11/13710 -7/13710 2/2285 }
{ 1/4570 -1/6855 -3/4570 7/6855 }
{ -1/4570 1/2742 -1/1371 1/6855 }
}"
}
}
;
HELP: >minors
{ $values { "matrix" square-matrix } { "matrix'" square-matrix } }
{ $description "Calculate the " { $emphasis "matrix of minors" } " of the input matrix. See " { $url URL" https://en.wikipedia.org/wiki/Minor_(linear_algebra)" "minor on Wikipedia" } "." }
{ $notelist
$keep-shape-note
$2d-only-note
{ $matrix-scalar-note determinant }
}
{ $errors { $link non-square-determinant } " if the input matrix is not a " { $link square-matrix } "." }
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{
{ -8 0 7 -11 }
{ 15 0 -3 -11 }
{ 1 -10 -4 6 }
{ 11 -15 3 -15 }
} >minors ."
"{
{ 1710 -130 2555 -1635 }
{ -690 -286 -2965 1385 }
{ 1650 -754 3795 -1215 }
{ 1100 416 2530 -810 }
}"
}
} ;
HELP: >cofactors
{ $values { "matrix" matrix } { "matrix'" matrix } }
{ $description "Calculate the " { $emphasis "matrix of cofactors" } " of the input matrix. See " { $url URL" https://en.wikipedia.org/wiki/Minor_(linear_algebra)#Inverse_of_a_matrix" "matrix of cofactors on Wikipedia" } ". Alternating elements of the input matrix have their signs inverted." $nl "On odd rows, the even elements have their signs inverted. On even rows, odd elements have their signs inverted." }
{ $notelist
$keep-shape-note
$2d-only-note
{ $matrix-scalar-note neg }
}
{ $examples
{ $example
"USING: math.matrices.extras prettyprint ;"
"{
{ 8 0 7 11 }
{ 15 0 3 11 }
{ 1 10 4 6 }
{ 11 15 3 15 }
} >cofactors ."
"{
{ 8 0 7 -11 }
{ -15 0 -3 11 }
{ 1 -10 4 -6 }
{ -11 15 -3 15 }
}"
}
{ $example
"USING: math.matrices.extras prettyprint ;"
"{
{ -8 0 7 -11 }
{ 15 0 -3 -11 }
{ 1 -10 -4 6 }
{ 11 -15 3 -15 }
} >cofactors ."
"{
{ -8 0 7 11 }
{ -15 0 3 -11 }
{ 1 10 -4 -6 }
{ -11 -15 -3 -15 }
}"
}
} ;
HELP: multiplicative-inverse
{ $values { "x" matrix } { "y" matrix } }
{ $description "Calculate the multiplicative inverse of the input." $nl "If the input is a " { $link square-matrix } ", this is done by multiplying the " { $link transpose } " of the " { $link2 >cofactors "cofactors" } " of the " { $link2 >minors "minors" } " of the input matrix by the " { $link2 1/det "inverse of the determinant" } " of the input matrix." }
{ $notelist
$keep-shape-note
$2d-only-note
{ $matrix-scalar-note determinant >cofactors 1/det }
}
{ $errors { $link non-square-determinant } " if the input matrix is not a " { $link square-matrix } "." } ;
HELP: covariance-matrix-ddof
{ $values { "matrix" matrix } { "ddof" object } { "cov" matrix } }
;
HELP: covariance-matrix
{ $values { "matrix" matrix } { "cov" matrix } }
;
HELP: sample-covariance-matrix
{ $values { "matrix" matrix } { "cov" matrix } }
;
HELP: population-covariance-matrix
{ $values { "matrix" matrix } { "cov" matrix } }
;

362
extra/math/matrices/extras/extras-tests.factor Normal file
View File

@ -0,0 +1,362 @@
! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens.
USING: arrays combinators.short-circuit grouping kernel math
math.matrices math.matrices.extras math.matrices.extras.private
math.statistics math.vectors sequences sequences.deep sets
tools.test ;
IN: math.matrices.extras
<PRIVATE
: call-eq? ( obj quots -- ? )
[ call( x -- x ) ] with map all-eq? ; ! inline
PRIVATE>
{ {
{ 4181 6765 }
{ 6765 10946 }
} } [
{ { 0 1 } { 1 1 } } 20 m^n
] unit-test
[ { { 0 1 } { 1 1 } } -20 m^n ] [ negative-power-matrix? ] must-fail-with
[ { { 0 1 } { 1 1 } } -8 m^n ] [ negative-power-matrix? ] must-fail-with
{ { 1 -2 3 -4 } } [ { 1 2 3 4 } t alternating-sign ] unit-test
{ { -1 2 -3 4 } } [ { 1 2 3 4 } f alternating-sign ] unit-test
{ t } [ 50 <box-matrix> dup transpose = ] unit-test
{ t } [ 50 <box-matrix> dup anti-transpose = ] unit-test
{ f } [ 4 <box-matrix> zero-matrix? ] unit-test
{ t } [ 2 4 15 <random-integer-matrix> mabs {
[ flatten [ 15 <= ] all? ]
[ regular-matrix? ]
[ length 2 = ]
[ first length 4 = ]
} 1&& ] unit-test
{ t } [ 4 4 -45 <random-integer-matrix> mabs {
[ flatten [ 45 <= ] all? ]
[ regular-matrix? ]
[ length 4 = ]
[ first length 4 = ]
} 1&& ] unit-test
{ t } [ 2 2 1 <random-integer-matrix> mabs {
[ flatten [ 1 <= ] all? ]
[ regular-matrix? ]
[ length 2 = ]
[ first length 2 = ]
} 1&& ] unit-test
{ t } [ 2 4 .89 <random-unit-matrix> mabs {
[ flatten [ .89 <= ] all? ]
[ regular-matrix? ]
[ length 2 = ]
[ first length 4 = ]
} 1&& ] unit-test
{ t } [ 2 4 -45.89 <random-unit-matrix> mabs {
[ flatten [ 45.89 <= ] all? ]
[ regular-matrix? ]
[ length 2 = ]
[ first length 4 = ]
} 1&& ] unit-test
{ t } [ 4 4 .89 <random-unit-matrix> mabs {
[ flatten [ .89 <= ] all? ]
[ regular-matrix? ]
[ length 4 = ]
[ first length 4 = ]
} 1&& ] unit-test
{ {
{ 1 1/2 1/3 1/4 }
{ 1/2 1/3 1/4 1/5 }
{ 1/3 1/4 1/5 1/6 }
} } [ 3 4 <hilbert-matrix> ] unit-test
{ {
{ 1 2 3 4 }
{ 2 1 2 3 }
{ 3 2 1 2 }
{ 4 3 2 1 }
} } [ 4 <toeplitz-matrix> ] unit-test
{ {
{ 1 2 3 4 }
{ 2 3 4 0 }
{ 3 4 0 0 }
{ 4 0 0 0 } }
} [ 4 <hankel-matrix> ] unit-test
{ {
{ 1 1 1 }
{ 4 2 1 }
{ 9 3 1 }
{ 25 5 1 } }
} [
{ 1 2 3 5 } 3 <vandermonde-matrix>
] unit-test
{ {
{ 0 5 0 10 }
{ 6 7 12 14 }
{ 0 15 0 20 }
{ 18 21 24 28 }
} } [ {
{ 1 2 }
{ 3 4 }
} {
{ 0 5 }
{ 6 7 }
} kronecker-product ] unit-test
{ {
{ 1 1 1 1 }
{ 1 -1 1 -1 }
{ 1 1 -1 -1 }
{ 1 -1 -1 1 }
} } [ {
{ 1 1 }
{ 1 -1 }
} dup kronecker-product ] unit-test
{ {
{ 1 1 1 1 1 1 1 1 }
{ 1 -1 1 -1 1 -1 1 -1 }
{ 1 1 -1 -1 1 1 -1 -1 }
{ 1 -1 -1 1 1 -1 -1 1 }
{ 1 1 1 1 -1 -1 -1 -1 }
{ 1 -1 1 -1 -1 1 -1 1 }
{ 1 1 -1 -1 -1 -1 1 1 }
{ 1 -1 -1 1 -1 1 1 -1 }
} } [ {
{ 1 1 }
{ 1 -1 }
} dup dup kronecker-product kronecker-product ] unit-test
{ {
{ 1 1 1 1 1 1 1 1 }
{ 1 -1 1 -1 1 -1 1 -1 }
{ 1 1 -1 -1 1 1 -1 -1 }
{ 1 -1 -1 1 1 -1 -1 1 }
{ 1 1 1 1 -1 -1 -1 -1 }
{ 1 -1 1 -1 -1 1 -1 1 }
{ 1 1 -1 -1 -1 -1 1 1 }
{ 1 -1 -1 1 -1 1 1 -1 }
} } [ {
{ 1 1 }
{ 1 -1 }
} dup dup kronecker-product swap kronecker-product ] unit-test
! kronecker-product is not generally commutative, make sure we have the right order
{ {
{ 1 2 3 4 5 1 2 3 4 5 }
{ 6 7 8 9 10 6 7 8 9 10 }
{ 1 2 3 4 5 -1 -2 -3 -4 -5 }
{ 6 7 8 9 10 -6 -7 -8 -9 -10 }
} } [ {
{ 1 1 }
{ 1 -1 }
} {
{ 1 2 3 4 5 }
{ 6 7 8 9 10 }
} kronecker-product ] unit-test
{ {
{ 1 1 2 2 3 3 4 4 5 5 }
{ 1 -1 2 -2 3 -3 4 -4 5 -5 }
{ 6 6 7 7 8 8 9 9 10 10 }
{ 6 -6 7 -7 8 -8 9 -9 10 -10 }
} } [ {
{ 1 1 }
{ 1 -1 }
} {
{ 1 2 3 4 5 }
{ 6 7 8 9 10 }
} swap kronecker-product ] unit-test
{ {
{ 5 10 15 }
{ 6 12 18 }
{ 7 14 21 }
} } [
{ 5 6 7 }
{ 1 2 3 }
outer-product
] unit-test
CONSTANT: test-points {
{ 80 27 89 } { 80 27 88 } { 75 25 90 }
{ 62 24 87 } { 62 22 87 } { 62 23 87 }
{ 62 24 93 } { 62 24 93 } { 58 23 87 }
{ 58 18 80 } { 58 18 89 } { 58 17 88 }
{ 58 18 82 } { 58 19 93 } { 50 18 89 }
{ 50 18 86 } { 50 19 72 } { 50 19 79 }
{ 50 20 80 } { 56 20 82 } { 70 20 91 }
}
{ {
{ 84+2/35 22+23/35 24+4/7 }
{ 22+23/35 9+104/105 6+87/140 }
{ 24+4/7 6+87/140 28+5/7 }
} } [ test-points sample-covariance-matrix ] unit-test
{ {
{ 80+8/147 21+85/147 23+59/147 }
{ 21+85/147 9+227/441 6+15/49 }
{ 23+59/147 6+15/49 27+17/49 }
} } [ test-points covariance-matrix ] unit-test
{
{
{ 80+8/147 21+85/147 23+59/147 }
{ 21+85/147 9+227/441 6+15/49 }
{ 23+59/147 6+15/49 27+17/49 }
}
} [
test-points population-covariance-matrix
] unit-test
{ t } [ { { 1 } }
{ [ drop 1 ] [ (1determinant) ] [ 1 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ 0 } [ { { 0 } } determinant ] unit-test
{ t } [ {
{ 4 6 } ! order is significant
{ 3 8 }
} { [ drop 14 ] [ (2determinant) ] [ 2 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 3 8 }
{ 4 6 }
} { [ drop -14 ] [ (2determinant) ] [ 2 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 2 5 }
{ 1 -3 }
} { [ drop -11 ] [ (2determinant) ] [ 2 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 1 -3 }
{ 2 5 }
} { [ drop 11 ] [ (2determinant) ] [ 2 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 3 0 -1 }
{ 2 -5 4 }
{ -3 1 3 }
} { [ drop -44 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 3 0 -1 }
{ -3 1 3 }
{ 2 -5 4 }
} { [ drop 44 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 2 -3 1 }
{ 4 2 -1 }
{ -5 3 -2 }
} { [ drop -19 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 2 -3 1 }
{ -5 3 -2 }
{ 4 2 -1 }
} { [ drop 19 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 4 2 -1 }
{ 2 -3 1 }
{ -5 3 -2 }
} { [ drop 19 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 5 1 -2 }
{ -1 0 4 }
{ 2 -3 3 }
} { [ drop 65 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 6 1 1 }
{ 4 -2 5 }
{ 2 8 7 }
} { [ drop -306 ] [ (3determinant) ] [ 3 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ -5 4 -3 2 }
{ -2 1 0 -1 }
{ -2 -3 -4 -5 }
{ 0 2 0 4 }
} { [ drop -24 ] [ 4 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ t } [ {
{ 2 4 2 2 }
{ 5 1 -6 10 }
{ 4 3 -1 7 }
{ 9 8 7 3 }
} { [ drop 272 ] [ 4 swap (ndeterminant) ] [ determinant ] }
call-eq?
] unit-test
{ {
{ 2 2 2 }
{ -2 3 3 }
{ 0 -10 0 }
} } [ {
{ 3 0 2 }
{ 2 0 -2 }
{ 0 1 1 }
} >minors ] unit-test
! i think this unit test is wrong
! { {
! { 1 -6 -13 }
! { 0 0 0 }
! { 1 -6 -13 }
! } } [ {
! { 1 2 1 }
! { 6 -1 0 }
! { 1 -2 -1 }
! } >minors ] unit-test
{ {
{ 1 6 -13 }
{ 0 0 0 }
{ 1 6 -13 }
} } [ {
{ 1 -6 -13 }
{ 0 0 0 }
{ 1 -6 -13 }
} >cofactors ] unit-test

338
extra/math/matrices/extras/extras.factor Normal file
View File

@ -0,0 +1,338 @@
! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens.
! See http://factorcode.org/license.txt for BSD license.
USING: accessors arrays combinators formatting fry kernel locals
math math.bits math.functions math.matrices
math.matrices.private math.order math.statistics math.text.english
math.vectors random sequences sequences.deep summary ;
IN: math.matrices.extras
! this is a questionable implementation
SINGLETONS: +full-rank+ +half-rank+ +zero-rank+ +deficient-rank+ +uncalculated-rank+ ;
UNION: rank-kind +full-rank+ +half-rank+ +zero-rank+ +deficient-rank+ +uncalculated-rank+ ;
ERROR: negative-power-matrix
{ m matrix } { n integer } ;
ERROR: non-square-determinant
{ m integer } { n integer } ;
ERROR: undefined-inverse
{ m integer } { n integer } { r rank-kind initial: +uncalculated-rank+ } ;
M: negative-power-matrix summary
n>> dup ordinal-suffix "%s%s power of a matrix is undefined" sprintf ;
M: non-square-determinant summary
[ m>> ] [ n>> ] bi "non-square %s x %s matrix has no determinant" sprintf ;
M: undefined-inverse summary
[ m>> ] [ n>> ] [ r>> name>> ] tri "%s x %s matrix of rank %s has no inverse" sprintf ;
<PRIVATE
DEFER: alternating-sign
: finish-randomizing-matrix ( matrix -- matrix' )
[ f alternating-sign randomize ] map randomize ; inline
PRIVATE>
: <random-integer-matrix> ( m n max -- matrix )
'[ _ _ 1 + random-integers ] replicate
finish-randomizing-matrix ; inline
: <random-unit-matrix> ( m n max -- matrix )
'[ _ random-units [ _ * ] map ] replicate
finish-randomizing-matrix ; inline
<PRIVATE
: (gram-schmidt) ( v seq -- newseq )
[ dupd proj v- ] each ;
PRIVATE>
: gram-schmidt ( matrix -- orthogonal )
[ V{ } clone [ over (gram-schmidt) suffix! ] reduce ] keep like ;
: gram-schmidt-normalize ( matrix -- orthonormal )
gram-schmidt [ normalize ] map ; inline
: kronecker-product ( m1 m2 -- m )
'[ [ _ n*m ] map ] map stitch stitch ;
: outer-product ( u v -- matrix )
'[ _ n*v ] map ;
! Special matrix constructors follow
: <hankel-matrix> ( n -- matrix )
[ <iota> dup ] keep '[ + abs 1 + dup _ > [ drop 0 ] when ] cartesian-map ;
: <hilbert-matrix> ( m n -- matrix )
[ <iota> ] bi@ [ + 1 + recip ] cartesian-map ;
: <toeplitz-matrix> ( n -- matrix )
<iota> dup [ - abs 1 + ] cartesian-map ;
: <box-matrix> ( r -- matrix )
2 * 1 + dup '[ _ 1 <array> ] replicate ;
: <vandermonde-matrix> ( u n -- matrix )
<iota> [ v^n ] with map reverse flip ;
! Transformation matrices
:: <rotation-matrix3> ( axis theta -- matrix )
theta cos :> c
theta sin :> s
axis first3 :> ( x y z )
x sq 1.0 x sq - c * + x y * 1.0 c - * z s * - x z * 1.0 c - * y s * + 3array
x y * 1.0 c - * z s * + y sq 1.0 y sq - c * + y z * 1.0 c - * x s * - 3array
x z * 1.0 c - * y s * - y z * 1.0 c - * x s * + z sq 1.0 z sq - c * + 3array
3array ;
:: <rotation-matrix4> ( axis theta -- matrix )
theta cos :> c
theta sin :> s
axis first3 :> ( x y z )
x sq 1.0 x sq - c * + x y * 1.0 c - * z s * - x z * 1.0 c - * y s * + 0 4array
x y * 1.0 c - * z s * + y sq 1.0 y sq - c * + y z * 1.0 c - * x s * - 0 4array
x z * 1.0 c - * y s * - y z * 1.0 c - * x s * + z sq 1.0 z sq - c * + 0 4array
{ 0.0 0.0 0.0 1.0 } 4array ;
:: <translation-matrix4> ( offset -- matrix )
offset first3 :> ( x y z )
{
{ 1.0 0.0 0.0 x }
{ 0.0 1.0 0.0 y }
{ 0.0 0.0 1.0 z }
{ 0.0 0.0 0.0 1.0 }
} ;
<PRIVATE
GENERIC: >scale-factors ( object -- x y z )
M: number >scale-factors
dup dup ;
M: sequence >scale-factors
first3 ;
PRIVATE>
:: <scale-matrix3> ( factors -- matrix )
factors >scale-factors :> ( x y z )
{
{ x 0.0 0.0 }
{ 0.0 y 0.0 }
{ 0.0 0.0 z }
} ;
:: <scale-matrix4> ( factors -- matrix )
factors >scale-factors :> ( x y z )
{
{ x 0.0 0.0 0.0 }
{ 0.0 y 0.0 0.0 }
{ 0.0 0.0 z 0.0 }
{ 0.0 0.0 0.0 1.0 }
} ;
: <ortho-matrix4> ( factors -- matrix )
[ recip ] map <scale-matrix4> ;
:: <frustum-matrix4> ( xy-dim near far -- matrix )
xy-dim first2 :> ( x y )
near x /f :> xf
near y /f :> yf
near far + near far - /f :> zf
2 near far * * near far - /f :> wf
{
{ xf 0.0 0.0 0.0 }
{ 0.0 yf 0.0 0.0 }
{ 0.0 0.0 zf wf }
{ 0.0 0.0 -1.0 0.0 }
} ;
:: <skew-matrix4> ( theta -- matrix )
theta tan :> zf
{
{ 1.0 0.0 0.0 0.0 }
{ 0.0 1.0 0.0 0.0 }
{ 0.0 zf 1.0 0.0 }
{ 0.0 0.0 0.0 1.0 }
} ;
! a simpler verison of this, like matrix-map -except, but map-index, should be possible
: cartesian-matrix-map ( matrix quot: ( ... pair matrix -- ... matrix' ) -- matrix-seq )
[ [ first length <cartesian-square-indices> ] keep ] dip
'[ _ @ ] matrix-map ; inline
: cartesian-column-map ( matrix quot: ( ... pair matrix -- ... matrix' ) -- matrix-seq )
[ cols first2 ] prepose cartesian-matrix-map ; inline
! -------------------------------------------------
! numerical analysis of matrices follows
<PRIVATE
: square-rank ( square-matrix -- rank ) ;
: nonsquare-rank ( matrix -- rank ) ;
PRIVATE>
GENERIC: rank ( matrix -- rank )
M: zero-matrix rank
drop +zero-rank+ ;
M: square-matrix rank
square-rank ;
M: matrix rank
nonsquare-rank ;
GENERIC: nullity ( matrix -- nullity )
! implementation details of determinant and inverse
<PRIVATE
: alternating-sign ( seq odd-elts? -- seq' )
'[ even? _ = [ neg ] unless ] map-index ;
! the determinant of a 1x1 matrix is the value itself
! this works for any-dimensional matrices too
: (1determinant) ( matrix -- 1det ) flatten first ; inline
! optimized to find the determinant of a 2x2 matrix
: (2determinant) ( matrix -- 2det )
! multiply the diagonals and subtract
[ main-diagonal ] [ anti-diagonal ] bi [ first2 * ] bi@ - ; inline
! optimized for 3x3
! https://www.mathsisfun.com/algebra/matrix-determinant.html
:: (3determinant) ( matrix-seq -- 3det )
! first 3 elements of row 1
matrix-seq first first3 :> ( a b c )
! last 2 rows, transposed to make the next step easier
matrix-seq rest transpose
! get the lower sub-matrices in reverse order of a b c columns
[ rest ] [ [ first ] [ third ] bi 2array ] [ 1 head* ] tri 3array
! find determinants
[ (2determinant) ] map
! negate odd elements of a b c and multiply by the new determinants
{ a b c } t alternating-sign v*
! sum the resulting sequence
sum ;
DEFER: (ndeterminant)
: make-determinants ( n matrix -- seq )
<repetition> [
cols-except [ length ] keep (ndeterminant) ! recurses here
] map-index ;
DEFER: (determinant)
! generalized to 4 and higher
: (ndeterminant) ( n matrix -- ndet )
! TODO? recurse for n < 3
over 4 < [ (determinant) ] [
[ nip first t alternating-sign ] [ rest make-determinants ] 2bi
v* sum
] if ;
! switches on dimensions only
: (determinant) ( n matrix -- determinant )
over {
{ 1 [ nip (1determinant) ] }
{ 2 [ nip (2determinant) ] }
{ 3 [ nip (3determinant) ] }
[ drop (ndeterminant) ]
} case ;
PRIVATE>
GENERIC: determinant ( matrix -- determinant )
M: zero-square-matrix determinant
drop 0 ;
M: square-matrix determinant
[ length ] keep (determinant) ;
! determinant is undefined for m =/= n, unlike inverse
M: matrix determinant
dimension first2 non-square-determinant ;
: 1/det ( matrix -- 1/det )
determinant recip ; inline
! -----------------------------------------------------
! inverse operations and implementations follow
ALIAS: multiplicative-inverse recip
! per element, find the determinant of all other elements except the element's row / col
! https://www.mathsisfun.com/algebra/matrix-inverse-minors-cofactors-adjugate.html
: >minors ( matrix -- matrix' )
matrix-except-all [ [ determinant ] map ] map ;
! alternately invert values of the matrix (see alternating-sign)
: >cofactors ( matrix -- matrix' )
[ even? alternating-sign ] map-index ;
! multiply a matrix by the inverse of its determinant
: m*1/det ( matrix -- matrix' )
[ 1/det ] keep n*m ; inline
! inverse implementation
<PRIVATE
! https://www.mathsisfun.com/algebra/matrix-inverse-minors-cofactors-adjugate.html
: (square-inverse) ( square-matrix -- inverted )
! inverse of the determinant of the input matrix
[ 1/det ]
! adjugate of the cofactors of the matrix of minors
[ >minors >cofactors transpose ]
! adjugate * 1/det
bi n*m ;
! TODO
: (left-inverse) ( matrix -- left-invert ) ;
: (right-inverse) ( matrix -- right-invert ) ;
! TODO update this when rank works properly
! only defined for rank(A) = rows(A) OR rank(A) = cols(M)
! https://en.wikipedia.org/wiki/Invertible_matrix
: (specialized-inverse) ( rect-matrix -- inverted )
dup [ rank ] [ dimension ] bi [ = ] with map {
{ { t f } [ (left-inverse) ] }
{ { f t } [ (right-inverse) ] }
[ no-case ]
} case ;
PRIVATE>
M: zero-square-matrix recip
; inline
M: square-matrix recip
(square-inverse) ; inline
M: zero-matrix recip
transpose ; inline ! TODO: error based on rankiness
M: matrix recip
(specialized-inverse) ; inline
! TODO: use the faster algorithm: [ determinant zero? ]
: invertible-matrix? ( matrix -- ? )
[ dimension first2 max <identity-matrix> ] keep
dup recip m. = ;
: linearly-independent-matrix? ( matrix -- ? ) ;
<PRIVATE
! this is the original definition of m^n as committed in 2012; it has not been lost
: (m^n) ( m n -- n )
make-bits over first length <identity-matrix>
[ [ dupd m. ] when [ dup m. ] dip ] reduce nip ;
PRIVATE>
! A^-1 is the inverse but other negative powers are nonsense
: m^n ( m n -- n ) {
{ [ dup -1 = ] [ drop recip ] }
{ [ dup 0 >= ] [ (m^n) ] }
[ negative-power-matrix ]
} cond ;
: n^m ( n m -- n ) swap m^n ; inline
: covariance-matrix-ddof ( matrix ddof -- cov )
'[ _ cov-ddof ] cartesian-column-map ; inline
: covariance-matrix ( matrix -- cov )
0 covariance-matrix-ddof ; inline
: sample-covariance-matrix ( matrix -- cov )
1 covariance-matrix-ddof ; inline
: population-covariance-matrix ( matrix -- cov ) 0 covariance-matrix-ddof ; inline

1
extra/math/matrices/extras/summary.txt Normal file
View File

@ -0,0 +1 @@
Matrix arithmetic - extra and miscellaneous words

2
extra/rosetta-code/bitmap/bitmap.factor
View File

@ -36,7 +36,7 @@ IN: rosetta-code.bitmap
! The storage functions ! The storage functions
: <raster-image> ( width height -- image ) : <raster-image> ( width height -- image )
zero-matrix [ drop { 0 0 0 } ] mmap ; <zero-matrix> [ drop { 0 0 0 } ] mmap ;
: fill-image ( {R,G,B} image -- image ) : fill-image ( {R,G,B} image -- image )
swap '[ drop _ ] mmap! ; swap '[ drop _ ] mmap! ;
: set-pixel ( {R,G,B} {i,j} image -- ) set-Mi,j ; inline : set-pixel ( {R,G,B} {i,j} image -- ) set-Mi,j ; inline

2
extra/rosetta-code/conjugate-transpose/conjugate-transpose.factor
View File

@ -38,4 +38,4 @@ IN: rosetta-code.conjugate-transpose
dup conj-t [ m. ] [ swap m. ] 2bi = ; dup conj-t [ m. ] [ swap m. ] 2bi = ;
: unitary-matrix? ( matrix -- ? ) : unitary-matrix? ( matrix -- ? )
[ dup conj-t m. ] [ length identity-matrix ] bi = ; [ dup conj-t m. ] [ length <identity-matrix> ] bi = ;