1 ! Copyright (C) 2006, 2007 Slava Pestov.
2 ! See http://factorcode.org/license.txt for BSD license.
3 USING: accessors 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 sets math.order ;
10 : -1^ ( m -- n ) odd? -1 1 ? ;
14 { [ dup not ] [ drop 0 >alt ] }
15 { [ dup number? ] [ { } associate ] }
16 { [ dup array? ] [ 1 swap associate ] }
17 { [ dup hashtable? ] [ ] }
21 : canonicalize ( assoc -- assoc' )
22 [ nip zero? not ] assoc-filter ;
26 : with-terms ( quot -- hash )
28 H{ } clone terms set call terms get canonicalize
32 : num-alt. ( n -- str )
36 [ number>string " + " prepend ]
39 : (alt.) ( basis n -- str )
44 swap [ name>> ] map "." join
52 [ (alt.) ] { } assoc>map concat " + " ?head drop print
57 terms get [ [ swap +@ ] assoc-each ] bind ;
60 [ >alt ] bi@ [ (alt+) (alt+) ] with-terms ;
63 : alt*n ( vec n -- vec )
70 : permutation ( seq -- perm )
71 [ natural-sort ] keep [ index ] curry map ;
73 : (inversions) ( n seq -- n )
74 [ > ] with filter length ;
76 : inversions ( seq -- n )
77 0 swap [ length ] keep [
78 [ nth ] 2keep swap 1+ tail-slice (inversions) +
81 : duplicates? ( seq -- ? )
82 dup prune [ length ] bi@ > ;
84 : (wedge) ( n basis1 basis2 -- n basis )
85 append dup duplicates? [
88 dup permutation inversions -1^ rot *
92 : wedge ( x y -- x.y )
97 swapd * -rot (wedge) +@
101 ] H{ } make-assoc canonicalize ;
106 : d= ( value basis -- )
107 boundaries [ ?set-at ] change ;
109 : ((d)) ( basis -- value ) boundaries get at ;
111 : dx.y ( x y -- vec ) [ ((d)) ] dip wedge ;
115 : x.dy ( x y -- vec ) (d) wedge -1 alt*n ;
117 : (d) ( product -- value )
118 [ H{ } ] [ unclip swap [ x.dy ] 2keep dx.y alt+ ] if-empty ;
120 : linear-op ( vec quot -- vec )
123 -rot [ swap call ] dip alt*n (alt+)
125 ] with-terms ; inline
128 >alt [ (d) ] linear-op ;
131 : (interior) ( y basis-elt -- i_y[basis-elt] )
133 -rot remove associate
138 : interior ( x y -- i_y[x] )
140 swap >alt [ dupd (interior) ] linear-op nip ;
143 : graded ( seq -- seq )
144 dup 0 [ length max ] reduce 1+ [ V{ } clone ] replicate
145 [ dup length pick nth push ] reduce ;
147 : nth-basis-elt ( generators n -- elt )
149 3dup bit? [ nth ] [ 2drop f ] if
152 : basis ( generators -- seq )
153 natural-sort dup length 2^ [ nth-basis-elt ] with map ;
155 : (tensor) ( seq1 seq2 -- seq )
157 [ prepend natural-sort ] curry map
160 : tensor ( graded-basis1 graded-basis2 -- bigraded-basis )
161 [ [ swap (tensor) ] curry map ] with map ;
163 ! Computing cohomology
164 : (op-matrix) ( range quot basis-elt -- row )
165 swap call [ at 0 or ] curry map ; inline
167 : op-matrix ( domain range quot -- matrix )
168 rot [ [ 2dup ] dip (op-matrix) ] map 2nip ; inline
170 : d-matrix ( domain range -- matrix )
173 : dim-im/ker-d ( domain range -- null/rank )
174 d-matrix null/rank 2array ;
177 : (graded-ker/im-d) ( n seq -- null/rank )
178 #! d: C(n) ---> C(n+1)
179 [ ?nth ] [ [ 1+ ] dip ?nth ] 2bi
182 : graded-ker/im-d ( graded-basis -- seq )
183 [ length ] keep [ (graded-ker/im-d) ] curry map ;
185 : graded-betti ( generators -- seq )
186 basis graded graded-ker/im-d unzip but-last 0 prefix v- ;
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
194 : bigraded-ker/im-d ( bigraded-basis -- seq )
197 [ 2dup ] dip spin (bigraded-ker/im-d)
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
209 : m.m' ( matrix -- matrix' ) dup flip m. ;
210 : m'.m ( matrix -- matrix' ) dup flip swap m. ;
212 : empty-matrix? ( matrix -- ? )
213 [ t ] [ first empty? ] if-empty ;
215 : ?m+ ( m1 m2 -- m3 )
226 : laplacian-matrix ( basis1 basis2 basis3 -- matrix )
227 dupd d-matrix m.m' [ d-matrix m'.m ] dip ?m+ ;
229 : laplacian-betti ( basis1 basis2 basis3 -- n )
230 laplacian-matrix null/rank drop ;
232 : laplacian-kernel ( basis1 basis2 basis3 -- basis )
234 laplacian-matrix dup empty-matrix? [
238 [ [ wedge (alt+) ] 2each ] with-terms
242 : graded-triple ( seq n -- triple )
243 3 [ 1- + ] with map swap [ ?nth ] curry map ;
245 : graded-triples ( seq -- triples )
246 dup length [ graded-triple ] with map ;
248 : graded-laplacian ( generators quot -- seq )
249 [ basis graded graded-triples [ first3 ] ] dip compose map ;
252 : graded-laplacian-betti ( generators -- seq )
253 [ laplacian-betti ] graded-laplacian ;
255 : graded-laplacian-kernel ( generators -- seq )
256 [ laplacian-kernel ] graded-laplacian ;
258 : graded-basis. ( seq -- )
260 "=== Degree " write pprint
261 ": dimension " write dup length .
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 ]
269 [ [ 2 + ] [ 1 - ] [ ] tri* ?nth ?nth ]
273 : bigraded-triples ( grid -- triples )
276 [ 2dup ] dip spin bigraded-triple
280 : bigraded-laplacian ( u-generators z-generators quot -- seq )
281 [ [ basis graded ] bi@ tensor bigraded-triples ] dip
282 [ [ first3 ] prepose map ] curry map ; inline
284 : bigraded-laplacian-betti ( u-generators z-generators -- seq )
285 [ laplacian-betti ] bigraded-laplacian ;
287 : bigraded-laplacian-kernel ( u-generators z-generators -- seq )
288 [ laplacian-kernel ] bigraded-laplacian ;
290 : bigraded-basis. ( seq -- )
292 "=== U-degree " write .
294 " === Z-degree " write pprint
295 ": dimension " write dup length .
296 [ " " write alt. ] each