:: run-length-uncompress-bitmap4 ( byte-array m n -- byte-array' )
byte-array <sequence-parser> :> sp
- m 1 + n zero-matrix :> matrix
+ m 1 + n <zero-matrix> :> matrix
n 4 mod n + :> stride
0 :> i!
0 :> j!
:: run-length-uncompress-bitmap8 ( byte-array m n -- byte-array' )
byte-array <sequence-parser> :> sp
- m 1 + n zero-matrix :> matrix
+ m 1 + n <zero-matrix> :> matrix
n 4 mod n + :> stride
0 :> i!
0 :> j!
Slava Pestov
+Joe Groff
+Doug Coleman
+Cat Stevens
+++ /dev/null
-Slava Pestov
+++ /dev/null
-USING: help.markup help.syntax math sequences ;
-
-IN: math.matrices.elimination
-
-HELP: inverse
-{ $values { "matrix" sequence } }
-{ $description "Computes the multiplicative inverse of a matrix. Assuming the matrix is invertible." }
-{ $examples
- "A matrix multiplied by its inverse is the identity matrix."
- { $example
- "USING: kernel math.matrices math.matrices.elimination prettyprint ;"
- "{ { 3 4 } { 7 9 } } dup inverse m. 2 identity-matrix = ."
- "t"
- }
-} ;
-
-HELP: echelon
-{ $values { "matrix" sequence } { "matrix'" sequence } }
-{ $description "Computes the reduced row-echelon form of the matrix." } ;
-
-HELP: nonzero-rows
-{ $values { "matrix" sequence } { "matrix'" sequence } }
-{ $description "Removes all all-zero rows from the matrix" }
-{ $examples
- { $example
- "USING: math.matrices.elimination prettyprint ;"
- "{ { 0 0 } { 5 6 } { 0 0 } { 4 0 } } nonzero-rows ."
- "{ { 5 6 } { 4 0 } }"
- }
-} ;
-
-HELP: leading
-{ $values
- { "seq" sequence }
- { "n" "the index of the first match, or " { $snippet f } "." }
- { "elt" "the first non-zero element, or " { $snippet f } "." }
-}
-{ $description "Find the first non-zero element of a sequence." } ;
+++ /dev/null
-IN: math.matrices.elimination.tests
-USING: kernel math.matrices math.matrices.elimination
-tools.test sequences ;
-
-{
- {
- { 1 0 0 0 }
- { 0 1 0 0 }
- { 0 0 1 0 }
- { 0 0 0 1 }
- }
-} [
- {
- { 1 0 0 0 }
- { 0 1 0 0 }
- { 0 0 1 0 }
- { 0 0 0 1 }
- } echelon
-] unit-test
-
-{
- {
- { 1 0 0 0 }
- { 0 1 0 0 }
- { 0 0 1 0 }
- { 0 0 0 1 }
- }
-} [
- {
- { 1 0 0 0 }
- { 1 1 0 0 }
- { 1 0 1 0 }
- { 1 0 0 1 }
- } echelon
-] unit-test
-
-{
- {
- { 1 0 0 0 }
- { 0 1 0 0 }
- { 0 0 1 0 }
- { 0 0 0 1 }
- }
-} [
- {
- { 1 0 0 0 }
- { 1 1 0 0 }
- { 1 0 1 0 }
- { 1 1 0 1 }
- } echelon
-] unit-test
-
-{
- {
- { 1 0 0 0 }
- { 0 1 0 0 }
- { 0 0 1 0 }
- { 0 0 0 1 }
- }
-} [
- {
- { 1 0 0 0 }
- { 1 1 0 0 }
- { 1 1 0 1 }
- { 1 0 1 0 }
- } echelon
-] unit-test
-
-{
- {
- { 1 0 0 0 }
- { 0 1 0 0 }
- { 0 0 0 0 }
- { 0 0 0 0 }
- }
-} [
- {
- { 0 1 0 0 }
- { 1 0 0 0 }
- { 1 0 0 0 }
- { 1 0 0 0 }
- } [
- [ 1 ] [ 0 0 pivot-row ] unit-test
- 1 0 do-row
- ] with-matrix
-] unit-test
-
-{
- {
- { 1 0 0 0 }
- { 0 1 0 0 }
- { 0 0 0 0 }
- { 0 0 0 0 }
- }
-} [
- {
- { 0 1 0 0 }
- { 1 0 0 0 }
- { 1 0 0 0 }
- { 1 0 0 0 }
- } echelon
-] unit-test
-
-{
- {
- { 1 0 0 0 }
- { 0 1 0 0 }
- { 0 0 0 1 }
- { 0 0 0 0 }
- }
-} [
- {
- { 1 0 0 0 }
- { 0 1 0 0 }
- { 1 0 0 1 }
- { 1 0 0 1 }
- } echelon
-] unit-test
-
-{
- {
- { 1 0 0 1 }
- { 0 1 0 1 }
- { 0 0 0 -1 }
- { 0 0 0 0 }
- }
-} [
- {
- { 0 1 0 1 }
- { 1 0 0 1 }
- { 1 0 0 0 }
- { 1 1 0 1 }
- } echelon
-] unit-test
-
-{
- 2
-} [
- {
- { 0 0 }
- { 0 0 }
- } nullspace length
-] unit-test
-
-{
- 1 3
-} [
- {
- { 0 1 0 1 }
- { 1 0 0 1 }
- { 1 0 0 0 }
- { 1 1 0 1 }
- } null/rank
-] unit-test
-
-{
- 1 3
-} [
- {
- { 0 0 0 0 0 1 0 1 }
- { 0 0 0 0 1 0 0 1 }
- { 0 0 0 0 1 0 0 0 }
- { 0 0 0 0 1 1 0 1 }
- } null/rank
-] unit-test
-
-{ { { 1 0 -1 } { 0 1 2 } } }
-[ { { 1 2 3 } { 4 5 6 } } solution ] unit-test
+++ /dev/null
-! Copyright (C) 2006, 2010 Slava Pestov.
-! See http://factorcode.org/license.txt for BSD license.
-USING: kernel locals math math.vectors math.matrices
-namespaces sequences fry sorting ;
-IN: math.matrices.elimination
-
-SYMBOL: matrix
-
-: with-matrix ( matrix quot -- )
- matrix swap [ matrix get ] compose with-variable ; inline
-
-: nth-row ( row# -- seq ) matrix get nth ;
-
-: change-row ( ..a row# quot: ( ..a seq -- ..b seq ) -- ..b )
- matrix get swap change-nth ; inline
-
-: exchange-rows ( row# row# -- ) matrix get exchange ;
-
-: rows ( -- n ) matrix get length ;
-
-: cols ( -- n ) 0 nth-row length ;
-
-: skip ( i seq quot -- n )
- over [ find-from drop ] dip swap [ ] [ length ] ?if ; inline
-
-: first-col ( row# -- n )
- ! First non-zero column
- 0 swap nth-row [ zero? not ] skip ;
-
-: clear-scale ( col# pivot-row i-row -- n )
- overd nth dup zero? [
- 3drop 0
- ] [
- [ nth dup zero? ] dip swap [
- 2drop 0
- ] [
- swap / neg
- ] if
- ] if ;
-
-: (clear-col) ( col# pivot-row i -- )
- [ [ clear-scale ] 2keep [ n*v ] dip v+ ] change-row ;
-
-: rows-from ( row# -- slice )
- rows dup <iota> <slice> ;
-
-: clear-col ( col# row# rows -- )
- [ nth-row ] dip [ [ 2dup ] dip (clear-col) ] each 2drop ;
-
-: do-row ( exchange-with row# -- )
- [ exchange-rows ] keep
- [ first-col ] keep
- dup 1 + rows-from clear-col ;
-
-: find-row ( row# quot -- i elt )
- [ rows-from ] dip find ; inline
-
-: pivot-row ( col# row# -- n )
- [ dupd nth-row nth zero? not ] find-row 2nip ;
-
-: (echelon) ( col# row# -- )
- over cols < over rows < and [
- 2dup pivot-row [ over do-row 1 + ] when*
- [ 1 + ] dip (echelon)
- ] [
- 2drop
- ] if ;
-
-: echelon ( matrix -- matrix' )
- [ 0 0 (echelon) ] with-matrix ;
-
-: nonzero-rows ( matrix -- matrix' )
- [ [ zero? ] all? ] reject ;
-
-: null/rank ( matrix -- null rank )
- echelon dup length swap nonzero-rows length [ - ] keep ;
-
-: leading ( seq -- n elt ) [ zero? not ] find ;
-
-: reduced ( matrix' -- matrix'' )
- [
- rows <iota> <reversed> [
- dup nth-row leading drop
- [ swap dup <iota> clear-col ] [ drop ] if*
- ] each
- ] with-matrix ;
-
-:: basis-vector ( row col# -- )
- row clone :> row'
- col# row' nth neg recip :> a
- 0 col# row' set-nth
- a row n*v col# matrix get set-nth ;
-
-: nullspace ( matrix -- seq )
- echelon reduced dup empty? [
- dup first length identity-matrix [
- [
- dup leading drop
- [ basis-vector ] [ drop ] if*
- ] each
- ] with-matrix flip nonzero-rows
- ] unless ;
-
-: 1-pivots ( matrix -- matrix )
- [ dup leading nip [ recip v*n ] when* ] map ;
-
-: solution ( matrix -- matrix )
- echelon nonzero-rows reduced 1-pivots ;
-
-: inverse ( matrix -- matrix ) ! Assumes an invertible matrix
- dup length
- [ identity-matrix [ append ] 2map solution ] keep
- [ tail ] curry map ;
+++ /dev/null
-Solving systems of linear equations
-USING: help.markup help.syntax math sequences ;
+! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens.
+USING: accessors arrays assocs generic.single formatting locals help.markup help.markup.private help.syntax io
+kernel math math.functions math.order math.ratios math.vectors opengl.gl prettyprint
+sequences sequences.generalizations urls ;
IN: math.matrices
+<PRIVATE
+! like $subsections but skip the extra blank line
+: $subs-nobl ( children -- )
+ [ $subsection* ] each ;
+
+! swapped and n-matrix variants
+: $equiv-word-note ( children -- )
+ [ "This word is the " ] dip
+ first2
+ " variant of " swap
+ [ { $link } ] dip suffix
+ "."
+ 5 narray print-element ;
+
+! words like <scale-matrix3> which have an array of inputs
+: $finite-input-note ( children -- )
+ [ "Only the first " ] dip
+ first2
+ " values in " swap
+ [ { $snippet } ] dip suffix
+ " are used."
+ 5 narray print-element ;
+
+! a note for when a word assumes a 2d matrix
+: $2d-only-note ( children -- )
+ drop { "This word is intended for use with \"flat\" (2-dimensional) matrices. "
+ ! "Using it with matrices of 3 or more dimensions may lead to unexpected results."
+ }
+ print-element ;
+
+! a note for numeric-specific operations
+: $matrix-scalar-note ( children -- )
+ \ $subs-nobl prefix
+ "This word assumes that elements of the input matrix are compatible with the following words:"
+ swap 2array
+ print-element ;
+
+: $keep-shape-note ( children -- )
+ drop { "The shape of the input matrix is preserved in the output." } print-element ;
+
+: $link2 ( children -- )
+ first2 swap [ write-link ] topic-span ;
+
+! so that we don't end up with multiple $notes calls leading to multiple Notes sections
+: $notelist ( children -- )
+ \ $list prefix $notes ;
+PRIVATE>
+
+ABOUT: "math.matrices"
+
+ARTICLE: "math.matrices" "Matrix operations"
+
+"The " { $vocab-link "math.matrices" } " vocabulary implements many ways of working with " { $emphasis "matrices" } " — sequences which have a minimum of 2 dimensions. Operations on 1-dimensional numeric vectors are implemented in " { $vocab-link "math.vectors" } ", upon which this vocabulary relies."
+$nl
+"In this vocabulary's documentation, " { $snippet "m" } " and " { $snippet "matrix" } " are the conventional names used for a given matrix object. " { $snippet "m" } " may also refer to a number."
+$nl
+"The " { $vocab-link "math.matrices.extras" } "vocabulary implements extensions to this one."
+$nl
+"Matrices are classified their mathematical properties, and by predicate words:"
+$nl
+! split up intentionally
+{ $subsections
+ matrix
+ irregular-matrix
+ square-matrix
+ zero-matrix
+ zero-square-matrix
+ null-matrix
+
+} { $subsections
+ matrix?
+ irregular-matrix?
+ square-matrix?
+ zero-matrix?
+ zero-square-matrix?
+ null-matrix?
+
+}
+
+"There are many ways to create 2-dimensional matrices:"
+{ $subsections
+ <matrix>
+ <matrix-by>
+ <matrix-by-indices>
+
+} { $subsections
+ <zero-matrix>
+ <zero-square-matrix>
+ <diagonal-matrix>
+ <anti-diagonal-matrix>
+ <identity-matrix>
+ <simple-eye>
+ <eye>
+
+} { $subsections
+ <coordinate-matrix>
+ <square-rows>
+ <square-cols>
+ <upper-matrix>
+ <lower-matrix>
+ <cartesian-square-indices>
+}
+
+"By-element mathematical operations on a matrix:"
+{ $subsections mneg m+n m-n m*n m/n n+m n-m n*m n/m }
+
+"By-element mathematical operations of two matricess:"
+{ $subsections m+ m- m* m/ m~ }
+
+"Dot product (multiplication) of vectors and matrices:"
+{ $subsections v.m m.v m. }
+
+"Transformations and elements of matrices:"
+{ $subsections
+ dimension
+ transpose anti-transpose
+ matrix-nth matrix-nths
+ matrix-set-nth matrix-set-nths
+
+} { $subsections
+ row rows rows-except
+ col cols cols-except
+
+} { $subsections
+ matrix-except matrix-except-all
+
+} { $subsections
+ matrix-map column-map stitch
+
+} { $subsections
+ main-diagonal
+ anti-diagonal
+
+} ;
+
+! PREDICATE CLASSES
+
+HELP: matrix
+{ $class-description "The class of regular, rectangular matrices. In mathematics and linear algebra, a matrix is a rectangular collection of scalar elements for the purpose of the uniform application of algorithms." }
+{ $notes "In Factor, any sequence with two or more dimensions (one or more layers of subsequences) can be a " { $link matrix } ", and the elements may be any " { $link object } "."
+$nl "A regular matrix is a sequence with two or more dimensions, whose subsequences are all of equal length. See " { $link regular-matrix? } "." }
+$nl "Irregular matrices are classified by " { $link irregular-matrix } "." ;
+
+HELP: irregular-matrix
+{ $class-description "The most common matrix, and most easily manipulated by this vocabulary, is rectangular. This predicate classifies irregular (non-rectangular) matrices." } ;
+
+HELP: square-matrix
+{ $class-description "The class of square matrices. A square matrix is a " { $link matrix } " which has the same number of rows and columns. In other words, its outermost two dimensions are of equal size." } ;
+
HELP: zero-matrix
-{ $values { "m" integer } { "n" integer } { "matrix" sequence } }
-{ $description "Creates a matrix of size " { $snippet "m x n" } ", filled with zeroes." } ;
+{ $class-description "The class of zero matrices. A zero matrix is a matrix whose only elements are the scalar " { $snippet "0" } "." }
+{ $notes "In mathematics, a zero-filled matrix is called a null matrix. In Factor, a "{ $link null-matrix } " is an empty matrix." } ;
-HELP: diagonal-matrix
-{ $values { "diagonal-seq" sequence } { "matrix" sequence } }
-{ $description "Creates a matrix with the specified diagonal values." } ;
+HELP: zero-square-matrix
+{ $class-description "The class of square zero matrices. This predicate is a composition of " { $link zero-matrix } " and " { $link square-matrix } "." } ;
-HELP: identity-matrix
-{ $values { "n" integer } { "matrix" sequence } }
-{ $description "Creates an identity matrix of size " { $snippet "n x n" } ", where the diagonal values are all ones." } ;
+HELP: null-matrix
+{ $class-description "The class of null matrices. A null matrix is an empty sequence, or a sequence which consists only of empty sequences." }
+{ $notes "In mathematics, a null matrix is a matrix full of zeroes. In Factor, such a matrix is called a " { $link zero-matrix } "." } ;
-HELP: m.v
-{ $values { "m" sequence } { "v" sequence } }
-{ $description "Computes the dot product between a matrix and a vector." }
+{ matrix irregular-matrix square-matrix zero-matrix null-matrix zero-square-matrix null-matrix } related-words
+
+! NON-PREDICATE TESTS
+
+HELP: regular-matrix?
+{ $values { "object" object } { "?" boolean } }
+{ $description "Tests if the object is a regular (well-formed, rectangular, etc) " { $link matrix } ". A regular matrix is a sequence with an equal number of elements in every row, and an equal number of elements in every column, such that there are no empty slots." }
+{ $notes "The " { $link null-matrix } " is considered regular, because of semantic requirements of the matrix implementation." }
{ $examples
- { $example
- "USING: math.matrices prettyprint ;"
- "{ { 1 -1 2 } { 0 -3 1 } } { 2 1 0 } m.v ."
- "{ 1 -3 }"
- }
+ "The example is an irregular matrix, because the rows have an unequal number of elements."
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 1 } { } } regular-matrix? ."
+ "f"
+ }
+ "The example is a regular matrix, because the rows have an equal number of elements."
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 1 } { 2 } } regular-matrix? ."
+ "t"
+ }
} ;
-HELP: m.
-{ $values { "m" sequence } }
-{ $description "Computes the dot product between two matrices, i.e multiplies them." }
+! BUILDERS
+HELP: <matrix>
+{ $values { "m" integer } { "n" integer } { "element" object } { "matrix" matrix } }
+{ $description "Creates a matrix of size " { $snippet "m x n" } ", filled with " { $snippet "element" } "." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "3 2 10 <matrix> ."
+ "{ { 10 10 } { 10 10 } { 10 10 } }"
+ }
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "4 1 \"¢\" <matrix> ."
+ "{ { \"¢\" } { \"¢\" } { \"¢\" } { \"¢\" } }"
+ }
+} ;
+
+HELP: <matrix-by>
+{ $values { "m" integer } { "n" integer } { "quot" { $quotation ( ... -- elt ) } } { "matrix" matrix } }
+{ $description "Creates a matrix of size " { $snippet "m x n" } " using elements given by " { $snippet "quot" } ", a quotation called to create each element." }
+{ $notes "The following are equivalent:"
+ { $code "m n [ 2drop foo ] <matrix-by-indices>" }
+ { $code "m n [ foo ] <matrix-by>" }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "4 5 [ 5 ] <matrix-by> ."
+ "{ { 5 5 5 5 5 } { 5 5 5 5 5 } { 5 5 5 5 5 } { 5 5 5 5 5 } }"
+ }
+} ;
+
+HELP: <matrix-by-indices>
+{ $values { "m" integer } { "n" integer } { "quot" { $quotation ( ... m' n' -- ... elt ) } } { "matrix" matrix } }
+{ $description "Creates an " { $snippet "m x n" } " " { $link matrix } " using elements given by " { $snippet "quot" } " . This word differs from " { $link <matrix-by> } " in that the indices are placed on the stack (in the same order) before " { $snippet "quot" } " runs. The output of the quotation will be the element at the given position in the matrix." }
+{ $notes "The following are equivalent:"
+ { $code "m n [ 2drop foo ] <matrix-by-indices>" }
+ { $code "m n [ foo ] <matrix-by>" }
+}
+{ $examples
+ { $example
+ "USING: math math.matrices prettyprint ;"
+ "3 4 [ * ] <matrix-by-indices> ."
+ "{ { 0 0 0 0 } { 0 1 2 3 } { 0 2 4 6 } }"
+ }
+} ;
+
+HELP: <zero-matrix>
+{ $values { "m" integer } { "n" integer } { "matrix" matrix } }
+{ $description "Creates a matrix of size " { $snippet "m x n" } ", filled with zeroes." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "2 3 <zero-matrix> ."
+ "{ { 0 0 0 } { 0 0 0 } }"
+ }
+} ;
+
+HELP: <zero-square-matrix>
+{ $values { "n" integer } { "matrix" matrix } }
+{ $description "Creates a matrix of size " { $snippet "n x n" } ", filled with zeroes. Shorthand for " { $code "n n <zero-matrix>" } "." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "2 <zero-square-matrix> ."
+ "{ { 0 0 } { 0 0 } }"
+ }
+} ;
+
+HELP: <diagonal-matrix>
+{ $values { "diagonal-seq" sequence } { "matrix" matrix } }
+{ $description "Creates a matrix with the specified main diagonal. This word has the opposite effect of " { $link anti-diagonal } "." }
+{ $notes "To use a diagonal starting in the lower right, reverse the input sequence before calling this word." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ 1 2 3 } <diagonal-matrix> ."
+ "{ { 1 0 0 } { 0 2 0 } { 0 0 3 } }"
+ }
+} ;
+
+HELP: <anti-diagonal-matrix>
+{ $values { "diagonal-seq" sequence } { "matrix" matrix } }
+{ $description "Creates a matrix with the specified anti-diagonal. This word has the opposite effect of " { $link main-diagonal } "." }
+{ $notes "To use a diagonal starting in the lower left, reverse the input sequence before calling this word." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ 1 2 3 } <anti-diagonal-matrix> ."
+ "{ { 0 0 1 } { 0 2 0 } { 3 0 0 } }"
+ }
+} ;
+
+HELP: <identity-matrix>
+{ $values { "n" integer } { "matrix" matrix } }
+{ $description "Creates an " { $url URL" http://enwp.org/Identity_matrix" "identity matrix" } " of size " { $snippet "n x n" } ", where the diagonal values are all ones." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "4 <identity-matrix> ."
+ "{ { 1 0 0 0 } { 0 1 0 0 } { 0 0 1 0 } { 0 0 0 1 } }"
+ }
+} ;
+
+HELP: <eye>
+{ $values { "m" integer } { "n" integer } { "k" integer } { "z" object } { "matrix" matrix } }
+{ $description "Creates an " { $snippet "m x n" } " matrix with a diagonal of " { $snippet "z" } " offset by " { $snippet "k" } " from the main diagonal. A positive value of " { $snippet "k" } " gives a diagonal above the main diagonal, whereas a negative value of " { $snippet "k" } " gives a diagonal below the main diagonal." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "5 6 0 4 <eye> ."
+ "{
+ { 4 0 0 0 0 0 }
+ { 0 4 0 0 0 0 }
+ { 0 0 4 0 0 0 }
+ { 0 0 0 4 0 0 }
+ { 0 0 0 0 4 0 }
+}"
+ }
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "5 5 2 2 <eye> ."
+ "{
+ { 0 0 2 0 0 }
+ { 0 0 0 2 0 }
+ { 0 0 0 0 2 }
+ { 0 0 0 0 0 }
+ { 0 0 0 0 0 }
+}"
+ }
+} ;
+
+HELP: <simple-eye>
+{ $values { "m" integer } { "n" integer } { "k" integer } { "matrix" matrix } }
+{ $description
+ "Creates an " { $snippet "m x n" } " matrix with a diagonal of ones offset by " { $snippet "k" } " from the main diagonal."
+ "The following are equivalent for any " { $snippet "m n k" } ":" { $code "m n k 1 <eye>" } { $code "m n k <simple-eye>" }
+ $nl
+ "Specify a different diagonal value with " { $link <eye> } "."
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "4 5 2 <simple-eye> ."
+ "{ { 0 0 1 0 0 } { 0 0 0 1 0 } { 0 0 0 0 1 } { 0 0 0 0 0 } }"
+ }
+} ;
+
+HELP: <coordinate-matrix>
+{ $values { "dim" pair } { "coordinates" matrix } }
+{ $description "Create a matrix in which each element is its own coordinate pair, also called a " { $link cartesian-product } "." }
+{ $notelist
+ { $equiv-word-note "non-square" <cartesian-square-indices> }
+ { $finite-input-note "two" "dim" }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ 2 4 } <coordinate-matrix> ."
+"{
+ { { 0 0 } { 0 1 } { 0 2 } { 0 3 } }
+ { { 1 0 } { 1 1 } { 1 2 } { 1 3 } }
+}"
+ }
+} ;
+
+HELP: <cartesian-indices>
+{ $values { "dim" pair } { "coordinates" matrix } }
+{ $description "An alias for " { $link <coordinate-matrix> } " which serves as the logical non-square companion to " { $link <cartesian-square-indices> } "." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ 2 4 } <cartesian-indices> ."
+"{
+ { { 0 0 } { 0 1 } { 0 2 } { 0 3 } }
+ { { 1 0 } { 1 1 } { 1 2 } { 1 3 } }
+}"
+ }
+} ;
+
+HELP: <cartesian-square-indices>
+{ $values { "n" integer } { "matrix" square-matrix } }
+{ $description "Create a " { $link square-matrix } " full of " { $link cartesian-product } "s. See " { $url URL" https://en.wikipedia.org/wiki/Cartesian_product" "cartesian product" } "." }
+{ $notes
+ { $equiv-word-note "square" <cartesian-indices> }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "1 <cartesian-square-indices> ."
+ "{ { { 0 0 } } }"
+ }
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "3 <cartesian-square-indices> ."
+"{
+ { { 0 0 } { 0 1 } { 0 2 } }
+ { { 1 0 } { 1 1 } { 1 2 } }
+ { { 2 0 } { 2 1 } { 2 2 } }
+}"
+ }
+} ;
+
+HELP: <square-rows>
+{ $values { "desc" { $or sequence integer matrix } } { "matrix" matrix } }
+{ $contract "Generate a " { $link square-matrix } " from a descriptor." }
+{ $description "If the descriptor is an " { $link integer } ", it is used to generate square rows within that range." $nl "If it is a 1-dimensional sequence, it is " { $link replicate } "d to create each row." $nl "If it is a " { $link matrix } ", it is cropped into a " { $link square-matrix } "." $nl "If it is a " { $link square-matrix } ", it is returned unchanged." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "3 <square-rows> ."
+ "{ { 0 1 2 } { 0 1 2 } { 0 1 2 } }"
+ }
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ 2 3 5 } <square-rows> ."
+ "{ { 2 3 5 } { 2 3 5 } { 2 3 5 } }"
+ }
+} ;
+
+HELP: <square-cols>
+{ $values { "desc" { $or sequence integer matrix } } { "matrix" matrix } }
+{ $contract "Generate a " { $link square-matrix } " from a descriptor." }
+{ $description "If the descriptor is an " { $link integer } ", it is used to generate square columns within that range." $nl "If it is a 1-dimensional sequence, it is " { $link replicate } "d to create each column." $nl "If it is a " { $link matrix } ", it is cropped into a " { $link square-matrix } "." $nl "If it is a " { $link square-matrix } ", it is returned unchanged." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "3 <square-cols> ."
+ "{ { 0 0 0 } { 1 1 1 } { 2 2 2 } }"
+ }
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ 2 3 5 } <square-cols> ."
+ "{ { 2 2 2 } { 3 3 3 } { 5 5 5 } }"
+ }
+} ;
+
+HELP: <lower-matrix>
+{ $values { "object" object } { "m" integer } { "n" integer } { "matrix" matrix } }
+{ $description "Make a lower triangular matrix, where all the values above the main diagonal are " { $snippet "0" } ". " { $snippet "object" } " will be used as the value for the nonzero part of the matrix, while " { $snippet "m" } " and " { $snippet "n" } " are used as the dimensions. The inverse of this word is " { $link <upper-matrix> } ". See " { $url URL" https://en.wikipedia.org/wiki/Triangular_matrix" "triangular matrix" } "." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "1 5 5 <lower-matrix> ."
+"{
+ { 1 0 0 0 0 }
+ { 1 1 0 0 0 }
+ { 1 1 1 0 0 }
+ { 1 1 1 1 0 }
+ { 1 1 1 1 1 }
+}"
+ }
+} ;
+
+HELP: <upper-matrix>
+{ $values { "object" object } { "m" integer } { "n" integer } { "matrix" matrix } }
+{ $description "Make an upper triangular matrix, where all the values below the main diagonal are " { $snippet "0" } ". " { $snippet "object" } " will be used as the value for the nonzero part of the matrix, while " { $snippet "m" } " and " { $snippet "n" } " are used as the dimensions. The inverse of this word is " { $link <lower-matrix> } ". See " { $url URL" https://en.wikipedia.org/wiki/Triangular_matrix" "triangular matrix" } "." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "1 5 5 <upper-matrix> ."
+"{
+ { 1 1 1 1 1 }
+ { 0 1 1 1 1 }
+ { 0 0 1 1 1 }
+ { 0 0 0 1 1 }
+ { 0 0 0 0 1 }
+}"
+ }
+} ;
+
+HELP: stitch
+{ $values { "m" matrix } { "m'" matrix } }
+{ $description "Folds an " { $snippet "n>2" } "-dimensional matrix onto itself." }
+{ $examples
+ { $unchecked-example
+ "USING: math.matrices prettyprint ;"
+"{
+ { { 0 5 } { 6 7 } { 0 15 } { 18 21 } }
+ { { 0 10 } { 12 14 } { 0 20 } { 24 28 } }
+} stitch ."
+"{
+ { 0 5 0 10 }
+ { 6 7 12 14 }
+ { 0 15 0 20 }
+ { 18 21 24 28 }
+}"
+ }
+} ;
+
+HELP: row
+{ $values { "n" integer } { "matrix" matrix } { "row" sequence } }
+{ $description "Get the nth row of the matrix." }
+{ $notes "Like most Factor sequences, indexing is 0-based. The first row is given by " { $snippet "m 0 row" } "." }
+{ $examples
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "{ { 1 2 } { 3 4 } } 1 swap row ."
+ "{ 3 4 }"
+ }
+} ;
+
+HELP: rows
+{ $values { "seq" sequence } { "matrix" matrix } { "rows" sequence } }
+{ $description "Get the rows from " { $snippet "matrix" } " listed by " { $snippet "seq" } "." }
+{ $notelist { $equiv-word-note "multiplexing" row } }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ 0 1 } { { 1 2 } { 3 4 } } rows ."
+ "{ { 1 2 } { 3 4 } }"
+ }
+} ;
+
+HELP: col
+{ $values { "n" integer } { "matrix" matrix } { "col" sequence } }
+{ $description "Get the " { $snippet "n" } "th column of the matrix." }
+{ $notes "Like most Factor sequences, indexing is 0-based. The first column is given by " { $snippet "m 0 col" } "." }
+{ $examples
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "{ { 1 2 } { 3 4 } } 1 swap col ."
+ "{ 2 4 }"
+ }
+} ;
+
+HELP: cols
+{ $values { "seq" sequence } { "matrix" matrix } { "cols" sequence } }
+{ $description "Get the columns from " { $snippet "matrix" } " listed by " { $snippet "seq" } "." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ 0 1 } { { 1 2 } { 3 4 } } cols ."
+ "{ { 1 3 } { 2 4 } }"
+ }
+} ;
+
+HELP: matrix-map
+{ $values { "matrix" matrix } { "quot" { $quotation ( ... elt -- ... elt' ) } } { "matrix'" matrix } }
+{ $description "Apply the quotation to every element of the matrix." }
+{ $notelist $2d-only-note }
+{ $examples
+ { $example
+ "USING: math.matrices kernel math prettyprint ;"
+ "3 <identity-matrix> [ zero? 15 -8 ? ] matrix-map ."
+ "{ { -8 15 15 } { 15 -8 15 } { 15 15 -8 } }"
+ }
+} ;
+
+HELP: column-map
+{ $values { "matrix" matrix } { "quot" { $quotation ( ... col -- ... col' ) } } { "matrix'" { $maybe sequence matrix } } }
+{ $description "Apply the quotation to every column of the matrix. The output of the quotation must be a sequence." }
+{ $notelist $2d-only-note { $equiv-word-note "transpose" map } }
{ $examples
- { $example
- "USING: math.matrices prettyprint ;"
- "{ { 1 -1 2 } { 0 -3 1 } } { { 3 7 } { 9 12 } } m. ."
- "{ { -6 -5 } { -27 -36 } }"
- }
+ { $example
+ "USING: sequences math.matrices prettyprint ;"
+ "3 <identity-matrix> [ reverse ] column-map ."
+ "{ { 0 0 1 } { 0 1 0 } { 1 0 0 } }"
+ }
+} ;
+
+HELP: matrix-nth
+{ $values { "pair" pair } { "matrix" matrix } { "elt" object } }
+{ $description "Retrieve the element in the matrix at the zero-indexed " { $snippet "row, column" } " pair." }
+{ $notelist { $equiv-word-note "two-dimensional" nth } $2d-only-note }
+{ $errors { $list
+ { { $link bounds-error } " if the first element in " { $snippet "pair" } " is greater than the maximum row index in " { $snippet "matrix" } }
+ { { $link bounds-error } " if the second element in " { $snippet "pair" } " is greater than the maximum column index in " { $snippet "matrix" } }
+} }
+{ $examples
+ "Get the entry at row 1, column 0."
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ 1 0 } { { 0 1 } { 2 3 } } matrix-nth ."
+ "2"
+ }
+} ;
+
+HELP: matrix-nths
+{ $values { "pairs" assoc } { "matrix" matrix } { "elts" sequence } }
+{ $description "Retrieve all the elements in the matrix at each of the zero-indexed " { $snippet "row, column" } " pairs in " { $snippet "pairs" } "." }
+{ $notelist { $equiv-word-note "two-dimensional" nths } $2d-only-note }
+{ $errors { $list
+ { { $link bounds-error } " if the first element of a pair in " { $snippet "pairs" } " is greater than the maximum row index in " { $snippet "matrix" } }
+ { { $link bounds-error } " if the second element of a pair in " { $snippet "pairs" } " is greater than the maximum column index in " { $snippet "matrix" } }
+} }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 1 0 } { 1 1 } } { { 0 1 } { 2 3 } } matrix-nths ."
+ "{ 2 3 }"
+ }
+} ;
+
+HELP: matrix-set-nth
+{ $values { "obj" object } { "pair" pair } { "matrix" matrix } }
+{ $description "Set the element in the matrix at the 2D index given by " { $snippet "pair" } " to " { $snippet "obj" } ". This operation is destructive." }
+{ $side-effects "matrix" }
+{ $notelist { $equiv-word-note "two-dimensional" set-nth } $2d-only-note }
+{ $errors { $list
+ { { $link bounds-error } " if the first element of a pair in " { $snippet "pairs" } " is greater than the maximum row index in " { $snippet "matrix" } }
+ { { $link bounds-error } " if the second element of a pair in " { $snippet "pairs" } " is greater than the maximum column index in " { $snippet "matrix" } }
+ "Throws an error if the sequence cannot hold elements of the given type."
+} }
+{ $examples
+ "Change the entry at row 1, column 0."
+ { $example
+ "USING: math.matrices kernel prettyprint ;"
+ "{ { 0 1 } { 2 3 } } \"a\" { 1 0 } pick matrix-set-nth ."
+ "{ { 0 1 } { \"a\" 3 } }"
+ }
+} ;
+
+HELP: matrix-set-nths
+{ $values { "obj" object } { "pairs" assoc } { "matrix" matrix } }
+{ $description "Applies " { $link matrix-set-nth } " to " { $snippet "matrix" } " for each " { $snippet "row, column" } " pair in " { $snippet "pairs" } ", setting the elements to " { $snippet "obj" } "." }
+{ $side-effects "matrix" }
+{ $notelist { $equiv-word-note "multiplexing" matrix-set-nth } $2d-only-note }
+{ $errors { $list
+ { { $link bounds-error } " if the first element of a pair in " { $snippet "pairs" } " is greater than the maximum row index in " { $snippet "matrix" } }
+ { { $link bounds-error } " if the second element of a pair in " { $snippet "pairs" } " is greater than the maximum column index in " { $snippet "matrix" } }
+ "Throws an error if the sequence cannot hold elements of the given type."
+} }
+{ $examples
+ "Change both entries on row 0."
+ { $example
+ "USING: math.matrices kernel prettyprint ;"
+ "{ { 0 1 } { 2 3 } } \"a\" { { 1 0 } { 1 1 } } pick matrix-set-nths ."
+ "{ { 0 1 } { \"a\" \"a\" } }"
+ }
+} ;
+
+
+HELP: mneg
+{ $values { "m" matrix } { "m'" matrix } }
+{ $description "Negate (invert the sign) of every element in the matrix. The resulting matrix is called the " { $emphasis "additive inverse" } " of the input matrix." }
+{ $notelist
+ { $equiv-word-note "companion" mabs }
+ $2d-only-note
+ { $matrix-scalar-note neg }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 5 9 } { 15 -17 } } mneg ."
+ "{ { -5 -9 } { -15 17 } }"
+ }
+} ;
+
+HELP: mabs
+{ $values { "m" matrix } { "m'" matrix } }
+{ $description "Compute the absolute value (" { $link abs } ") of each element in the matrix." }
+{ $notelist
+ { $equiv-word-note "companion" mneg }
+ $2d-only-note
+ { $matrix-scalar-note abs }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { -5 -9 } { -15 17 } } mabs ."
+ "{ { 5 9 } { 15 17 } }"
+ }
+} ;
+
+HELP: n+m
+{ $values { "n" object } { "m" matrix } }
+{ $description { $snippet "n" } " is treated as a scalar and added to each element of the matrix " { $snippet "m" } "." }
+{ $notelist
+ { $equiv-word-note "swapped" m+n }
+ $2d-only-note
+ { $matrix-scalar-note + }
+}
+{ $examples
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "1 3 <identity-matrix> n+m ."
+ "{ { 2 1 1 } { 1 2 1 } { 1 1 2 } }"
+ }
+} ;
+
+HELP: m+n
+{ $values { "m" matrix } { "n" object } }
+{ $description { $snippet "n" } " is treated as a scalar and added to each element of the matrix " { $snippet "m" } "." }
+{ $notelist
+ { $equiv-word-note "swapped" n+m }
+ $2d-only-note
+ { $matrix-scalar-note + }
+}
+{ $examples
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "3 <identity-matrix> 1 m+n ."
+ "{ { 2 1 1 } { 1 2 1 } { 1 1 2 } }"
+ }
+} ;
+
+HELP: n-m
+{ $values { "n" object } { "m" matrix } }
+{ $description { $snippet "n" } " is treated as a scalar and subtracted from each element of the matrix " { $snippet "m" } "." }
+{ $notelist
+ { $equiv-word-note "swapped" m-n }
+ $2d-only-note
+ { $matrix-scalar-note - }
+}
+{ $examples
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "1 3 <identity-matrix> n-m ."
+ "{ { 0 1 1 } { 1 0 1 } { 1 1 0 } }"
+ }
+} ;
+
+HELP: m-n
+{ $values { "m" matrix } { "n" object } }
+{ $description { $snippet "n" } " is treated as a scalar and subtracted from each element of the matrix " { $snippet "m" } "." }
+{ $notelist
+ { $equiv-word-note "swapped" n-m }
+ $2d-only-note
+ { $matrix-scalar-note - }
+}
+{ $examples
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "3 <identity-matrix> 1 m-n ."
+ "{ { 0 -1 -1 } { -1 0 -1 } { -1 -1 0 } }"
+ }
+} ;
+
+HELP: n*m
+{ $values { "n" object } { "m" matrix } }
+{ $description "Every element in the input matrix " { $snippet "m" } " is multiplied by the scalar "{ $snippet "n" } "." }
+{ $notelist
+ $keep-shape-note
+ { $equiv-word-note "swapped" m*n }
+ $2d-only-note
+ { $matrix-scalar-note * }
+}
+{ $examples
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "3 3 <identity-matrix> n*m ."
+ "{ { 3 0 0 } { 0 3 0 } { 0 0 3 } }"
+ }
+} ;
+
+HELP: m*n
+{ $values { "m" matrix } { "n" object } }
+{ $description "Every element in the input matrix " { $snippet "m" } " is multiplied by the scalar "{ $snippet "n" } "." }
+{ $notelist
+ $keep-shape-note
+ { $equiv-word-note "swapped" n*m }
+ $2d-only-note
+ { $matrix-scalar-note * }
+}
+
+{ $examples
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "3 <identity-matrix> 3 m*n ."
+ "{ { 3 0 0 } { 0 3 0 } { 0 0 3 } }"
+ }
+} ;
+
+HELP: n/m
+{ $values { "n" object } { "m" matrix } }
+{ $description "Every element in the input matrix " { $snippet "m" } " is divided by the scalar "{ $snippet "n" } "." }
+{ $notelist
+ $keep-shape-note
+ { $equiv-word-note "swapped" m/n }
+ $2d-only-note
+ { $matrix-scalar-note / }
+}
+{ $examples
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "2 { { 4 5 } { 2 1 } } n/m ."
+ "{ { 1/2 2/5 } { 1 2 } }"
+ }
+} ;
+
+HELP: m/n
+{ $values { "m" matrix } { "n" object } }
+{ $description "Every element in the input matrix " { $snippet "m" } " is divided by the scalar "{ $snippet "n" } "." }
+{ $notelist
+ $keep-shape-note
+ { $equiv-word-note "swapped" n/m }
+ $2d-only-note
+ { $matrix-scalar-note / }
+}
+{ $examples
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "{ { 4 5 } { 2 1 } } 2 m/n ."
+ "{ { 2 2+1/2 } { 1 1/2 } }"
+ }
} ;
HELP: m+
-{ $values { "m" sequence } }
-{ $description "Adds the matrices component-wise." }
+{ $values { "m1" matrix } { "m2" matrix } { "m" matrix } }
+{ $description "Adds two matrices element-wise." }
+{ $notelist
+ $2d-only-note
+ { $matrix-scalar-note + }
+}
{ $examples
- { $example
- "USING: math.matrices prettyprint ;"
- "{ { 1 2 } { 3 4 } } { { 5 6 } { 7 8 } } m+ ."
- "{ { 6 8 } { 10 12 } }"
- }
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 1 2 3 } { 3 2 1 } } { { 4 5 6 } { 6 5 4 } } m+ ."
+ "{ { 5 7 9 } { 9 7 5 } }"
+ }
} ;
HELP: m-
-{ $values { "m" sequence } }
-{ $description "Subtracts the matrices component-wise." }
+{ $values { "m1" matrix } { "m2" matrix } { "m" matrix } }
+{ $description "Subtracts two matrices element-wise." }
+{ $notelist
+ $2d-only-note
+ { $matrix-scalar-note - }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 4 5 6 } { 6 5 4 } } { { 1 2 3 } { 3 2 1 } } m- ."
+ "{ { 3 3 3 } { 3 3 3 } }"
+ }
+} ;
+
+HELP: m*
+{ $values { "m1" matrix } { "m2" matrix } { "m" matrix } }
+{ $description "Multiplies two matrices element-wise." }
+{ $notelist
+ $2d-only-note
+ { $matrix-scalar-note * }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 5 9 } { 15 17 } } { { 3 2 } { 4 9 } } m* ."
+ "{ { 15 18 } { 60 153 } }"
+ }
+} ;
+
+HELP: m/
+{ $values { "m1" matrix } { "m2" matrix } { "m" matrix } }
+{ $description "Divides two matrices element-wise." }
+{ $notelist
+ $2d-only-note
+ { $matrix-scalar-note / }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 5 9 } { 15 17 } } { { 3 2 } { 4 9 } } m/ ."
+ "{ { 1+2/3 4+1/2 } { 3+3/4 1+8/9 } }"
+ }
+} ;
+
+HELP: m.v
+{ $values { "m" matrix } { "v" sequence } { "p" matrix } }
+{ $description "Computes the dot product of a matrix and a vector." }
+{ $notelist
+ { $equiv-word-note "swapped" v.m }
+ $2d-only-note
+ { $matrix-scalar-note * + }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 1 -1 2 } { 0 -3 1 } } { 2 1 0 } m.v ."
+ "{ 1 -3 }"
+ }
+} ;
+
+HELP: v.m
+{ $values { "v" sequence } { "m" matrix } { "p" matrix } }
+{ $description "Computes the dot product of a vector and a matrix." }
+{ $notelist
+ { $equiv-word-note "swapped" m.v }
+ $2d-only-note
+ { $matrix-scalar-note * + }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ 2 1 0 } { { 1 -1 2 } { 0 -3 1 } } v.m ."
+ "{ 2 -5 5 }"
+ }
+} ;
+
+HELP: m.
+{ $values { "m" matrix } }
+{ $description "Computes the dot product of two matrices, i.e multiplies them." }
+{ $notelist
+ $2d-only-note
+ { $matrix-scalar-note * + }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 1 -1 2 } { 0 -3 1 } } { { 3 7 } { 9 12 } } m. ."
+ "{ { -6 -5 } { -27 -36 } }"
+ }
+} ;
+
+HELP: m~
+{ $values { "m1" matrix } { "m2" matrix } { "epsilon" number } { "?" boolean } }
+{ $description "Compares the matrices like " { $link ~ } ", using the " { $snippet "epsilon" } "." }
+{ $notelist
+ $2d-only-note
+ { $matrix-scalar-note ~ }
+}
{ $examples
- { $example
- "USING: math.matrices prettyprint ;"
- "{ { 5 9 } { 15 17 } } { { 3 2 } { 4 9 } } m- ."
- "{ { 2 7 } { 11 8 } }"
- }
+ { "In the example, only " { $snippet ".01" } " was added to each element, so the new matrix is within the epsilon " { $snippet ".1" } "of the original." }
+ { $example
+ "USING: kernel math math.matrices prettyprint ;"
+ "{ { 5 9 } { 15 17 } } dup [ .01 + ] matrix-map .1 m~ ."
+ "t"
+ }
} ;
-HELP: kron
-{ $values { "m1" sequence } { "m2" sequence } { "m" sequence } }
-{ $description "Calculates the Kronecker product of two matrices." }
+HELP: mmin
+{ $values { "m" matrix } { "n" object } }
+{ $description "Determine the minimum value of the matrix." }
+{ $notelist
+ $2d-only-note
+ { $matrix-scalar-note min }
+}
{ $examples
- { $example "USING: math.matrices prettyprint ;"
- "{ { 1 2 } { 3 4 } } { { 0 5 } { 6 7 } } kron ."
- "{ { 0 5 0 10 } { 6 7 12 14 } { 0 15 0 20 } { 18 21 24 28 } }" }
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 5 9 } { 15 17 } } mmin ."
+ "5"
+ }
} ;
-HELP: outer
-{ $values { "u" sequence } { "v" sequence } { "m" sequence } }
-{ $description "Computers the outer product of " { $snippet "u" } " and " { $snippet "v" } "." }
+HELP: mmax
+{ $values { "m" matrix } { "n" object } }
+{ $description "Determine the maximum value of the matrix." }
+{ $notelist
+ $2d-only-note
+ { $matrix-scalar-note max }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 5 9 } { 15 17 } } mmax ."
+ "17"
+ }
+} ;
+
+HELP: mnorm
+{ $values { "m" matrix } { "m'" matrix } }
+{ $description "Normalize a matrix. Each element from the input matrix is computed as a fraction of the maximum element. The maximum element becomes " { $snippet "1/1" } "." }
+{ $notelist
+ $2d-only-note
+ { $matrix-scalar-note max abs / }
+}
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 5 9 } { 15 17 } } mnorm ."
+ "{ { 5/17 9/17 } { 15/17 1 } }"
+ }
+} ;
+
+HELP: m-infinity-norm
+{ $values { "m" matrix } { "n" number } } ;
+
+HELP: m-1norm
+{ $values { "m" matrix } { "n" number } } ;
+
+HELP: frobenius-norm
+{ $values { "m" matrix } { "n" number } }
+{ $notes "Also known as the Hilbert-Schmidt norm." } ;
+
+HELP: >square-matrix
+{ $values { "m" matrix } { "subset" square-matrix } }
+{ $description "Find only the " { $link2 square-matrix "square" } " subset of the input matrix." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 0 2 4 6 } { 1 3 5 7 } } >square-matrix ."
+ "{ { 0 2 } { 1 3 } }"
+ }
+} ;
+
+HELP: main-diagonal
+{ $values { "matrix" matrix } { "seq" sequence } }
+{ $description "Find the main diagonal of a matrix." $nl "This diagonal begins in the upper left of the matrix at index " { $snippet "{ 0 0 }" } ", continuing downward and rightward for all indices " { $snippet "{ n n }" } " in the " { $link square-matrix } " subset of the input (see " { $link <square-rows> } ")." }
+{ $notelist
+ { "If the number of rows in the square subset of the input is even, then this diagonal will not contain elements found in the " { $link anti-diagonal } ". However, if the size of the square subset is odd, then this diagonal will share at most one element with " { $link anti-diagonal } "." }
+ { "This diagonal is sometimes called the " { $emphasis "first diagonal" } "." }
+ { $equiv-word-note "opposite" anti-diagonal }
+}
+{ $examples
+ { "The operation is simple on a " { $link square-matrix } ":" }
+ { $example
+ "USING: math.matrices prettyprint ;"
+"{
+ { 7 2 11 }
+ { 9 7 7 }
+ { 1 8 0 }
+} main-diagonal ."
+ "{ 7 7 0 }"
+ }
+ "The square subset of the following input matrix consists of all rows but the last. The main diagonal does not include the last row because it has no fourth element."
+ { $example
+ "USING: math.matrices prettyprint ;"
+"{
+ { 6 5 0 }
+ { 7 2 6 }
+ { 4 3 9 }
+ { 3 3 3 }
+} main-diagonal ."
+ "{ 6 2 9 }"
+ }
+} ;
+
+HELP: anti-diagonal
+{ $values { "matrix" matrix } { "seq" sequence } }
+{ $description "Find the anti-diagonal of a matrix." $nl "This diagonal begins in the upper right of the matrix, continuing downward and leftward for all indices in the " { $link square-matrix } " subset of the input (see " { $link <square-rows> } ")." }
+{ $notelist
+ { "If the number of rows in the square subset of the input is even, then this diagonal will not contain elements found in the " { $link main-diagonal } ". However, if the size of the square subset is odd, then this diagonal will share at most one element with " { $link main-diagonal } "." }
+ { "This diagonal is sometimes called the " { $emphasis "second diagonal" } "." }
+ { $equiv-word-note "opposite" main-diagonal }
+}
+{ $examples
+ { "The operation is simple on a " { $link square-matrix } ":" }
+ { $example
+ "USING: math.matrices prettyprint ;"
+"{
+ { 7 2 11 }
+ { 9 7 7 }
+ { 1 8 0 }
+} anti-diagonal ."
+ "{ 11 7 1 }"
+ }
+ "The square subset of the following input matrix consists of all rows but the last. The anti-diagonal does not include the last row because it has no fourth element."
+ { $example
+ "USING: math.matrices prettyprint ;"
+"{
+ { 6 5 0 }
+ { 7 2 6 }
+ { 4 3 9 }
+ { 3 3 3 }
+} anti-diagonal ."
+ "{ 0 2 4 }"
+ }
+} ;
+
+
+HELP: transpose
+{ $values { "matrix" matrix } { "newmatrix" matrix } }
+{ $description "Transpose the input matrix over its " { $link main-diagonal } ". The main diagonal itself is preserved, whereas the anti-diagonal is reversed." }
+{ $notelist
+ { "This word is an alias for " { $link flip } ", so that it may be recognised as the common mathematical operation." }
+ { $equiv-word-note "opposite" anti-transpose }
+}
+{ $examples
+ { $example
+ "USING: math.matrices sequences prettyprint ;"
+ "5 <iota> <anti-diagonal-matrix> transpose ."
+"{
+ { 0 0 0 0 4 }
+ { 0 0 0 3 0 }
+ { 0 0 2 0 0 }
+ { 0 1 0 0 0 }
+ { 0 0 0 0 0 }
+}"
+ }
+} ;
+
+HELP: anti-transpose
+{ $values { "matrix" matrix } { "newmatrix" matrix } }
+{ $description "Like " { $link transpose } " except that the matrix is transposed over the " { $link anti-diagonal } ", so that the anti-diagonal itself is preserved and the " { $link main-diagonal } " is reversed." }
+{ $notes { $equiv-word-note "opposite" transpose } }
+{ $examples
+ { $example
+ "USING: math.matrices sequences prettyprint ;"
+ "5 <iota> <diagonal-matrix> anti-transpose ."
+"{
+ { 4 0 0 0 0 }
+ { 0 3 0 0 0 }
+ { 0 0 2 0 0 }
+ { 0 0 0 1 0 }
+ { 0 0 0 0 0 }
+}"
+ }
+} ;
+
+HELP: rows-except
+{ $values { "matrix" matrix } { "desc" { $or integer sequence } } { "others" matrix } }
+{ $contract "Get all the rows from " { $snippet "matrix" } " " { $emphasis "not" } " described by " { $snippet "desc" } "." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+"{
+ { 2 7 12 2 }
+ { 8 9 10 0 }
+ { 1 3 3 5 }
+ { 8 13 7 12 }
+} { 1 3 } rows-except ."
+ "{ { 2 7 12 2 } { 1 3 3 5 } }"
+ }
+} ;
+
+HELP: cols-except
+{ $values { "matrix" matrix } { "desc" { $or integer sequence } } { "others" matrix } }
+{ $contract "Get all the columns from " { $snippet "matrix" } " " { $emphasis "not" } " described by " { $snippet "desc" } "." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+"{
+ { 2 7 12 2 }
+ { 8 9 10 0 }
+ { 1 3 3 5 }
+ { 8 13 7 12 }
+} { 1 3 } cols-except . "
+ "{ { 2 12 } { 8 10 } { 1 3 } { 8 7 } }"
+ }
+} ;
+HELP: matrix-except
+{ $values { "matrix" matrix } { "exclude-pair" pair } { "submatrix" matrix } }
+{ $description "Get all the rows and columns from " { $snippet "matrix" } " except the row and column given in " { $snippet "exclude-pair" } ". The result is the " { $snippet "submatrix" } " containing no values from the given row and column." }
+{ $examples
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 0 1 } { 2 3 } } { 0 1 } matrix-except ."
+ "{ { 2 } }"
+ }
+} ;
+
+HELP: submatrix-excluding
+{ $values { "matrix" matrix } { "exclude-pair" pair } { "submatrix" matrix } }
+{ $description "A possibly more obvious word for " { $link matrix-except } "." } ;
+
+HELP: matrix-except-all
+{ $values { "matrix" matrix } { "submatrices" { $sequence matrix } } }
+{ $description "Find every possible submatrix of " { $snippet "matrix" } " by using " { $link matrix-except } " for every value's row-column pair." }
+{ $examples
+ "There are 9 possible 2x2 submatrices of a 3x3 matrix with 9 indices, because there are 9 indices to exclude creating a new submatrix."
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ { 0 1 2 } { 3 4 5 } { 6 7 8 } } matrix-except-all ."
+ "{
+ {
+ { { 4 5 } { 7 8 } }
+ { { 3 5 } { 6 8 } }
+ { { 3 4 } { 6 7 } }
+ }
+ {
+ { { 1 2 } { 7 8 } }
+ { { 0 2 } { 6 8 } }
+ { { 0 1 } { 6 7 } }
+ }
+ {
+ { { 1 2 } { 4 5 } }
+ { { 0 2 } { 3 5 } }
+ { { 0 1 } { 3 4 } }
+ }
+}"
+ }
+} ;
+
+HELP: all-submatrices
+{ $values { "matrix" matrix } { "submatrices" { $sequence matrix } } }
+{ $description "A possibly more obvious name for " { $link matrix-except-all } "." } ;
+
+HELP: dimension
+{ $values { "matrix" matrix } { "dimension" pair } }
+{ $description "Find the dimension of the input matrix, in the order of " { $snippet "{ rows cols }"} "." }
+{ $notelist $2d-only-note "Not to be confused with dimensionality, or the number of dimension scalars needed to describe a matrix." }
{ $examples
- { $example "USING: math.matrices prettyprint ;"
- "{ 5 6 7 } { 1 2 3 } outer ."
- "{ { 5 10 15 } { 6 12 18 } { 7 14 21 } }" }
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "4 30 1 <matrix> dimension ."
+ "{ 4 30 }"
+ }
+ { $example
+ "USING: math.matrices prettyprint ;"
+ "{ } dimension ."
+ "{ 0 0 }"
+ }
} ;
-USING: math.matrices math.vectors tools.test math kernel ;
-IN: math.matrices.tests
+! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens.
+USING: arrays combinators.short-circuit grouping kernel math math.matrices math.matrices.private
+math.statistics math.vectors sequences sequences.deep sets tools.test ;
+IN: math.matrices
+
+<PRIVATE
+: call-eq? ( obj quots -- ? )
+ [ call( x -- x ) ] with map all-eq? ; ! inline
+PRIVATE>
+! ------------------------
+! predicates
+
+{ t } [ { } regular-matrix? ] unit-test
+{ t } [ { { } } regular-matrix? ] unit-test
+{ t } [ { { 1 2 } } regular-matrix? ] unit-test
+{ t } [ { { 1 2 } { 3 4 } } regular-matrix? ] unit-test
+{ t } [ { { 1 } { 3 } } regular-matrix? ] unit-test
+{ f } [ { { 1 2 } { 3 } } regular-matrix? ] unit-test
+{ f } [ { { 1 } { 3 2 } } regular-matrix? ] unit-test
-{
- { { 0 } { 0 } { 0 } }
-} [
- 3 1 zero-matrix
-] unit-test
-{
- { { 1 0 0 }
- { 0 1 0 }
- { 0 0 1 } }
-} [
- 3 identity-matrix
-] unit-test
+{ t } [ { } square-matrix? ] unit-test
+{ t } [ { { 1 } } square-matrix? ] unit-test
+{ t } [ { { 1 2 } { 3 4 } } square-matrix? ] unit-test
+{ f } [ { { 1 } { 2 3 } } square-matrix? ] unit-test
+{ f } [ { { 1 2 } } square-matrix? ] unit-test
-{
- { { 1 0 0 }
- { 0 2 0 }
- { 0 0 3 } }
-} [
- { 1 2 3 } diagonal-matrix
+! any deep-empty matrix is null
+! it doesn't make any sense for { } to be null while { { } } to be considered nonnull
+{ t } [ {
+ { }
+ { { } }
+ { { { } } }
+ { { } { } { } }
+ { { { } } { { { } } } }
+} [ null-matrix? ] map [ ] all?
] unit-test
-{
- { { 1 1 1 }
- { 4 2 1 }
- { 9 3 1 }
- { 25 5 1 } }
-} [
- { 1 2 3 5 } 3 vandermonde-matrix
+{ f } [ {
+ { 1 2 }
+ { { 1 2 } }
+ { { 1 } { 2 } }
+ { { { 1 } } { 2 } { } }
+} [ null-matrix? ] map [ ] any?
] unit-test
-{
- {
- { 1 0 0 }
- { 0 1 0 }
- { 0 0 1 }
- }
-} [
- 3 3 0 eye
+{ t } [ 10 dup <zero-matrix> zero-matrix? ] unit-test
+{ t } [ 10 10 15 <simple-eye> zero-matrix? ] unit-test
+{ t } [ 0 dup <zero-matrix> null-matrix? ] unit-test
+{ f } [ 0 dup <zero-matrix> zero-matrix? ] unit-test
+{ f } [ 4 <identity-matrix> zero-matrix? ] unit-test
+! make sure we're not using the sum-to-zero strategy
+{ f } [ { { 0 -2 } { 1 -1 } } zero-matrix? ] unit-test
+{ f } [ { { 0 0 } { 1 -1 } } zero-matrix? ] unit-test
+{ f } [ { { 0 1 } { 0 -1 } } zero-matrix? ] unit-test
+
+! nth etc
+
+{ 3 } [ { 1 2 3 } 0 swap nth-end ] unit-test
+{ 2 } [ { 1 2 3 } 1 swap nth-end ] unit-test
+{ 1 } [ { 1 2 3 } 2 swap nth-end ] unit-test
+
+[ { 1 2 3 } -1 swap nth-end ] [ bounds-error? ] must-fail-with
+[ { 1 2 3 } 3 swap nth-end ] [ bounds-error? ] must-fail-with
+[ { 1 2 3 } 4 swap nth-end ] [ bounds-error? ] must-fail-with
+
+{ { 0 0 1 } } [ { 0 0 0 } dup 1 0 rot set-nth-end ] unit-test
+{ { 0 2 0 } } [ { 0 0 0 } dup 2 1 rot set-nth-end ] unit-test
+{ { 3 0 0 } } [ { 0 0 0 } dup 3 2 rot set-nth-end ] unit-test
+
+[ { 0 0 0 } dup 1 -1 rot set-nth-end ] [ bounds-error? ] must-fail-with
+[ { 0 0 0 } dup 2 3 rot set-nth-end ] [ bounds-error? ] must-fail-with
+[ { 0 0 0 } dup 3 4 rot set-nth-end ] [ bounds-error? ] must-fail-with
+
+! constructors
+
+{ {
+ { 5 5 }
+ { 5 5 }
+} } [ 2 2 5 <matrix> ] unit-test
+! a matrix-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 } }
+} } } [ 2 3 2 2 -1 <matrix> <matrix> ] unit-test
+
+{ {
+ { 5 5 }
+ { 5 5 }
+} } [ 2 2 [ 5 ] <matrix-by> ] unit-test
+{ {
+ { 6 6 }
+ { 6 6 }
+} } [ 2 2 [ 3 2 * ] <matrix-by> ] unit-test
+
+{ {
+ { 0 1 2 }
+ { 1 2 3 }
+} } [ 2 3 [ + ] <matrix-by-indices> ] unit-test
+{ {
+ { 0 0 0 }
+ { 0 1 2 }
+ { 0 2 4 }
+} } [ 3 3 [ * ] <matrix-by-indices> ] unit-test
+
+{ t } [ 3 3 <zero-matrix> zero-square-matrix? ] unit-test
+{ t } [ 3 <zero-square-matrix> zero-square-matrix? ] unit-test
+{ t f } [ 3 1 <zero-matrix> [ zero-matrix? ] [ square-matrix? ] bi ] unit-test
+
+{ {
+ { 1 0 0 }
+ { 0 2 0 }
+ { 0 0 3 }
+} } [
+ { 1 2 3 } <diagonal-matrix>
] unit-test
-{
- {
- { 0 1 0 }
- { 0 0 1 }
- { 0 0 0 }
- }
-} [
- 3 3 1 eye
+{ {
+ { -11 0 0 0 }
+ { 0 -12 0 0 }
+ { 0 0 -33 0 }
+ { 0 0 0 -14 }
+} } [ { -11 -12 -33 -14 } <diagonal-matrix> ] unit-test
+
+{ {
+ { 0 0 1 }
+ { 0 2 0 }
+ { 3 0 0 }
+} } [ { 1 2 3 } <anti-diagonal-matrix> ] unit-test
+
+{ {
+ { 0 0 0 -11 }
+ { 0 0 -12 0 }
+ { 0 -33 0 0 }
+ { -14 0 0 0 }
+} } [ { -11 -12 -33 -14 } <anti-diagonal-matrix> ] unit-test
+
+{ {
+ { 1 0 0 }
+ { 0 1 0 }
+ { 0 0 1 }
+} } [
+ 3 <identity-matrix>
] unit-test
-{
- {
- { 0 0 0 }
- { 1 0 0 }
- { 0 1 0 }
- }
-} [
- 3 3 -1 eye
+{ {
+ { 2 0 0 }
+ { 0 2 0 }
+ { 0 0 2 }
+} } [
+ 3 3 0 2 <eye>
] unit-test
-{
- {
- { 1 0 0 0 }
- { 0 1 0 0 }
- { 0 0 1 0 }
- }
-} [
- 3 4 0 eye
+{ {
+ { 0 2 0 }
+ { 0 0 2 }
+ { 0 0 0 }
+} } [
+ 3 3 1 2 <eye>
] unit-test
-{
- {
- { 0 1 0 }
- { 0 0 1 }
- { 0 0 0 }
- { 0 0 0 }
- }
-} [
- 4 3 1 eye
-] unit-test
-
-{
- {
- { 0 0 0 }
- { 1 0 0 }
- { 0 1 0 }
- { 0 0 1 }
- }
-} [
- 4 3 -1 eye
+{ {
+ { 0 0 0 0 }
+ { 2 0 0 0 }
+ { 0 2 0 0 }
+} } [
+ 3 4 -1 2 <eye>
] 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 0 4 }
- { 0 7 0 }
- { 6 0 3 } }
-} [
- { { 1 0 0 }
- { 0 2 0 }
- { 0 0 3 } }
-
- { { 0 0 4 }
- { 0 5 0 }
- { 6 0 0 } }
+{ {
+ { 1 0 0 }
+ { 0 1 0 }
+ { 0 0 1 }
+} } [
+ 3 3 0 <simple-eye>
+] unit-test
- m+
+{ {
+ { 0 1 0 }
+ { 0 0 1 }
+ { 0 0 0 }
+} } [
+ 3 3 1 <simple-eye>
] unit-test
-{
- { { 1 0 4 }
- { 0 7 0 }
- { 6 0 3 } }
-} [
- { { 1 0 0 }
- { 0 2 0 }
- { 0 0 3 } }
+{ {
+ { 0 0 0 }
+ { 1 0 0 }
+ { 0 1 0 }
+} } [
+ 3 3 -1 <simple-eye>
+] unit-test
- { { 0 0 -4 }
- { 0 -5 0 }
- { -6 0 0 } }
+{ {
+ { 1 0 0 0 }
+ { 0 1 0 0 }
+ { 0 0 1 0 }
+} } [
+ 3 4 0 <simple-eye>
+] unit-test
- m-
+{ {
+ { 0 1 0 }
+ { 0 0 1 }
+ { 0 0 0 }
+ { 0 0 0 }
+} } [
+ 4 3 1 <simple-eye>
] unit-test
-{
- { 10 20 30 }
-} [
- 10 { 1 2 3 } n*v
+{ {
+ { 0 0 0 }
+ { 1 0 0 }
+ { 0 1 0 }
+ { 0 0 1 }
+} } [
+ 4 3 -1 <simple-eye>
] unit-test
-{
+{ {
+ { { 0 0 } { 0 1 } { 0 2 } }
+ { { 1 0 } { 1 1 } { 1 2 } }
+ { { 2 0 } { 2 1 } { 2 2 } }
+ { { 3 0 } { 3 1 } { 3 2 } }
+} } [ { 4 3 } <coordinate-matrix> ] unit-test
+
+{ {
+ { 0 1 }
+ { 0 1 }
+} } [ 2 <square-rows> ] unit-test
+
+{ {
+ { 0 0 }
+ { 1 1 }
+} } [ 2 <square-cols> ] unit-test
+
+{ {
+ { 5 6 }
+ { 5 6 }
+} } [ { 5 6 } <square-rows> ] unit-test
+
+{ {
+ { 5 5 }
+ { 6 6 }
+} } [ { 5 6 } <square-cols> ] unit-test
+
+{ {
+ { 1 }
+} } [ {
+ { 1 2 }
+} <square-rows> ] unit-test
+
+{ {
+ { 1 2 }
{ 3 4 }
-} [
- { { 1 0 }
- { 0 1 } }
+} } [ {
+ { 1 2 5 }
+ { 3 4 6 }
+} <square-rows> ] unit-test
+{ {
+ { 1 2 }
{ 3 4 }
-
- m.v
-] unit-test
-
-{
- { 4 3 }
-} [
- { { 0 1 }
- { 1 0 } }
-
+} } [ {
+ { 1 2 }
{ 3 4 }
-
- m.v
-] unit-test
-
-{
- { { 6 } }
-} [
- { { 3 } } { { 2 } } m.
+ { 5 6 }
+} <square-rows> ] unit-test
+
+{ {
+ { 1 0 4 }
+ { 0 7 0 }
+ { 6 0 3 } }
+} [ {
+ { 1 0 0 }
+ { 0 2 0 }
+ { 0 0 3 }
+} {
+ { 0 0 4 }
+ { 0 5 0 }
+ { 6 0 0 }
+}
+ m+
] unit-test
-{
- { { 11 } }
-} [
- { { 1 3 } } { { 5 } { 2 } } m.
+{ {
+ { 1 0 4 }
+ { 0 7 0 }
+ { 6 0 3 }
+} } [ {
+ { 1 0 0 }
+ { 0 2 0 }
+ { 0 0 3 }
+} {
+ { 0 0 -4 }
+ { 0 -5 0 }
+ { -6 0 0 }
+}
+ m-
] unit-test
-{
- { { 28 } }
-} [
- { { 2 4 6 } }
+{ { 3 4 } } [ { { 1 0 } { 0 1 } } { 3 4 } m.v ] unit-test
+{ { 4 3 } } [ { { 0 1 } { 1 0 } } { 3 4 } m.v ] unit-test
- { { 1 }
- { 2 }
- { 3 } }
+{ { { 6 } } } [ { { 3 } } { { 2 } } m. ] unit-test
+{ { { 11 } } } [ { { 1 3 } } { { 5 } { 2 } } m. ] unit-test
+{ { { 28 } } } [
+ { { 2 4 6 } }
+ { { 1 } { 2 } { 3 } }
m.
] 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
-
-{ { { 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 5 0 10 } { 6 7 12 14 } { 0 15 0 20 } { 18 21 24 28 } }
-}
-[ { { 1 2 } { 3 4 } } { { 0 5 } { 6 7 } } kron ] unit-test
+! TODO: note: merge conflict from HEAD contained the following
+! ------------------------
+! predicates
-{
- {
- { 1 1 1 1 }
- { 1 -1 1 -1 }
- { 1 1 -1 -1 }
- { 1 -1 -1 1 }
- }
-} [ { { 1 1 } { 1 -1 } } dup kron ] 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 kron kron ] unit-test
+{ t } [ { } square-matrix? ] unit-test
+{ t } [ { { 1 } } square-matrix? ] unit-test
+{ t } [ { { 1 2 } { 3 4 } } square-matrix? ] unit-test
+{ f } [ { { 1 } { 2 3 } } square-matrix? ] unit-test
+{ f } [ { { 1 2 } } square-matrix? ] 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 kron swap kron ] unit-test
+{ 9 }
+[ { { 2 -2 1 } { 1 3 -1 } { 2 -4 2 } } m-1norm ] unit-test
+{ 8 }
+[ { { 2 -2 1 } { 1 3 -1 } { 2 -4 2 } } m-infinity-norm ] unit-test
-! kron 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 } } kron
+{ 2.0 }
+[ { { 1 1 } { 1 1 } } frobenius-norm ] unit-test
+! from "intermediate commit"
+! any deep-empty matrix is null
+! it doesn't make any sense for { } to be null while { { } } to be considered nonnull
+{ t } [ {
+ { }
+ { { } }
+ { { { } } }
+ { { } { } { } }
+ { { { } } { { { } } } }
+} [ null-matrix? ] map [ ] all?
] 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 kron
+{ f } [ {
+ { 1 2 }
+ { { 1 2 } }
+ { { 1 } { 2 } }
+ { { { 1 } } { 2 } { } }
+} [ null-matrix? ] map [ ] any?
] unit-test
-{
- { { 5 10 15 }
- { 6 12 18 }
- { 7 14 21 } }
-} [ { 5 6 7 } { 1 2 3 } outer ] 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-cov-matrix
+{ t } [ 10 dup <zero-matrix> zero-matrix? ] unit-test
+{ t } [ 10 10 15 <simple-eye> zero-matrix? ] unit-test
+{ t } [ 0 dup <zero-matrix> null-matrix? ] unit-test
+{ f } [ 0 dup <zero-matrix> zero-matrix? ] unit-test
+{ f } [ 4 <identity-matrix> zero-matrix? ] unit-test
+
+{ t } [ { } regular-matrix? ] unit-test
+{ t } [ { { } } regular-matrix? ] unit-test
+{ t } [ { { 1 2 } } regular-matrix? ] unit-test
+{ t } [ { { 1 2 } { 3 4 } } regular-matrix? ] unit-test
+{ t } [ { { 1 } { 3 } } regular-matrix? ] unit-test
+{ f } [ { { 1 2 } { 3 } } regular-matrix? ] unit-test
+{ f } [ { { 1 } { 3 2 } } regular-matrix? ] unit-test
+! TODO: note: lines since last HEAD comment were deleted in "fix more code and add more rigorous tests"
+
+! diagonals
+
+! diagonal getters
+{ { 1 1 1 1 } } [ 4 <identity-matrix> main-diagonal ] unit-test
+{ { 0 0 0 0 } } [ 4 <identity-matrix> anti-diagonal ] unit-test
+{ { 4 8 } } [ { { 4 6 } { 3 8 } } main-diagonal ] unit-test
+{ { 6 3 } } [ { { 4 6 } { 3 8 } } anti-diagonal ] unit-test
+{ { 1 2 3 } } [ { { 0 0 1 } { 0 2 0 } { 3 0 0 } } anti-diagonal ] unit-test
+{ { 1 2 3 4 } } [ { 1 2 3 4 } <diagonal-matrix> main-diagonal ] unit-test
+
+! transposition
+{ { 1 2 3 4 } } [ { 1 2 3 4 } <diagonal-matrix> transpose main-diagonal ] unit-test
+{ t } [ 50 <identity-matrix> dup transpose = ] unit-test
+{ { 4 3 2 1 } } [ { 1 2 3 4 } <anti-diagonal-matrix> transpose anti-diagonal ] unit-test
+
+{ {
+ { 1 4 7 }
+ { 2 5 8 }
+ { 3 6 9 }
+} } [ {
+ { 1 2 3 }
+ { 4 5 6 }
+ { 7 8 9 }
+} transpose ] unit-test
+
+! anti transposition
+{ { 1 2 3 4 } } [ { 1 2 3 4 } <anti-diagonal-matrix> anti-transpose anti-diagonal ] unit-test
+{ t } [ 50 <iota> <anti-diagonal-matrix> dup anti-transpose = ] unit-test
+{ { 4 3 2 1 } } [ { 1 2 3 4 } <diagonal-matrix> anti-transpose main-diagonal ] unit-test
+
+{ {
+ { 9 6 3 }
+ { 8 5 2 }
+ { 7 4 1 }
+} } [ {
+ { 1 2 3 }
+ { 4 5 6 }
+ { 7 8 9 }
+} anti-transpose ] unit-test
+
+<PRIVATE
+SYMBOLS: A B C D E F G H I J K L M N O P ;
+PRIVATE>
+{ { {
+ { E F G H }
+ { I J K L }
+ { M N O P }
+} {
+ { A B C D }
+ { I J K L }
+ { M N O P }
+} {
+ { A B C D }
+ { E F G H }
+ { M N O P }
+} {
+ { A B C D }
+ { E F G H }
+ { I J K L }
+} } } [
+ 4 {
+ { A B C D }
+ { E F G H }
+ { I J K L }
+ { M N O P }
+ } <repetition>
+ [ rows-except ] map-index
] 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-cov-matrix
+{ { { 2 } } } [ { { 1 } { 2 } } 0 rows-except ] unit-test
+{ { { 1 } } } [ { { 1 } { 2 } } 1 rows-except ] unit-test
+{ { } } [ { { 1 } } 0 rows-except ] unit-test
+{ { { 1 } } } [ { { 1 } } 1 rows-except ] unit-test
+{ {
+ { 2 7 12 2 } ! 0
+ { 1 3 3 5 } ! 2
+} } [ {
+ { 2 7 12 2 }
+ { 8 9 10 0 }
+ { 1 3 3 5 }
+ { 8 13 7 12 }
+} { 1 3 } rows-except ] unit-test
+
+{ { {
+ { B C D }
+ { F G H }
+ { J K L }
+ { N O P }
+} {
+ { A C D }
+ { E G H }
+ { I K L }
+ { M O P }
+} {
+ { A B D }
+ { E F H }
+ { I J L }
+ { M N P }
+} {
+ { A B C }
+ { E F G }
+ { I J K }
+ { M N O }
+} } } [
+ 4 {
+ { A B C D }
+ { E F G H }
+ { I J K L }
+ { M N O P }
+ } <repetition>
+ [ cols-except ] map-index
] unit-test
-{
- {
- { 5 5 }
- { 5 5 }
- }
-} [
- 2 2 5 <matrix>
+{ { } } [ { { 1 } { 2 } } 0 cols-except ] unit-test
+{ { { 1 } { 2 } } } [ { { 1 } { 2 } } 1 cols-except ] unit-test
+{ { } } [ { { 1 } } 0 cols-except ] unit-test
+{ { { 1 } } } [ { { 1 } } 1 cols-except ] unit-test
+{ { { 2 } { 4 } } } [ { { 1 2 } { 3 4 } } 0 cols-except ] unit-test
+{ { { 1 } { 3 } } } [ { { 1 2 } { 3 4 } } 1 cols-except ] unit-test
+{ {
+ { 2 12 }
+ { 8 10 }
+ { 1 3 }
+ { 8 7 }
+} } [ {
+ { 2 7 12 2 }
+ { 8 9 10 0 }
+ { 1 3 3 5 }
+ { 8 13 7 12 }
+} { 1 3 } cols-except ] unit-test
+
+{ { {
+ { F G H }
+ { J K L }
+ { N O P }
+} {
+ { A C D }
+ { I K L }
+ { M O P }
+} {
+ { A B D }
+ { E F H }
+ { M N P }
+} {
+ { A B C }
+ { E F G }
+ { I J K }
+} } } [
+ 4 {
+ { A B C D }
+ { E F G H }
+ { I J K L }
+ { M N O P }
+ } <repetition>
+ [ dup 2array matrix-except ] map-index
] unit-test
-{
- {
- { 5 5 }
- { 5 5 }
- }
-} [
- 2 2 [ 5 ] make-matrix
-] unit-test
+! prepare for bracket hell
+! going to test the Matrix of Minors permutation strategy
-{
- {
- { 0 1 2 }
- { 1 2 3 }
- }
-} [
- 2 3 [ + ] make-matrix-with-indices
-] unit-test
+! going to test 1x2 inputs
+! the input had 2 elements, the output has 2 0-matrices across 2 arrays ;)
+{ { { { } { } } } } [ { { 1 2 } } matrix-except-all ] unit-test
-{
- {
- { 0 1 }
- { 0 1 }
- }
-} [
- 2 square-rows
-] unit-test
+! any matrix with a 1 in its dimensions will give a void matrix output
+{ t } [ { { 1 2 } } matrix-except-all null-matrix? ] unit-test
+{ t } [ { { 1 } { 2 } } matrix-except-all null-matrix? ] unit-test
-{
- {
- { 0 0 }
- { 1 1 }
- }
-} [
- 2 square-cols
-] unit-test
+! going to test 2x2 inputs
+! these 1x1 output matrices have omitted a row and column from the 2x2 input
+! the input had 4 elements, the output has 4 1-matrices across 2 arrays
+! the permutations of indices 0 1 are: 0 0, 0 1, 1 0, 1 1
{
- {
- { 5 6 }
- { 5 6 }
+ { ! output array
+ { ! item #1: excluding row 0...
+ { { 3 } } ! and col 0 = 0 0
+ { { 2 } } ! and col 1 = 0 1
+ }
+ { ! item #2: excluding row 1...
+ { { 1 } } ! and col 0 = 1 0
+ { { 0 } } ! and col 1 = 1 1
+ }
}
} [
- { 5 6 } square-rows
+ ! the input to the function is a simple 2x2
+ { { 0 1 } { 2 3 } } matrix-except-all
] unit-test
+! we are going to ensure that "duplicate" matrices are not omitted in the output
{
{
- { 5 5 }
- { 6 6 }
+ { ! item 1
+ { { 0 } }
+ { { 0 } }
+ }
+ { ! item 2
+ { { 0 } }
+ { { 0 } }
+ }
}
-} [
- { 5 6 } square-cols
+} [ { { 0 0 } { 0 0 } } matrix-except-all ] unit-test
+! the output only has elements from the input
+{ t } [ 44 <zero-square-matrix> matrix-except-all zero-matrix? ] unit-test
+
+! going to test 2x3 and 3x2 inputs
+{
+ { ! output array
+ { ! excluding row 0
+ { { 2 } { 3 } } ! and col 0
+ { { 1 } { 2 } } ! and col 1
+ }
+ { ! excluding row 1
+ { { 1 } { 3 } } ! and col 0
+ { { 0 } { 2 } } ! and col 1
+ }
+ { ! excluding row 2
+ { { 1 } { 2 } } ! col 0
+ { { 0 } { 1 } } ! col 1
+ }
+ }
+} [ {
+ { 0 1 }
+ { 1 2 }
+ { 2 3 }
+} matrix-except-all ] unit-test
+
+{
+ { ! output array
+ { ! excluding row 0
+ { { 2 3 } } ! col 0
+ { { 1 3 } } ! col 1
+ { { 1 2 } } ! col 2
+ }
+ { ! row 1
+ { { 1 2 } } ! col 0
+ { { 0 2 } } ! col 1
+ { { 0 1 } } ! col 2
+ }
+ }
+} [ {
+ { 0 1 2 }
+ { 1 2 3 }
+} matrix-except-all ] unit-test
+
+! going to test 3x3 inputs
+
+! the input had 9 elements, the output has 9 2-matrices across 3 arrays
+! every element from the input is represented 4 times in the output
+! the number of copies of each element found in the output is the side length of the next smaller square matrix
+! 3x3 input gives 4 copies of each element; (N-1) ^ 2 = 4 where N=3
+! the permutations of indices 0 1 2 are: 0 0, 0 1, 0 2; 1 0, 1 1, 1 2; 2 0, 2 1, 2 2
+{
+ { ! output array
+ { ! item #1: excluding row 0...
+ { ! and col 0 = 0 0
+ { 4 5 }
+ { 7 8 }
+ }
+ { ! and col 1 = 0 1
+ { 3 5 }
+ { 6 8 }
+ }
+ { ! and col 2 = 0 2
+ { 3 4 }
+ { 6 7 }
+ }
+ }
+
+ { ! item #2: excluding row 1...
+ { ! and col 0 = 1 0
+ { 1 2 }
+ { 7 8 }
+ }
+ { ! and col 1 = 1 1
+ { 0 2 }
+ { 6 8 }
+ }
+ { ! and col 2 = 1 2
+ { 0 1 }
+ { 6 7 }
+ }
+ }
+
+ { ! item #2: excluding row 2...
+ { ! and col 0 = 2 0
+ { 1 2 }
+ { 4 5 }
+ }
+ { ! and col 1 = 2 1
+ { 0 2 }
+ { 3 5 }
+ }
+ { ! and col 2 = 2 2
+ { 0 1 }
+ { 3 4 }
+ }
+ }
+ }
+ t ! note this
+} [ {
+ { 0 1 2 }
+ { 3 4 5 }
+ { 6 7 8 }
+} matrix-except-all dup flatten sorted-histogram [ second ] map
+ { [ length 9 = ] [ [ 4 = ] all? ] }
+ 1&&
] unit-test
-{ t } [ { } square-matrix? ] unit-test
-{ t } [ { { 1 } } square-matrix? ] unit-test
-{ t } [ { { 1 2 } { 3 4 } } square-matrix? ] unit-test
-{ f } [ { { 1 } { 2 3 } } square-matrix? ] unit-test
-{ f } [ { { 1 2 } } square-matrix? ] unit-test
-
-{ 9 }
-[ { { 2 -2 1 } { 1 3 -1 } { 2 -4 2 } } m-1norm ] unit-test
-
-{ 8 }
-[ { { 2 -2 1 } { 1 3 -1 } { 2 -4 2 } } m-infinity-norm ] unit-test
-
-{ 2.0 }
-[ { { 1 1 } { 1 1 } } frobenius-norm ] unit-test
+! going to test 4x4 inputs
+
+! don't feel like handwriting this right now, so a sanity check test instead
+! the input contains 4 rows and 4 columns for 16 elements
+! 4x4 input gives 9 copies of each element; (N-1) ^ 2 = 9 where N = 4
+{ t } [ {
+ { 0 1 2 3 }
+ { 4 5 6 7 }
+ { 8 9 10 11 }
+ { 12 13 14 15 }
+} matrix-except-all flatten sorted-histogram [ second ] map
+ { [ length 16 = ] [ [ 9 = ] all? ] }
+ 1&&
+] unit-test
-! 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.
-USING: accessors arrays columns kernel locals math math.bits
-math.functions math.order math.vectors sequences
-sequences.private fry math.statistics grouping
-combinators.short-circuit math.ranges combinators.smart ;
+USING: accessors arrays classes.singleton columns combinators
+combinators.short-circuit combinators.smart formatting fry
+grouping kernel locals math math.bits math.functions math.order
+math.private math.ranges math.statistics math.vectors
+math.vectors.private sequences sequences.deep sequences.private
+slots.private summary ;
IN: math.matrices
-! Matrices
-: make-matrix ( m n quot -- matrix )
- '[ _ _ replicate ] replicate ; inline
+! defined here because of issue #1943
+DEFER: regular-matrix?
+: 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 )
'[ _ _ <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
-: diagonal-matrix ( diagonal-seq -- matrix )
- dup length dup zero-matrix
- [ '[ dup _ nth set-nth ] each-index ] keep ; inline
-
-: identity-matrix ( n -- matrix )
- 1 <repetition> diagonal-matrix ; inline
-
-: eye ( m n k -- matrix )
- [ [ <iota> ] bi@ ] dip neg '[ _ + = 1 0 ? ] cartesian-map ;
-
-: hilbert-matrix ( m n -- matrix )
- [ <iota> ] bi@ [ + 1 + recip ] cartesian-map ;
-
-: toeplitz-matrix ( n -- matrix )
- <iota> dup [ - abs 1 + ] cartesian-map ;
-
-: hankel-matrix ( n -- matrix )
- [ <iota> dup ] keep '[ + abs 1 + dup _ > [ drop 0 ] when ] cartesian-map ;
-
-: box-matrix ( r -- matrix )
- 2 * 1 + dup '[ _ 1 <array> ] replicate ;
-
-: vandermonde-matrix ( u n -- matrix )
- <iota> [ v^n ] with map reverse flip ;
-
-:: 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 }
- } ;
-
-: >scale-factors ( number/sequence -- x y z )
- dup number? [ dup dup ] [ first3 ] if ;
-
-:: 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 ( dim -- 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 }
- } ;
-
-! Matrix operations
-: mneg ( m -- m ) [ vneg ] map ;
-
-: n+m ( n m -- m ) [ n+v ] with map ;
-: m+n ( m n -- m ) [ v+n ] curry map ;
-: n-m ( n m -- m ) [ n-v ] with map ;
-: m-n ( m n -- m ) [ v-n ] curry map ;
-: n*m ( n m -- m ) [ n*v ] with map ;
-: m*n ( m n -- m ) [ v*n ] curry map ;
-: n/m ( n m -- m ) [ n/v ] with map ;
-: m/n ( m n -- m ) [ v/n ] curry map ;
+: <zero-square-matrix> ( n -- matrix )
+ dup <zero-matrix> ; inline
-: m+ ( m m -- m ) [ v+ ] 2map ;
-: m- ( m m -- m ) [ v- ] 2map ;
-: m* ( m m -- m ) [ v* ] 2map ;
-: m/ ( m m -- m ) [ v/ ] 2map ;
+<PRIVATE
+: (nth-from-end) ( n seq -- n )
+ length 1 - swap - ; inline flushable
-: v.m ( v m -- v ) flip [ v. ] with map ;
-: m.v ( m v -- v ) [ v. ] curry map ;
-: m. ( m m -- m ) flip [ swap m.v ] curry map ;
+: nth-end ( n seq -- elt )
+ [ (nth-from-end) ] keep nth ; inline flushable
-: m~ ( m m epsilon -- ? ) [ v~ ] curry 2all? ;
+: nth-end-unsafe ( n seq -- elt )
+ [ (nth-from-end) ] keep nth-unsafe ; inline flushable
-: mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ;
-: mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ;
-: mnorm ( m -- n ) dup mmax abs m/n ;
-: m-infinity-norm ( m -- n ) [ [ abs ] map-sum ] map supremum ;
-: m-1norm ( m -- n ) flip m-infinity-norm ;
-: frobenius-norm ( m -- n ) [ [ sq ] map-sum ] map-sum sqrt ;
+: array-nth-end-unsafe ( n seq -- elt )
+ [ (nth-from-end) ] keep swap 2 fixnum+fast slot ; inline flushable
-: 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
+: set-nth-end ( elt n seq -- )
+ [ (nth-from-end) ] keep set-nth ; inline
-:: normal ( vec1 vec2 vec3 -- vec4 )
- vec2 vec1 v- vec3 vec1 v- cross normalize ; inline
+: set-nth-end-unsafe ( elt n seq -- )
+ [ (nth-from-end) ] keep set-nth-unsafe ; inline
+PRIVATE>
-: proj ( v u -- w )
- [ [ v. ] [ norm-sq ] bi / ] keep n*v ;
+! main-diagonal matrix
+: <diagonal-matrix> ( diagonal-seq -- matrix )
+ [ length <zero-square-matrix> ] keep over
+ '[ dup _ nth set-nth-unsafe ] each-index ; inline
-: perp ( v u -- w )
- dupd proj v- ;
+! could also be written slower as: <diagonal-matrix> [ reverse ] map
+: <anti-diagonal-matrix> ( diagonal-seq -- matrix )
+ [ length <zero-square-matrix> ] keep over
+ '[ dup _ nth set-nth-end-unsafe ] each-index ; inline
-: angle-between ( v u -- a )
- [ normalize ] bi@ h. acos ;
+: <identity-matrix> ( n -- matrix )
+ 1 <repetition> <diagonal-matrix> ; inline
-: (gram-schmidt) ( v seq -- newseq )
- [ dupd proj v- ] each ;
+: <eye> ( m n k z -- matrix )
+ [ [ <iota> ] bi@ ] 2dip
+ '[ _ neg + = _ 0 ? ]
+ cartesian-map ; inline
-: gram-schmidt ( seq -- orthogonal )
- V{ } clone [ over (gram-schmidt) suffix! ] reduce ;
+! if m = n and k = 0 then <identity-matrix> is (possibly) more efficient
+:: <simple-eye> ( m n k -- matrix )
+ m n = k 0 = and
+ [ n <identity-matrix> ]
+ [ m n k 1 <eye> ] if ; inline
-: norm-gram-schmidt ( seq -- orthonormal )
- gram-schmidt [ normalize ] map ;
+: <coordinate-matrix> ( dim -- coordinates )
+ first2 [ <iota> ] bi@ cartesian-product ; inline
-ERROR: negative-power-matrix m n ;
+ALIAS: <cartesian-indices> <coordinate-matrix>
-: (m^n) ( m n -- n )
- make-bits over first length identity-matrix
- [ [ dupd m. ] when [ dup m. ] dip ] reduce nip ;
+: <cartesian-square-indices> ( n -- matrix )
+ dup 2array <cartesian-indices> ; inline
-: m^n ( m n -- n )
- dup 0 >= [ (m^n) ] [ negative-power-matrix ] if ;
+ALIAS: transpose flip
-: stitch ( m -- m' )
- [ ] [ [ append ] 2map ] map-reduce ;
+<PRIVATE
+: array-matrix? ( matrix -- ? )
+ [ array? ]
+ [ [ array? ] all? ] bi and ; inline foldable flushable
+
+: matrix-cols-iota ( matrix -- cols-iota )
+ first-unsafe length <iota> ; inline
-: kron ( m1 m2 -- m )
- '[ [ _ n*m ] map ] map stitch stitch ;
+: unshaped-cols-iota ( matrix -- cols-iota )
+ [ first-unsafe length 1 ] keep
+ [ length min ] (each) (each-integer) <iota> ; inline
-: outer ( u v -- m )
- [ n*v ] curry map ;
+: generic-anti-transpose-unsafe ( cols-iota matrix -- newmatrix )
+ [ <reversed> [ nth-end-unsafe ] with { } map-as ] curry { } map-as ; inline
-: row ( n matrix -- col )
+: 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 -- cols )
+: rows ( seq matrix -- rows )
'[ _ row ] map ; inline
: col ( n matrix -- col )
: cols ( seq matrix -- cols )
'[ _ col ] map ; inline
-: set-index ( object pair matrix -- )
- [ first2 swap ] dip nth set-nth ; 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 ;
-: set-indices ( object sequence matrix -- )
- '[ _ set-index ] with each ; inline
+<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 ;
-: matrix-map ( matrix quot -- )
+! 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 -- seq )
- [ [ first length <iota> ] keep ] dip '[ _ col @ ] map ; inline
+: column-map ( matrix quot: ( ... col -- ... col' ) -- matrix' )
+ [ transpose ] dip map transpose ; inline
-: cartesian-square-indices ( n -- matrix )
- <iota> dup cartesian-product ; inline
+: matrix-nth ( pair matrix -- elt )
+ [ first2 swap ] dip nth nth ; inline
-: cartesian-matrix-map ( matrix quot -- matrix' )
- [ [ first length cartesian-square-indices ] keep ] dip
- '[ _ @ ] matrix-map ; inline
+: matrix-nths ( pairs matrix -- elts )
+ '[ _ matrix-nth ] map ; inline
-: cartesian-matrix-column-map ( matrix quot -- matrix' )
- [ cols first2 ] prepose cartesian-matrix-map ; inline
+: matrix-set-nth ( obj pair matrix -- )
+ [ first2 swap ] dip nth set-nth ; inline
-: cov-matrix-ddof ( matrix ddof -- cov )
- '[ _ cov-ddof ] cartesian-matrix-column-map ; inline
+: matrix-set-nths ( obj pairs matrix -- )
+ '[ _ matrix-set-nth ] with each ; inline
-: population-cov-matrix ( matrix -- cov ) 0 cov-matrix-ddof ; inline
+! -------------------------------------------
+! simple math of matrices follows
+: mneg ( m -- m' ) [ vneg ] map ;
+: mabs ( m -- m' ) [ vabs ] map ;
-: sample-cov-matrix ( matrix -- cov ) 1 cov-matrix-ddof ; inline
+: n+m ( n m -- m ) [ n+v ] with map ;
+: m+n ( m n -- m ) [ v+n ] curry map ;
+: n-m ( n m -- m ) [ n-v ] with map ;
+: m-n ( m n -- m ) [ v-n ] curry map ;
+: n*m ( n m -- m ) [ n*v ] with map ;
+: m*n ( m n -- m ) [ v*n ] curry map ;
+: n/m ( n m -- m ) [ n/v ] with map ;
+: m/n ( m n -- m ) [ v/n ] curry map ;
-GENERIC: square-rows ( object -- matrix )
-M: integer square-rows <iota> square-rows ;
-M: sequence square-rows
- [ length ] keep >array '[ _ clone ] { } replicate-as ;
+: m+ ( m1 m2 -- m ) [ v+ ] 2map ;
+: m- ( m1 m2 -- m ) [ v- ] 2map ;
+: m* ( m1 m2 -- m ) [ v* ] 2map ;
+: m/ ( m1 m2 -- m ) [ v/ ] 2map ;
-GENERIC: square-cols ( object -- matrix )
-M: integer square-cols <iota> square-cols ;
-M: sequence square-cols
- [ length ] keep [ <array> ] with { } map-as ;
+: v.m ( v m -- p ) flip [ v. ] with map ;
+: m.v ( m v -- p ) [ v. ] curry map ;
+: m. ( m m -- m ) flip [ swap m.v ] curry map ;
-: make-matrix-with-indices ( m n quot -- matrix )
- [ [ <iota> ] bi@ ] dip cartesian-map ; inline
+: m~ ( m1 m2 epsilon -- ? ) [ v~ ] curry 2all? ;
-: null-matrix? ( matrix -- ? ) empty? ; inline
+: mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ;
+: mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ;
+: mnorm ( m -- m' ) dup mmax abs m/n ;
+: m-infinity-norm ( m -- n ) [ [ abs ] map-sum ] map supremum ;
+: m-1norm ( m -- n ) flip m-infinity-norm ;
+: frobenius-norm ( m -- n ) [ [ sq ] map-sum ] map-sum sqrt ;
-: well-formed-matrix? ( matrix -- ? )
- [ t ] [
- [ ] [ first length ] bi
- '[ length _ = ] all?
- ] if-empty ;
+! well-defined for square matrices; but works on nonsquare too
+: main-diagonal ( matrix -- seq )
+ >square-matrix [ swap nth-unsafe ] map-index ; inline
-: dim ( matrix -- pair/f )
- [ 2 0 <array> ]
- [ [ length ] [ first length ] bi 2array ] if-empty ;
+! top right to bottom left; reverse the result if you expected it to start in the lower left
+: anti-diagonal ( matrix -- seq )
+ >square-matrix [ swap nth-end-unsafe ] map-index ; inline
-: square-matrix? ( matrix -- ? )
- { [ well-formed-matrix? ] [ dim all-eq? ] } 1&& ;
+<PRIVATE
+: (rows-iota) ( matrix -- rows-iota )
+ dimension first-unsafe <iota> ;
+: (cols-iota) ( matrix -- cols-iota )
+ dimension second-unsafe <iota> ;
-: matrix-coordinates ( dim -- coordinates )
- first2 [ <iota> ] bi@ cartesian-product ; inline
+: simple-rows-except ( matrix desc quot -- others )
+ curry [ dup (rows-iota) ] dip
+ pick reject-as swap rows ; inline
-: dimension-range ( matrix -- dim range )
- dim [ matrix-coordinates ] [ first [1,b] ] bi ;
+: simple-cols-except ( matrix desc quot -- others )
+ curry [ dup (cols-iota) ] dip
+ pick reject-as swap cols transpose ; inline ! need to un-transpose the result of cols
-: upper-matrix-indices ( matrix -- matrix' )
- dimension-range <reversed> [ tail-slice* >array ] 2map concat ;
+CONSTANT: scalar-except-quot [ = ]
+CONSTANT: sequence-except-quot [ member? ]
+PRIVATE>
-: lower-matrix-indices ( matrix -- matrix' )
- dimension-range [ head-slice >array ] 2map concat ;
+GENERIC: rows-except ( matrix desc -- others )
+M: integer rows-except scalar-except-quot simple-rows-except ;
+M: sequence rows-except sequence-except-quot simple-rows-except ;
+
+GENERIC: cols-except ( matrix desc -- others )
+M: integer cols-except scalar-except-quot simple-cols-except ;
+M: sequence cols-except sequence-except-quot simple-cols-except ;
+
+! well-defined for any regular matrix
+: matrix-except ( matrix exclude-pair -- submatrix )
+ first2 [ rows-except ] dip cols-except ;
+
+ALIAS: submatrix-excluding matrix-except
+
+:: matrix-except-all ( matrix -- submatrices )
+ matrix dimension [ <iota> ] map first2-unsafe cartesian-product
+ [ [ matrix swap matrix-except ] map ] map ;
-: make-lower-matrix ( object m n -- matrix )
- zero-matrix [ lower-matrix-indices ] [ set-indices ] [ ] tri ;
+ALIAS: all-submatrices matrix-except-all
-: make-upper-matrix ( object m n -- matrix )
- zero-matrix [ upper-matrix-indices ] [ set-indices ] [ ] tri ;
+: dimension ( matrix -- dimension )
+ [ { 0 0 } ]
+ [ [ length ] [ first-unsafe length ] bi 2array ] if-empty ;
} ;
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." } ;
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." } ;
HELP: vs+
kernel.private alien.c-types math.functions ;
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 } } [ { 2 4 6 } 1/2 v*n ] unit-test
{ { 1 2 3 } } [ { 2 4 6 } 2 v/n ] unit-test
{ { 0 30 100 } } [
{ -10 30 120 } { 0 0 0 } { 100 100 100 } vclamp
] 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
: vclamp ( v min max -- w )
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 ;
-USING: kernel locals math math.matrices math.order math.vectors
-prettyprint sequences ;
+USING: kernel locals math math.matrices math.matrices.extras
+math.order math.vectors prettyprint sequences ;
IN: benchmark.3d-matrix-scalar
:: p-matrix ( dim fov near far -- matrix )
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 )
- { 1.0 0.0 0.0 } pitch rotation-matrix4
- { 0.0 1.0 0.0 } yaw rotation-matrix4
- location vneg translation-matrix4 m. m. ;
+ { 1.0 0.0 0.0 } pitch <rotation-matrix4>
+ { 0.0 1.0 0.0 } yaw <rotation-matrix4>
+ location vneg <translation-matrix4> m. m. ;
:: 3d-matrix-scalar-benchmark ( -- )
f :> result!
-USING: kernel locals math math.matrices.simd math.order math.vectors
-math.vectors.simd prettyprint sequences typed ;
+USING: kernel locals math math.matrices math.matrices.simd
+math.order math.vectors math.vectors.simd prettyprint sequences
+typed ;
QUALIFIED-WITH: alien.c-types c
IN: benchmark.3d-matrix-vector
-USING: locals math math.combinatorics math.matrices
+USING: locals math math.combinatorics math.matrices math.matrices.extras
prettyprint sequences typed ;
IN: benchmark.matrix-exponential-scalar
:: matrix-exponential-scalar-benchmark ( -- )
f :> result!
- 4 identity-matrix :> i4
+ 4 <identity-matrix> :> i4
10000 [
i4 20 e^m result!
] times
! See http://factorcode.org/license.txt for BSD license.
USING: accessors colors.constants game.debug game.loop
game.worlds gpu gpu.framebuffers gpu.util.wasd kernel literals
-locals make math math.matrices math.parser math.trig sequences
-specialized-arrays ui.gadgets.worlds ui.pixel-formats ;
+locals make math math.matrices math.matrices.extras math.parser
+math.trig sequences specialized-arrays ui.gadgets.worlds
+ui.pixel-formats ;
FROM: alien.c-types => float ;
SPECIALIZED-ARRAY: float
IN: game.debug.tests
{ 0 0 0 } { 1 0 0 } COLOR: red debug-line
{ 0 0 0 } { 0 1 0 } COLOR: green 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
{ 0 9 0 } 8 2 COLOR: blue debug-cylinder
] float-array{ } make
! Copyright (C) 2010 Slava Pestov.
-USING: arrays kernel math.matrices math.vectors.simd.cords
-math.trig gml.runtime ;
+USING: arrays gml.runtime kernel math.matrices
+math.matrices.extras math.trig math.vectors.simd.cords ;
IN: gml.geometry
GML: rot_vec ( v n alpha -- v )
! Inefficient!
- deg>rad rotation-matrix4 swap >array m.v >double-4 ;
+ deg>rad <rotation-matrix4> swap >array m.v >double-4 ;
! See http://factorcode.org/license.txt for BSD license.
USING: accessors arrays combinators.tuple game.loop game.worlds
generalizations gpu gpu.render gpu.shaders gpu.util gpu.util.wasd
-kernel literals math math.libm math.matrices math.order math.vectors
-method-chains sequences ui ui.gadgets ui.gadgets.worlds
-ui.pixel-formats audio.engine audio.loader locals ;
+kernel literals math math.libm math.matrices math.matrices.extras
+math.order math.vectors method-chains sequences ui ui.gadgets
+ui.gadgets.worlds ui.pixel-formats audio.engine audio.loader locals ;
IN: gpu.demos.raytrace
GLSL-SHADER-FILE: raytrace-vertex-shader vertex-shader "raytrace.v.glsl"
dup dtheta>> [ + ] curry change-theta drop ;
: sphere-center ( sphere -- center )
- [ [ axis>> ] [ theta>> ] bi rotation-matrix4 ]
+ [ [ axis>> ] [ theta>> ] bi <rotation-matrix4> ]
[ home>> ] bi m.v ;
M: sphere audio-position sphere-center ; inline
game.input.scancodes game.loop game.worlds
gpu.render gpu.state kernel literals
locals math math.constants math.functions math.matrices
-math.order math.vectors opengl.gl sequences
-ui ui.gadgets.worlds specialized-arrays audio.engine ;
+math.matrices.extras math.order math.vectors opengl.gl
+sequences ui ui.gadgets.worlds specialized-arrays audio.engine ;
FROM: alien.c-types => float ;
SPECIALIZED-ARRAY: float
IN: gpu.util.wasd
M: wasd-world wasd-fly-vertically? drop t ;
: wasd-mv-matrix ( world -- matrix )
- [ { 1.0 0.0 0.0 } swap pitch>> rotation-matrix4 ]
- [ { 0.0 1.0 0.0 } swap yaw>> rotation-matrix4 ]
- [ location>> vneg translation-matrix4 ] tri m. m. ;
+ [ { 1.0 0.0 0.0 } swap pitch>> <rotation-matrix4> ]
+ [ { 0.0 1.0 0.0 } swap yaw>> <rotation-matrix4> ]
+ [ location>> vneg <translation-matrix4> ] tri m. m. ;
: wasd-mv-inv-matrix ( world -- matrix )
- [ location>> translation-matrix4 ]
- [ { 0.0 -1.0 0.0 } swap yaw>> rotation-matrix4 ]
- [ { -1.0 0.0 0.0 } swap pitch>> rotation-matrix4 ] tri m. m. ;
+ [ location>> <translation-matrix4> ]
+ [ { 0.0 -1.0 0.0 } swap yaw>> <rotation-matrix4> ]
+ [ { -1.0 0.0 0.0 } swap pitch>> <rotation-matrix4> ] tri m. m. ;
: wasd-p-matrix ( world -- matrix )
p-matrix>> ;
world wasd-far-plane :> far-plane
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 )
loc world dim>> [ /f 0.5 - 2.0 * ] 2map
--- /dev/null
+Slava Pestov
--- /dev/null
+USING: help.markup help.syntax math sequences ;
+
+IN: math.matrices.elimination
+
+HELP: inverse
+{ $values { "matrix" sequence } }
+{ $description "Computes the multiplicative inverse of a matrix. Assuming the matrix is invertible." }
+{ $examples
+ "A matrix multiplied by its inverse is the identity matrix."
+ { $example
+ "USING: kernel math.matrices prettyprint ;"
+ "FROM: math.matrices.elimination => inverse ;"
+ "{ { 3 4 } { 7 9 } } dup inverse m. 2 <identity-matrix> = ."
+ "t"
+ }
+} ;
+
+HELP: echelon
+{ $values { "matrix" sequence } { "matrix'" sequence } }
+{ $description "Computes the reduced row-echelon form of the matrix." } ;
+
+HELP: nonzero-rows
+{ $values { "matrix" sequence } { "matrix'" sequence } }
+{ $description "Removes all all-zero rows from the matrix" }
+{ $examples
+ { $example
+ "USING: math.matrices.elimination prettyprint ;"
+ "{ { 0 0 } { 5 6 } { 0 0 } { 4 0 } } nonzero-rows ."
+ "{ { 5 6 } { 4 0 } }"
+ }
+} ;
+
+HELP: leading
+{ $values
+ { "seq" sequence }
+ { "n" "the index of the first match, or " { $snippet f } "." }
+ { "elt" "the first non-zero element, or " { $snippet f } "." }
+}
+{ $description "Find the first non-zero element of a sequence." } ;
--- /dev/null
+IN: math.matrices.elimination.tests
+USING: kernel math.matrices math.matrices.elimination
+tools.test sequences ;
+
+{
+ {
+ { 1 0 0 0 }
+ { 0 1 0 0 }
+ { 0 0 1 0 }
+ { 0 0 0 1 }
+ }
+} [
+ {
+ { 1 0 0 0 }
+ { 0 1 0 0 }
+ { 0 0 1 0 }
+ { 0 0 0 1 }
+ } echelon
+] unit-test
+
+{
+ {
+ { 1 0 0 0 }
+ { 0 1 0 0 }
+ { 0 0 1 0 }
+ { 0 0 0 1 }
+ }
+} [
+ {
+ { 1 0 0 0 }
+ { 1 1 0 0 }
+ { 1 0 1 0 }
+ { 1 0 0 1 }
+ } echelon
+] unit-test
+
+{
+ {
+ { 1 0 0 0 }
+ { 0 1 0 0 }
+ { 0 0 1 0 }
+ { 0 0 0 1 }
+ }
+} [
+ {
+ { 1 0 0 0 }
+ { 1 1 0 0 }
+ { 1 0 1 0 }
+ { 1 1 0 1 }
+ } echelon
+] unit-test
+
+{
+ {
+ { 1 0 0 0 }
+ { 0 1 0 0 }
+ { 0 0 1 0 }
+ { 0 0 0 1 }
+ }
+} [
+ {
+ { 1 0 0 0 }
+ { 1 1 0 0 }
+ { 1 1 0 1 }
+ { 1 0 1 0 }
+ } echelon
+] unit-test
+
+{
+ {
+ { 1 0 0 0 }
+ { 0 1 0 0 }
+ { 0 0 0 0 }
+ { 0 0 0 0 }
+ }
+} [
+ {
+ { 0 1 0 0 }
+ { 1 0 0 0 }
+ { 1 0 0 0 }
+ { 1 0 0 0 }
+ } [
+ [ 1 ] [ 0 0 pivot-row ] unit-test
+ 1 0 do-row
+ ] with-matrix
+] unit-test
+
+{
+ {
+ { 1 0 0 0 }
+ { 0 1 0 0 }
+ { 0 0 0 0 }
+ { 0 0 0 0 }
+ }
+} [
+ {
+ { 0 1 0 0 }
+ { 1 0 0 0 }
+ { 1 0 0 0 }
+ { 1 0 0 0 }
+ } echelon
+] unit-test
+
+{
+ {
+ { 1 0 0 0 }
+ { 0 1 0 0 }
+ { 0 0 0 1 }
+ { 0 0 0 0 }
+ }
+} [
+ {
+ { 1 0 0 0 }
+ { 0 1 0 0 }
+ { 1 0 0 1 }
+ { 1 0 0 1 }
+ } echelon
+] unit-test
+
+{
+ {
+ { 1 0 0 1 }
+ { 0 1 0 1 }
+ { 0 0 0 -1 }
+ { 0 0 0 0 }
+ }
+} [
+ {
+ { 0 1 0 1 }
+ { 1 0 0 1 }
+ { 1 0 0 0 }
+ { 1 1 0 1 }
+ } echelon
+] unit-test
+
+{
+ 2
+} [
+ {
+ { 0 0 }
+ { 0 0 }
+ } nullspace length
+] unit-test
+
+{
+ 1 3
+} [
+ {
+ { 0 1 0 1 }
+ { 1 0 0 1 }
+ { 1 0 0 0 }
+ { 1 1 0 1 }
+ } null/rank
+] unit-test
+
+{
+ 1 3
+} [
+ {
+ { 0 0 0 0 0 1 0 1 }
+ { 0 0 0 0 1 0 0 1 }
+ { 0 0 0 0 1 0 0 0 }
+ { 0 0 0 0 1 1 0 1 }
+ } null/rank
+] unit-test
+
+{ { { 1 0 -1 } { 0 1 2 } } }
+[ { { 1 2 3 } { 4 5 6 } } solution ] unit-test
--- /dev/null
+! Copyright (C) 2006, 2010 Slava Pestov.
+! See http://factorcode.org/license.txt for BSD license.
+USING: kernel locals math math.vectors math.matrices
+namespaces sequences fry sorting ;
+IN: math.matrices.elimination
+
+SYMBOL: matrix
+
+: with-matrix ( matrix quot -- )
+ matrix swap [ matrix get ] compose with-variable ; inline
+
+: nth-row ( row# -- seq ) matrix get nth ;
+
+: change-row ( ..a row# quot: ( ..a seq -- ..b seq ) -- ..b )
+ matrix get swap change-nth ; inline
+
+: exchange-rows ( row# row# -- ) matrix get exchange ;
+
+: rows ( -- n ) matrix get length ;
+
+: cols ( -- n ) 0 nth-row length ;
+
+: skip ( i seq quot -- n )
+ over [ find-from drop ] dip swap [ ] [ length ] ?if ; inline
+
+: first-col ( row# -- n )
+ ! First non-zero column
+ 0 swap nth-row [ zero? not ] skip ;
+
+: clear-scale ( col# pivot-row i-row -- n )
+ overd nth dup zero? [
+ 3drop 0
+ ] [
+ [ nth dup zero? ] dip swap [
+ 2drop 0
+ ] [
+ swap / neg
+ ] if
+ ] if ;
+
+: (clear-col) ( col# pivot-row i -- )
+ [ [ clear-scale ] 2keep [ n*v ] dip v+ ] change-row ;
+
+: rows-from ( row# -- slice )
+ rows dup <iota> <slice> ;
+
+: clear-col ( col# row# rows -- )
+ [ nth-row ] dip [ [ 2dup ] dip (clear-col) ] each 2drop ;
+
+: do-row ( exchange-with row# -- )
+ [ exchange-rows ] keep
+ [ first-col ] keep
+ dup 1 + rows-from clear-col ;
+
+: find-row ( row# quot -- i elt )
+ [ rows-from ] dip find ; inline
+
+: pivot-row ( col# row# -- n )
+ [ dupd nth-row nth zero? not ] find-row 2nip ;
+
+: (echelon) ( col# row# -- )
+ over cols < over rows < and [
+ 2dup pivot-row [ over do-row 1 + ] when*
+ [ 1 + ] dip (echelon)
+ ] [
+ 2drop
+ ] if ;
+
+: echelon ( matrix -- matrix' )
+ [ 0 0 (echelon) ] with-matrix ;
+
+: nonzero-rows ( matrix -- matrix' )
+ [ [ zero? ] all? ] reject ;
+
+: null/rank ( matrix -- null rank )
+ echelon dup length swap nonzero-rows length [ - ] keep ;
+
+: leading ( seq -- n elt ) [ zero? not ] find ;
+
+: reduced ( matrix' -- matrix'' )
+ [
+ rows <iota> <reversed> [
+ dup nth-row leading drop
+ [ swap dup <iota> clear-col ] [ drop ] if*
+ ] each
+ ] with-matrix ;
+
+:: basis-vector ( row col# -- )
+ row clone :> row'
+ col# row' nth neg recip :> a
+ 0 col# row' set-nth
+ a row n*v col# matrix get set-nth ;
+
+: nullspace ( matrix -- seq )
+ echelon reduced dup empty? [
+ dup first length <identity-matrix> [
+ [
+ dup leading drop
+ [ basis-vector ] [ drop ] if*
+ ] each
+ ] with-matrix flip nonzero-rows
+ ] unless ;
+
+: 1-pivots ( matrix -- matrix )
+ [ dup leading nip [ recip v*n ] when* ] map ;
+
+: solution ( matrix -- matrix )
+ echelon nonzero-rows reduced 1-pivots ;
+
+: inverse ( matrix -- matrix ) ! Assumes an invertible matrix
+ dup length
+ [ <identity-matrix> [ append ] 2map solution ] keep
+ [ tail ] curry map ;
--- /dev/null
+Solving systems of linear equations
--- /dev/null
+Slava Pestov
+Joe Groff
+Doug Coleman
+Cat Stevens
--- /dev/null
+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 } }
+;
--- /dev/null
+! 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
--- /dev/null
+! 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
--- /dev/null
+Matrix arithmetic - extra and miscellaneous words
\ No newline at end of file
! The storage functions
: <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 )
swap '[ drop _ ] mmap! ;
: set-pixel ( {R,G,B} {i,j} image -- ) set-Mi,j ; inline
dup conj-t [ m. ] [ swap m. ] 2bi = ;
: unitary-matrix? ( matrix -- ? )
- [ dup conj-t m. ] [ length identity-matrix ] bi = ;
+ [ dup conj-t m. ] [ length <identity-matrix> ] bi = ;
projects / factor.git / commitdiff