libs/topology/matrix.factor

   1 ! Copyright (C) 2006 Slava Pestov.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 IN: matrices
   4 USING: kernel math namespaces parser sequences ;
   5
   6 SYMBOL: matrix
   7
   8 : with-matrix ( matrix quot -- )
   9     [ swap matrix set call matrix get ] with-scope ; inline
  10
  11 : nth-row ( row# -- seq ) matrix get nth ;
  12
  13 : nth-col ( col# ignore-rows -- seq )
  14     matrix get tail-slice [ nth ] map-with ;
  15
  16 : change-row ( row# quot -- | quot: seq -- seq )
  17     matrix get swap change-nth ; inline
  18
  19 : exchange-rows ( row# row# -- ) matrix get exchange ;
  20
  21 : rows ( -- n ) matrix get length ;
  22
  23 : cols ( -- n ) 0 nth-row length ;
  24
  25 : first-col ( row# -- n )
  26     #! First non-zero column
  27     0 swap nth-row [ zero? not ] skip ;
  28
  29 : clear-scale ( col# pivot-row i-row -- n )
  30     >r over r> nth dup zero? [
  31         3drop 0
  32     ] [
  33         >r nth dup zero? [
  34             r> 2drop 0
  35         ] [
  36             r> swap / neg
  37         ] if
  38     ] if ;
  39
  40 : (clear-col) ( col# pivot-row i -- )
  41     [ [ clear-scale ] 2keep >r n*v r> v+ ] change-row ;
  42
  43 : (each-row) ( row# -- slice )
  44     rows dup <slice> ;
  45
  46 : each-row ( row# quot -- )
  47     >r (each-row) r> each ; inline
  48
  49 : clear-col ( col# row# -- )
  50     [ nth-row ] keep 1+
  51     [ >r 2dup r> (clear-col) ] each-row
  52     2drop ;
  53
  54 : do-row ( exchange-with row# -- )
  55     [ exchange-rows ] keep
  56     [ first-col ] keep
  57     clear-col ;
  58
  59 : find-row ( row# quot -- i elt )
  60     >r (each-row) r> find ; inline
  61
  62 : pivot-row ( col# row# -- n )
  63     [ dupd nth-row nth zero? not ] find-row 2nip ;
  64
  65 : (row-reduce) ( col# row# -- )
  66     over cols < over rows < and [
  67         2dup pivot-row [ over do-row 1+ ] when* >r 1+ r>
  68         (row-reduce)
  69     ] [
  70         2drop
  71     ] if ;
  72
  73 : row-reduce ( matrix -- matrix' )
  74     [ 0 0 (row-reduce) ] with-matrix ;
  75
  76 : null/rank ( matrix -- null rank )
  77     row-reduce [ [ [ zero? ] all? ] subset ] keep
  78     [ length ] 2apply over - ;