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 symbols 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 ) >r ((d)) r> wedge ;
115 : x.dy ( x y -- vec ) (d) wedge -1 alt*n ;
117 : (d) ( product -- value )
119 [ drop H{ } ] [ unclip swap [ x.dy ] 2keep dx.y alt+ ] if ;
121 : linear-op ( vec quot -- vec )
124 -rot >r swap call r> alt*n (alt+)
126 ] with-terms ; inline
129 >alt [ (d) ] linear-op ;
132 : (interior) ( y basis-elt -- i_y[basis-elt] )
134 -rot remove associate
139 : interior ( x y -- i_y[x] )
141 swap >alt [ dupd (interior) ] linear-op nip ;
144 : graded ( seq -- seq )
145 dup 0 [ length max ] reduce 1+ [ V{ } clone ] replicate
146 [ dup length pick nth push ] reduce ;
148 : nth-basis-elt ( generators n -- elt )
150 3dup bit? [ nth ] [ 2drop f ] if
153 : basis ( generators -- seq )
154 natural-sort dup length 2^ [ nth-basis-elt ] with map ;
156 : (tensor) ( seq1 seq2 -- seq )
158 [ prepend natural-sort ] curry map
161 : tensor ( graded-basis1 graded-basis2 -- bigraded-basis )
162 [ [ swap (tensor) ] curry map ] with map ;
164 ! Computing cohomology
165 : (op-matrix) ( range quot basis-elt -- row )
166 swap call [ at 0 or ] curry map ; inline
168 : op-matrix ( domain range quot -- matrix )
169 rot [ >r 2dup r> (op-matrix) ] map 2nip ; inline
171 : d-matrix ( domain range -- matrix )
174 : dim-im/ker-d ( domain range -- null/rank )
175 d-matrix null/rank 2array ;
178 : (graded-ker/im-d) ( n seq -- null/rank )
179 #! d: C(n) ---> C(n+1)
180 [ ?nth ] 2keep >r 1+ r> ?nth
183 : graded-ker/im-d ( graded-basis -- seq )
184 [ length ] keep [ (graded-ker/im-d) ] curry map ;
186 : graded-betti ( generators -- seq )
187 basis graded graded-ker/im-d flip first2 but-last 0 prefix v- ;
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
195 : bigraded-ker/im-d ( bigraded-basis -- seq )
198 >r 2dup r> spin (bigraded-ker/im-d)
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
210 : m.m' ( matrix -- matrix' ) dup flip m. ;
211 : m'.m ( matrix -- matrix' ) dup flip swap m. ;
213 : empty-matrix? ( matrix -- ? )
214 dup empty? [ drop t ] [ first empty? ] if ;
216 : ?m+ ( m1 m2 -- m3 )
227 : laplacian-matrix ( basis1 basis2 basis3 -- matrix )
228 dupd d-matrix m.m' >r d-matrix m'.m r> ?m+ ;
230 : laplacian-betti ( basis1 basis2 basis3 -- n )
231 laplacian-matrix null/rank drop ;
233 : laplacian-kernel ( basis1 basis2 basis3 -- basis )
235 laplacian-matrix dup empty-matrix? [
239 [ [ wedge (alt+) ] 2each ] with-terms
243 : graded-triple ( seq n -- triple )
244 3 [ 1- + ] with map swap [ ?nth ] curry map ;
246 : graded-triples ( seq -- triples )
247 dup length [ graded-triple ] with map ;
249 : graded-laplacian ( generators quot -- seq )
250 >r basis graded graded-triples [ first3 ] r> compose map ;
253 : graded-laplacian-betti ( generators -- seq )
254 [ laplacian-betti ] graded-laplacian ;
256 : graded-laplacian-kernel ( generators -- seq )
257 [ laplacian-kernel ] graded-laplacian ;
259 : graded-basis. ( seq -- )
261 "=== Degree " write pprint
262 ": dimension " write dup length .
266 : bigraded-triple ( u-deg z-deg bigraded-basis -- triple )
267 #! d: C(u,z) ---> C(u+2,z-1)
268 [ [ 2 - ] [ 1 + ] [ ] tri* ?nth ?nth ]
270 [ [ 2 + ] [ 1 - ] [ ] tri* ?nth ?nth ]
274 : bigraded-triples ( grid -- triples )
277 >r 2dup r> spin bigraded-triple
281 : bigraded-laplacian ( u-generators z-generators quot -- seq )
282 >r [ basis graded ] bi@ tensor bigraded-triples r>
283 [ [ first3 ] prepose map ] curry map ; inline
285 : bigraded-laplacian-betti ( u-generators z-generators -- seq )
286 [ laplacian-betti ] bigraded-laplacian ;
288 : bigraded-laplacian-kernel ( u-generators z-generators -- seq )
289 [ laplacian-kernel ] bigraded-laplacian ;
291 : bigraded-basis. ( seq -- )
293 "=== U-degree " write .
295 " === Z-degree " write pprint
296 ": dimension " write dup length .
297 [ " " write alt. ] each