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 locals 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 [ (op-matrix) ] with with map ; 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 ( basis -- seq )
195 basis length iota [| z |
196 basis first length iota [| u |
197 u z basis (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 )
233 basis1 basis2 basis3 laplacian-matrix :> lap
234 lap empty-matrix? [ f ] [
236 basis2 x [ [ wedge (alt+) ] 2each ] with-terms
240 : graded-triple ( seq n -- triple )
241 3 [ 1 - + ] with map swap [ ?nth ] curry map ;
243 : graded-triples ( seq -- triples )
244 dup length [ graded-triple ] with map ;
246 : graded-laplacian ( generators quot -- seq )
247 [ basis graded graded-triples [ first3 ] ] dip compose map ;
250 : graded-laplacian-betti ( generators -- seq )
251 [ laplacian-betti ] graded-laplacian ;
253 : graded-laplacian-kernel ( generators -- seq )
254 [ laplacian-kernel ] graded-laplacian ;
256 : graded-basis. ( seq -- )
258 "=== Degree " write pprint
259 ": dimension " write dup length .
263 : bigraded-triple ( u-deg z-deg bigraded-basis -- triple )
264 #! d: C(u,z) ---> C(u+2,z-1)
265 [ [ 2 - ] [ 1 + ] [ ] tri* ?nth ?nth ]
267 [ [ 2 + ] [ 1 - ] [ ] tri* ?nth ?nth ]
271 :: bigraded-triples ( grid -- triples )
273 grid first length [| u |
274 u z grid bigraded-triple
278 : bigraded-laplacian ( u-generators z-generators quot -- seq )
279 [ [ basis graded ] bi@ tensor bigraded-triples ] dip
280 [ [ first3 ] prepose map ] curry map ; inline
282 : bigraded-laplacian-betti ( u-generators z-generators -- seq )
283 [ laplacian-betti ] bigraded-laplacian ;
285 : bigraded-laplacian-kernel ( u-generators z-generators -- seq )
286 [ laplacian-kernel ] bigraded-laplacian ;
288 : bigraded-basis. ( seq -- )
290 "=== U-degree " write .
292 " === Z-degree " write pprint
293 ": dimension " write dup length .
294 [ " " write alt. ] each