extra/koszul/koszul.factor

   1 ! Copyright (C) 2006, 2007 Slava Pestov.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: accessors arrays assocs combinators fry hashtables io
   4 kernel locals make math math.matrices math.matrices.elimination
   5 math.order math.parser math.vectors namespaces prettyprint
   6 sequences sets shuffle sorting splitting ;
   7 FROM: namespaces => set ;
   8 IN: koszul
   9
  10 ! Utilities
  11 : -1^ ( m -- n ) odd? -1 1 ? ;
  12
  13 : >alt ( obj -- vec )
  14     {
  15         { [ dup not ] [ drop 0 >alt ] }
  16         { [ dup number? ] [ { } associate ] }
  17         { [ dup array? ] [ 1 swap associate ] }
  18         { [ dup hashtable? ] [ ] }
  19         [ 1array >alt ]
  20     } cond ;
  21
  22 : canonicalize ( assoc -- assoc' )
  23     [ nip zero? ] assoc-reject ;
  24
  25 SYMBOL: terms
  26
  27 : with-terms ( quot -- hash )
  28     [
  29         H{ } clone terms set call terms get canonicalize
  30     ] with-scope ; inline
  31
  32 ! Printing elements
  33 : num-alt. ( n -- str )
  34     {
  35         { 1 [ " + " ] }
  36         { -1 [ " - " ] }
  37         [ number>string " + " prepend ]
  38     } case ;
  39
  40 : (alt.) ( basis n -- str )
  41     over empty? [
  42         nip number>string
  43     ] [
  44         num-alt.
  45         swap [ name>> ] map "." join
  46         append
  47     ] if ;
  48
  49 : alt. ( assoc -- )
  50     dup assoc-empty? [
  51         drop 0 .
  52     ] [
  53         [ (alt.) ] { } assoc>map concat " + " ?head drop print
  54     ] if ;
  55
  56 ! Addition
  57 : (alt+) ( x -- )
  58     terms get [ [ swap +@ ] assoc-each ] with-variables ;
  59
  60 : alt+ ( x y -- x+y )
  61     [ >alt ] bi@ [ (alt+) (alt+) ] with-terms ;
  62
  63 ! Multiplication
  64 : alt*n ( vec n -- vec )
  65     dup zero? [
  66         2drop H{ }
  67     ] [
  68         [ * ] curry assoc-map
  69     ] if ;
  70
  71 : permutation ( seq -- perm )
  72     [ natural-sort ] keep [ index ] curry map ;
  73
  74 : (inversions) ( n seq -- n )
  75     [ > ] with count ;
  76
  77 : inversions ( seq -- n )
  78     0 swap [ length iota ] keep [
  79         [ nth ] 2keep swap 1 + tail-slice (inversions) +
  80     ] curry each ;
  81
  82 : (wedge) ( n basis1 basis2 -- n basis )
  83     append dup all-unique? not [
  84         2drop 0 { }
  85     ] [
  86         dup permutation inversions -1^ rot *
  87         swap natural-sort
  88     ] if ;
  89
  90 : wedge ( x y -- x.y )
  91     [ >alt ] bi@ [
  92         swap building get '[
  93             [
  94                 2swap [
  95                     swapd * -rot (wedge) _ at+
  96                 ] 2keep
  97             ] assoc-each 2drop
  98         ] curry assoc-each
  99     ] H{ } make canonicalize ;
 100
 101 ! Differential
 102 SYMBOL: boundaries
 103
 104 : d= ( value basis -- )
 105     boundaries [ ?set-at ] change ;
 106
 107 : ((d)) ( basis -- value ) boundaries get at ;
 108
 109 : dx.y ( x y -- vec ) [ ((d)) ] dip wedge ;
 110
 111 DEFER: (d)
 112
 113 : x.dy ( x y -- vec ) (d) wedge -1 alt*n ;
 114
 115 : (d) ( product -- value )
 116     [ H{ } ] [ unclip swap [ x.dy ] 2keep dx.y alt+ ] if-empty ;
 117
 118 : linear-op ( vec quot -- vec )
 119         [
 120         [
 121             -rot [ swap call ] dip alt*n (alt+)
 122         ] curry assoc-each
 123     ] with-terms ; inline
 124
 125 : d ( x -- dx )
 126     >alt [ (d) ] linear-op ;
 127
 128 ! Interior product
 129 : (interior) ( y basis-elt -- i_y[basis-elt] )
 130     2dup index dup [
 131         -rot remove associate
 132     ] [
 133         3drop 0
 134     ] if ;
 135
 136 : interior ( x y -- i_y[x] )
 137     #! y is a generator
 138     swap >alt [ dupd (interior) ] linear-op nip ;
 139
 140 ! Computing a basis
 141 : graded ( seq -- seq )
 142     dup 0 [ length max ] reduce 1 + [ V{ } clone ] replicate
 143     [ dup length pick nth push ] reduce ;
 144
 145 : nth-basis-elt ( generators n -- elt )
 146     over length iota [
 147         3dup bit? [ nth ] [ 2drop f ] if
 148     ] map sift 2nip ;
 149
 150 : basis ( generators -- seq )
 151     natural-sort dup length 2^ iota [ nth-basis-elt ] with map ;
 152
 153 : (tensor) ( seq1 seq2 -- seq )
 154     [
 155         [ prepend natural-sort ] curry map
 156     ] with map concat ;
 157
 158 : tensor ( graded-basis1 graded-basis2 -- bigraded-basis )
 159     [ [ swap (tensor) ] curry map ] with map ;
 160
 161 ! Computing cohomology
 162 : (op-matrix) ( range quot basis-elt -- row )
 163     swap call [ at 0 or ] curry map ; inline
 164
 165 : op-matrix ( domain range quot -- matrix )
 166     rot [ (op-matrix) ] 2with map ; inline
 167
 168 : d-matrix ( domain range -- matrix )
 169     [ (d) ] op-matrix ;
 170
 171 : dim-im/ker-d ( domain range -- null/rank )
 172     d-matrix null/rank 2array ;
 173
 174 ! Graded by degree
 175 : (graded-ker/im-d) ( n seq -- null/rank )
 176     #! d: C(n) ---> C(n+1)
 177     [ ?nth ] [ [ 1 + ] dip ?nth ] 2bi
 178     dim-im/ker-d ;
 179
 180 : graded-ker/im-d ( graded-basis -- seq )
 181     [ length iota ] keep [ (graded-ker/im-d) ] curry map ;
 182
 183 : graded-betti ( generators -- seq )
 184     basis graded graded-ker/im-d unzip but-last 0 prefix v- ;
 185
 186 ! Bi-graded for two-step complexes
 187 : (bigraded-ker/im-d) ( u-deg z-deg bigraded-basis -- null/rank )
 188     #! d: C(u,z) ---> C(u+2,z-1)
 189     [ ?nth ?nth ] 3keep [ [ 2 + ] dip 1 - ] dip ?nth ?nth
 190     dim-im/ker-d ;
 191
 192 :: bigraded-ker/im-d ( basis -- seq )
 193     basis length iota [| z |
 194          basis first length iota [| u |
 195             u z basis (bigraded-ker/im-d)
 196         ] map
 197     ] map ;
 198
 199 : bigraded-betti ( u-generators z-generators -- seq )
 200     [ basis graded ] bi@ tensor bigraded-ker/im-d
 201     [ [ keys ] map ] keep
 202     [ values 2 head* { 0 0 } prepend ] map
 203     rest dup first length 0 <array> suffix
 204     [ v- ] 2map ;
 205
 206 ! Laplacian
 207 : m.m' ( matrix -- matrix' ) dup flip m. ;
 208 : m'.m ( matrix -- matrix' ) dup flip swap m. ;
 209
 210 : empty-matrix? ( matrix -- ? )
 211     [ t ] [ first empty? ] if-empty ;
 212
 213 : ?m+ ( m1 m2 -- m3 )
 214     over empty-matrix? [
 215         nip
 216     ] [
 217         dup empty-matrix? [
 218             drop
 219         ] [
 220             m+
 221         ] if
 222     ] if ;
 223
 224 : laplacian-matrix ( basis1 basis2 basis3 -- matrix )
 225     dupd d-matrix m.m' [ d-matrix m'.m ] dip ?m+ ;
 226
 227 : laplacian-betti ( basis1 basis2 basis3 -- n )
 228     laplacian-matrix null/rank drop ;
 229
 230 :: laplacian-kernel ( basis1 basis2 basis3 -- basis )
 231     basis1 basis2 basis3 laplacian-matrix :> lap
 232     lap empty-matrix? [ f ] [
 233         lap nullspace [| x |
 234             basis2 x [ [ wedge (alt+) ] 2each ] with-terms
 235         ] map
 236     ] if ;
 237
 238 : graded-triple ( seq n -- triple )
 239     3 [ 1 - + ] with map swap [ ?nth ] curry map ;
 240
 241 : graded-triples ( seq -- triples )
 242     dup length [ graded-triple ] with map ;
 243
 244 : graded-laplacian ( generators quot -- seq )
 245     [ basis graded graded-triples [ first3 ] ] dip compose map ; inline
 246
 247 : graded-laplacian-betti ( generators -- seq )
 248     [ laplacian-betti ] graded-laplacian ;
 249
 250 : graded-laplacian-kernel ( generators -- seq )
 251     [ laplacian-kernel ] graded-laplacian ;
 252
 253 : graded-basis. ( seq -- )
 254     [
 255         "=== Degree " write pprint
 256         ": dimension " write dup length .
 257         [ alt. ] each
 258     ] each-index ;
 259
 260 : bigraded-triple ( u-deg z-deg bigraded-basis -- triple )
 261     #! d: C(u,z) ---> C(u+2,z-1)
 262     [ [ 2 - ] [ 1 + ] [ ] tri* ?nth ?nth ]
 263     [ ?nth ?nth ]
 264     [ [ 2 + ] [ 1 - ] [ ] tri* ?nth ?nth ]
 265     3tri
 266     3array ;
 267
 268 :: bigraded-triples ( grid -- triples )
 269     grid length iota [| z |
 270         grid first length iota [| u |
 271             u z grid bigraded-triple
 272         ] map
 273     ] map ;
 274
 275 : bigraded-laplacian ( u-generators z-generators quot -- seq )
 276     [ [ basis graded ] bi@ tensor bigraded-triples ] dip
 277     [ [ first3 ] prepose map ] curry map ; inline
 278
 279 : bigraded-laplacian-betti ( u-generators z-generators -- seq )
 280     [ laplacian-betti ] bigraded-laplacian ;
 281
 282 : bigraded-laplacian-kernel ( u-generators z-generators -- seq )
 283     [ laplacian-kernel ] bigraded-laplacian ;
 284
 285 : bigraded-basis. ( seq -- )
 286     [
 287         "=== U-degree " write .
 288         [
 289             "  === Z-degree " write pprint
 290             ": dimension " write dup length .
 291             [ "  " write alt. ] each
 292         ] each-index
 293     ] each-index ;