extra/koszul/koszul.factor

   1 ! Copyright (C) 2006, 2007 Slava Pestov.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: arrays assocs 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 symbols 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 [ word-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 ) >r ((d)) r> wedge ;
 112
 113 DEFER: (d)
 114
 115 : x.dy ( x y -- vec ) (d) wedge -1 alt*n ;
 116
 117 : (d) ( product -- value )
 118     dup empty?
 119     [ drop H{ } ] [ unclip swap [ x.dy ] 2keep dx.y alt+ ] if ;
 120
 121 : linear-op ( vec quot -- vec )
 122         [
 123         [
 124             -rot >r swap call r> alt*n (alt+)
 125         ] curry assoc-each
 126     ] with-terms ; inline
 127
 128 : d ( x -- dx )
 129     >alt [ (d) ] linear-op ;
 130
 131 ! Interior product
 132 : (interior) ( y basis-elt -- i_y[basis-elt] )
 133     2dup index dup [
 134         -rot remove associate
 135     ] [
 136         3drop 0
 137     ] if ;
 138
 139 : interior ( x y -- i_y[x] )
 140     #! y is a generator
 141     swap >alt [ dupd (interior) ] linear-op nip ;
 142
 143 ! Computing a basis
 144 : graded ( seq -- seq )
 145     dup 0 [ length max ] reduce 1+ [ drop V{ } clone ] map
 146     [ dup length pick nth push ] reduce ;
 147
 148 : nth-basis-elt ( generators n -- elt )
 149     over length [
 150         3dup bit? [ nth ] [ 2drop f ] if
 151     ] map sift 2nip ;
 152
 153 : basis ( generators -- seq )
 154     natural-sort dup length 2^ [ nth-basis-elt ] with map ;
 155
 156 : (tensor) ( seq1 seq2 -- seq )
 157     [
 158         [ prepend natural-sort ] curry map
 159     ] with map concat ;
 160
 161 : tensor ( graded-basis1 graded-basis2 -- bigraded-basis )
 162     [ [ swap (tensor) ] curry map ] with map ;
 163
 164 ! Computing cohomology
 165 : (op-matrix) ( range quot basis-elt -- row )
 166     swap call [ at 0 or ] curry map ; inline
 167
 168 : op-matrix ( domain range quot -- matrix )
 169     rot [ >r 2dup r> (op-matrix) ] map 2nip ; inline
 170
 171 : d-matrix ( domain range -- matrix )
 172     [ (d) ] op-matrix ;
 173
 174 : dim-im/ker-d ( domain range -- null/rank )
 175     d-matrix null/rank 2array ;
 176
 177 ! Graded by degree
 178 : (graded-ker/im-d) ( n seq -- null/rank )
 179     #! d: C(n) ---> C(n+1)
 180     [ ?nth ] 2keep >r 1+ r> ?nth
 181     dim-im/ker-d ;
 182
 183 : graded-ker/im-d ( graded-basis -- seq )
 184     [ length ] keep [ (graded-ker/im-d) ] curry map ;
 185
 186 : graded-betti ( generators -- seq )
 187     basis graded graded-ker/im-d flip first2 but-last 0 prefix v- ;
 188
 189 ! Bi-graded for two-step complexes
 190 : (bigraded-ker/im-d) ( u-deg z-deg bigraded-basis -- null/rank )
 191     #! d: C(u,z) ---> C(u+2,z-1)
 192     [ ?nth ?nth ] 3keep >r >r 2 + r> 1 - r> ?nth ?nth
 193     dim-im/ker-d ;
 194
 195 : bigraded-ker/im-d ( bigraded-basis -- seq )
 196     dup length [
 197         over first length [
 198             >r 2dup r> spin (bigraded-ker/im-d)
 199         ] map 2nip
 200     ] with map ;
 201
 202 : bigraded-betti ( u-generators z-generators -- seq )
 203     [ basis graded ] bi@ tensor bigraded-ker/im-d
 204     [ [ [ first ] map ] map ] keep
 205     [ [ second ] map 2 head* { 0 0 } prepend ] map
 206     rest dup first length 0 <array> suffix
 207     [ v- ] 2map ;
 208
 209 ! Laplacian
 210 : m.m' ( matrix -- matrix' ) dup flip m. ;
 211 : m'.m ( matrix -- matrix' ) dup flip swap m. ;
 212
 213 : empty-matrix? ( matrix -- ? )
 214     dup empty? [ drop t ] [ first empty? ] if ;
 215
 216 : ?m+ ( m1 m2 -- m3 )
 217     over empty-matrix? [
 218         nip
 219     ] [
 220         dup empty-matrix? [
 221             drop
 222         ] [
 223             m+
 224         ] if
 225     ] if ;
 226
 227 : laplacian-matrix ( basis1 basis2 basis3 -- matrix )
 228     dupd d-matrix m.m' >r d-matrix m'.m r> ?m+ ;
 229
 230 : laplacian-betti ( basis1 basis2 basis3 -- n )
 231     laplacian-matrix null/rank drop ;
 232
 233 : laplacian-kernel ( basis1 basis2 basis3 -- basis )
 234     >r tuck r>
 235     laplacian-matrix dup empty-matrix? [
 236         2drop f
 237     ] [
 238         nullspace [
 239             [ [ wedge (alt+) ] 2each ] with-terms
 240         ] with map
 241     ] if ;
 242
 243 : graded-triple ( seq n -- triple )
 244     3 [ 1- + ] with map swap [ ?nth ] curry map ;
 245
 246 : graded-triples ( seq -- triples )
 247     dup length [ graded-triple ] with map ;
 248
 249 : graded-laplacian ( generators quot -- seq )
 250     >r basis graded graded-triples [ first3 ] r> compose map ;
 251     inline
 252
 253 : graded-laplacian-betti ( generators -- seq )
 254     [ laplacian-betti ] graded-laplacian ;
 255
 256 : graded-laplacian-kernel ( generators -- seq )
 257     [ laplacian-kernel ] graded-laplacian ;
 258
 259 : graded-basis. ( seq -- )
 260     dup length [
 261         "=== Degree " write pprint
 262         ": dimension " write dup length .
 263         [ alt. ] each
 264     ] 2each ;
 265
 266 : bigraded-triple ( u-deg z-deg bigraded-basis -- triple )
 267     #! d: C(u,z) ---> C(u+2,z-1)
 268     [ >r >r 2 - r> 1 + r> ?nth ?nth ] 3keep
 269     [ ?nth ?nth ] 3keep
 270     >r >r 2 + r> 1 - r> ?nth ?nth
 271     3array ;
 272
 273 : bigraded-triples ( grid -- triples )
 274     dup length [
 275         over first length [
 276             >r 2dup r> spin bigraded-triple
 277         ] map 2nip
 278     ] with map ;
 279
 280 : bigraded-laplacian ( u-generators z-generators quot -- seq )
 281     >r [ basis graded ] bi@ tensor bigraded-triples r>
 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     dup length [
 292         "=== U-degree " write .
 293         dup length [
 294             "  === Z-degree " write pprint
 295             ": dimension " write dup length .
 296             [ "  " write alt. ] each
 297         ] 2each
 298     ] 2each ;