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