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 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? ] assoc-reject ;
  23
  24 SYMBOL: terms
  25
  26 : with-terms ( quot -- hash )
  27     [
  28         H{ } clone terms namespaces: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 ] with-variables ;
  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 count ;
  75
  76 : inversions ( seq -- n )
  77     0 swap [ length <iota> ] keep [
  78         [ nth ] 2keep swap 1 + tail-slice (inversions) +
  79     ] curry each ;
  80
  81 : (wedge) ( n basis1 basis2 -- n basis )
  82     append dup all-unique? not [
  83         2drop 0 { }
  84     ] [
  85         dup permutation inversions -1^ rot *
  86         swap natural-sort
  87     ] if ;
  88
  89 : wedge ( x y -- x.y )
  90     [ >alt ] bi@ [
  91         swap building get '[
  92             [
  93                 2swap [
  94                     swapd * -rot (wedge) _ at+
  95                 ] 2keep
  96             ] assoc-each 2drop
  97         ] curry assoc-each
  98     ] H{ } make canonicalize ;
  99
 100 ! Differential
 101 SYMBOL: boundaries
 102
 103 : d= ( value basis -- )
 104     boundaries [ ?set-at ] change ;
 105
 106 : get-boundary ( basis -- value ) boundaries get at ;
 107
 108 : dx.y ( x y -- vec ) [ get-boundary ] dip wedge ;
 109
 110 DEFER: (d)
 111
 112 : x.dy ( x y -- vec ) (d) wedge -1 alt*n ;
 113
 114 : (d) ( product -- value )
 115     [ H{ } ] [ unclip swap [ x.dy ] 2keep dx.y alt+ ] if-empty ;
 116
 117 : linear-op ( vec quot -- vec )
 118         [
 119         [
 120             -rot [ swap call ] dip alt*n (alt+)
 121         ] curry assoc-each
 122     ] with-terms ; inline
 123
 124 : d ( x -- dx )
 125     >alt [ (d) ] linear-op ;
 126
 127 ! Interior product
 128 : (interior) ( y basis-elt -- i_y[basis-elt] )
 129     2dup index dup [
 130         -rot remove associate
 131     ] [
 132         3drop 0
 133     ] if ;
 134
 135 : interior ( x y -- i_y[x] )
 136     ! y is a generator
 137     swap >alt [ dupd (interior) ] linear-op nip ;
 138
 139 ! Computing a basis
 140 : graded ( seq -- seq )
 141     dup 0 [ length max ] reduce 1 + [ V{ } clone ] replicate
 142     [ dup length pick nth push ] reduce ;
 143
 144 : nth-basis-elt ( generators n -- elt )
 145     over length <iota> [
 146         3dup bit? [ nth ] [ 2drop f ] if
 147     ] map sift 2nip ;
 148
 149 : basis ( generators -- seq )
 150     natural-sort dup length 2^ <iota> [ nth-basis-elt ] with map ;
 151
 152 : (tensor) ( seq1 seq2 -- seq )
 153     [
 154         [ prepend natural-sort ] curry map
 155     ] with map concat ;
 156
 157 : tensor ( graded-basis1 graded-basis2 -- bigraded-basis )
 158     [ [ swap (tensor) ] curry map ] with map ;
 159
 160 ! Computing cohomology
 161 : (op-matrix) ( range quot basis-elt -- row )
 162     swap call [ at 0 or ] curry map ; inline
 163
 164 : op-matrix ( domain range quot -- matrix )
 165     rot [ (op-matrix) ] 2with map ; inline
 166
 167 : d-matrix ( domain range -- matrix )
 168     [ (d) ] op-matrix ;
 169
 170 : dim-im/ker-d ( domain range -- null/rank )
 171     d-matrix null/rank 2array ;
 172
 173 ! Graded by degree
 174 : (graded-ker/im-d) ( n seq -- null/rank )
 175     ! d: C(n) ---> C(n+1)
 176     [ ?nth ] [ [ 1 + ] dip ?nth ] 2bi
 177     dim-im/ker-d ;
 178
 179 : graded-ker/im-d ( graded-basis -- seq )
 180     [ length <iota> ] keep [ (graded-ker/im-d) ] curry map ;
 181
 182 : graded-betti ( generators -- seq )
 183     basis graded graded-ker/im-d unzip but-last 0 prefix v- ;
 184
 185 ! Bi-graded for two-step complexes
 186 : (bigraded-ker/im-d) ( u-deg z-deg bigraded-basis -- null/rank )
 187     ! d: C(u,z) ---> C(u+2,z-1)
 188     [ ?nth ?nth ] 3keep [ [ 2 + ] dip 1 - ] dip ?nth ?nth
 189     dim-im/ker-d ;
 190
 191 :: bigraded-ker/im-d ( basis -- seq )
 192     basis length <iota> [| z |
 193          basis first length <iota> [| u |
 194             u z basis (bigraded-ker/im-d)
 195         ] map
 196     ] map ;
 197
 198 : bigraded-betti ( u-generators z-generators -- seq )
 199     [ basis graded ] bi@ tensor bigraded-ker/im-d
 200     [ [ keys ] map ] keep
 201     [ values 2 head* { 0 0 } prepend ] map
 202     rest dup first length 0 <array> suffix
 203     [ v- ] 2map ;
 204
 205 ! Laplacian
 206 : m.m' ( matrix -- matrix' ) dup flip m. ;
 207 : m'.m ( matrix -- matrix' ) dup flip swap m. ;
 208
 209 : empty-matrix? ( matrix -- ? )
 210     [ t ] [ first empty? ] if-empty ;
 211
 212 : ?m+ ( m1 m2 -- m3 )
 213     over empty-matrix? [
 214         nip
 215     ] [
 216         dup empty-matrix? [
 217             drop
 218         ] [
 219             m+
 220         ] if
 221     ] if ;
 222
 223 : laplacian-matrix ( basis1 basis2 basis3 -- matrix )
 224     dupd d-matrix m.m' [ d-matrix m'.m ] dip ?m+ ;
 225
 226 : laplacian-betti ( basis1 basis2 basis3 -- n )
 227     laplacian-matrix null/rank drop ;
 228
 229 :: laplacian-kernel ( basis1 basis2 basis3 -- basis )
 230     basis1 basis2 basis3 laplacian-matrix :> lap
 231     lap empty-matrix? [ f ] [
 232         lap nullspace [| x |
 233             basis2 x [ [ wedge (alt+) ] 2each ] with-terms
 234         ] map
 235     ] if ;
 236
 237 : graded-triple ( seq n -- triple )
 238     3 [ 1 - + ] with map swap [ ?nth ] curry map ;
 239
 240 : graded-triples ( seq -- triples )
 241     dup length [ graded-triple ] with map ;
 242
 243 : graded-laplacian ( generators quot -- seq )
 244     [ basis graded graded-triples [ first3 ] ] dip compose map ; inline
 245
 246 : graded-laplacian-betti ( generators -- seq )
 247     [ laplacian-betti ] graded-laplacian ;
 248
 249 : graded-laplacian-kernel ( generators -- seq )
 250     [ laplacian-kernel ] graded-laplacian ;
 251
 252 : graded-basis. ( seq -- )
 253     [
 254         "=== Degree " write pprint
 255         ": dimension " write dup length .
 256         [ alt. ] each
 257     ] each-index ;
 258
 259 : bigraded-triple ( u-deg z-deg bigraded-basis -- triple )
 260     ! d: C(u,z) ---> C(u+2,z-1)
 261     [ [ 2 - ] [ 1 + ] [ ] tri* ?nth ?nth ]
 262     [ ?nth ?nth ]
 263     [ [ 2 + ] [ 1 - ] [ ] tri* ?nth ?nth ]
 264     3tri
 265     3array ;
 266
 267 :: bigraded-triples ( grid -- triples )
 268     grid length <iota> [| z |
 269         grid first length <iota> [| u |
 270             u z grid bigraded-triple
 271         ] map
 272     ] map ;
 273
 274 : bigraded-laplacian ( u-generators z-generators quot -- seq )
 275     [ [ basis graded ] bi@ tensor bigraded-triples ] dip
 276     [ [ first3 ] prepose map ] curry map ; inline
 277
 278 : bigraded-laplacian-betti ( u-generators z-generators -- seq )
 279     [ laplacian-betti ] bigraded-laplacian ;
 280
 281 : bigraded-laplacian-kernel ( u-generators z-generators -- seq )
 282     [ laplacian-kernel ] bigraded-laplacian ;
 283
 284 : bigraded-basis. ( seq -- )
 285     [
 286         "=== U-degree " write .
 287         [
 288             "  === Z-degree " write pprint
 289             ": dimension " write dup length .
 290             [ "  " write alt. ] each
 291         ] each-index
 292     ] each-index ;