Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/interfaces/chomp.py
4045 views
1
r"""
2
Interface to CHomP
3
4
CHomP stands for "Computation Homology Program", and is good at
5
computing homology of simplicial complexes, cubical complexes, and
6
chain complexes. It can also compute homomorphisms induced on
7
homology by maps. See the CHomP web page http://chomp.rutgers.edu/
8
for more information.
9
10
AUTHOR:
11
12
- John H. Palmieri
13
"""
14
import re
15
16
_have_chomp = {}
17
def have_chomp(program='homsimpl'):
18
"""
19
Return True if this computer has ``program`` installed.
20
21
The first time it is run, this function caches its result in the
22
variable ``_have_chomp`` -- a dictionary indexed by program name
23
-- and any subsequent time, it just checks the value of the
24
variable.
25
26
This program is used in the routine CHomP.__call__.
27
28
If this computer doesn't have CHomP installed, you may obtain it
29
from http://chomp.rutgers.edu/.
30
31
EXAMPLES::
32
33
sage: from sage.interfaces.chomp import have_chomp
34
sage: have_chomp() # random -- depends on whether CHomP is installed
35
True
36
sage: 'homsimpl' in sage.interfaces.chomp._have_chomp
37
True
38
sage: sage.interfaces.chomp._have_chomp['homsimpl'] == have_chomp()
39
True
40
"""
41
global _have_chomp
42
if program not in _have_chomp:
43
from sage.misc.sage_ostools import have_program
44
_have_chomp[program] = have_program(program)
45
return _have_chomp[program]
46
47
class CHomP:
48
"""
49
Interface to the CHomP package.
50
51
:param program: which CHomP program to use
52
:type program: string
53
:param complex: a simplicial or cubical complex
54
:param subcomplex: a subcomplex of ``complex`` or None (the default)
55
:param base_ring: ring over which to perform computations -- must be `\ZZ` or `\GF{p}`.
56
:type base_ring: ring; optional, default `\ZZ`
57
:param generators: if True, also return list of generators
58
:type generators: boolean; optional, default False
59
:param verbose: if True, print helpful messages as the computation
60
progresses
61
:type verbose: boolean; optional, default False
62
:param extra_opts: options passed directly to ``program``
63
:type extra_opts: string
64
:return: homology groups as a dictionary indexed by dimension
65
66
The programs ``homsimpl``, ``homcubes``, and ``homchain`` are
67
available through this interface. ``homsimpl`` computes the
68
relative or absolute homology groups of simplicial complexes.
69
``homcubes`` computes the relative or absolute homology groups of
70
cubical complexes. ``homchain`` computes the homology groups of
71
chain complexes. For consistency with Sage's other homology
72
computations, the answers produced by ``homsimpl`` and
73
``homcubes`` in the absolute case are converted to reduced
74
homology.
75
76
Note also that CHomP can only compute over the integers or
77
`\GF{p}`. CHomP is fast enough, though, that if you want
78
rational information, you should consider using CHomP with integer
79
coefficients, or with mod `p` coefficients for a sufficiently
80
large `p`, rather than using Sage's built-in homology algorithms.
81
82
See also the documentation for the functions :func:`homchain`,
83
:func:`homcubes`, and :func:`homsimpl` for more examples,
84
including illustrations of some of the optional parameters.
85
86
EXAMPLES::
87
88
sage: from sage.interfaces.chomp import CHomP
89
sage: T = cubical_complexes.Torus()
90
sage: CHomP()('homcubes', T) # optional - CHomP
91
{0: 0, 1: Z x Z, 2: Z}
92
93
Relative homology of a segment relative to its endpoints::
94
95
sage: edge = simplicial_complexes.Simplex(1)
96
sage: ends = edge.n_skeleton(0)
97
sage: CHomP()('homsimpl', edge) # optional - CHomP
98
{0: 0}
99
sage: CHomP()('homsimpl', edge, ends) # optional - CHomP
100
{0: 0, 1: Z}
101
102
Homology of a chain complex::
103
104
sage: C = ChainComplex({3: 2 * identity_matrix(ZZ, 2)}, degree=-1)
105
sage: CHomP()('homchain', C) # optional - CHomP
106
{2: C2 x C2}
107
"""
108
def __repr__(self):
109
"""
110
String representation.
111
112
EXAMPLES::
113
114
sage: from sage.interfaces.chomp import CHomP
115
sage: CHomP().__repr__()
116
'CHomP interface'
117
"""
118
return "CHomP interface"
119
120
def __call__(self, program, complex, subcomplex=None, **kwds):
121
"""
122
Call a CHomP program to compute the homology of a chain
123
complex, simplicial complex, or cubical complex.
124
125
See :class:`CHomP` for full documentation.
126
127
EXAMPLES::
128
129
sage: from sage.interfaces.chomp import CHomP
130
sage: T = cubical_complexes.Torus()
131
sage: CHomP()('homcubes', T) # indirect doctest, optional - CHomP
132
{0: 0, 1: Z x Z, 2: Z}
133
"""
134
from sage.misc.misc import tmp_filename
135
from sage.homology.all import CubicalComplex, cubical_complexes
136
from sage.homology.all import SimplicialComplex, Simplex
137
from sage.homology.cubical_complex import Cube
138
from sage.homology.chain_complex import HomologyGroup, ChainComplex
139
from subprocess import Popen, PIPE
140
from sage.rings.all import QQ, ZZ
141
from sage.modules.all import VectorSpace, vector
142
from sage.combinat.free_module import CombinatorialFreeModule
143
144
if not have_chomp(program):
145
raise OSError, "Program %s not found" % program
146
147
verbose = kwds.get('verbose', False)
148
generators = kwds.get('generators', False)
149
extra_opts = kwds.get('extra_opts', '')
150
base_ring = kwds.get('base_ring', ZZ)
151
152
if extra_opts:
153
extra_opts = extra_opts.split()
154
else:
155
extra_opts = []
156
157
# type of complex:
158
cubical = False
159
simplicial = False
160
chain = False
161
# CHomP seems to have problems with cubical complexes if the
162
# first interval in the first cube defining the complex is
163
# degenerate. So replace the complex X with [0,1] x X.
164
if isinstance(complex, CubicalComplex):
165
cubical = True
166
edge = cubical_complexes.Cube(1)
167
original_complex = complex
168
complex = edge.product(complex)
169
if verbose:
170
print "Cubical complex"
171
elif isinstance(complex, SimplicialComplex):
172
simplicial = True
173
if verbose:
174
print "Simplicial complex"
175
else:
176
chain = True
177
base_ring = kwds.get('base_ring', complex.base_ring())
178
if verbose:
179
print "Chain complex over %s" % base_ring
180
181
if base_ring == QQ:
182
raise ValueError, "CHomP doesn't compute over the rationals, only over Z or F_p."
183
if base_ring.is_prime_field():
184
p = base_ring.characteristic()
185
extra_opts.append('-p%s' % p)
186
mod_p = True
187
else:
188
mod_p = False
189
190
#
191
# complex
192
#
193
try:
194
data = complex._chomp_repr_()
195
except AttributeError:
196
raise AttributeError, "Complex can not be converted to use with CHomP."
197
198
datafile = tmp_filename()
199
f = open(datafile, 'w')
200
f.write(data)
201
f.close()
202
203
#
204
# subcomplex
205
#
206
if subcomplex is None:
207
if cubical:
208
subcomplex = CubicalComplex([complex.n_cells(0)[0]])
209
elif simplicial:
210
m = re.search(r'\(([^,]*),', data)
211
v = int(m.group(1))
212
subcomplex = SimplicialComplex([v], [[v]])
213
else:
214
# replace subcomplex with [0,1] x subcomplex.
215
if cubical:
216
subcomplex = edge.product(subcomplex)
217
#
218
# generators
219
#
220
if generators:
221
genfile = tmp_filename()
222
extra_opts.append('-g%s' % genfile)
223
224
#
225
# call program
226
#
227
if subcomplex is not None:
228
try:
229
sub = subcomplex._chomp_repr_()
230
except AttributeError:
231
raise AttributeError, "Subcomplex can not be converted to use with CHomP."
232
subfile = tmp_filename()
233
f = open(subfile, 'w')
234
f.write(sub)
235
f.close()
236
else:
237
subfile = ''
238
if verbose:
239
print "Popen called with arguments",
240
print [program, datafile, subfile] + extra_opts
241
print
242
print "CHomP output:"
243
print
244
# output = Popen([program, datafile, subfile, extra_opts],
245
cmd = [program, datafile]
246
if subfile:
247
cmd.append(subfile)
248
if extra_opts:
249
cmd.extend(extra_opts)
250
output = Popen(cmd, stdout=PIPE).communicate()[0]
251
if verbose:
252
print output
253
print "End of CHomP output"
254
print
255
if generators:
256
gens = open(genfile, 'r').read()
257
if verbose:
258
print "Generators:"
259
print gens
260
261
#
262
# process output
263
#
264
# output contains substrings of one of the forms
265
# "H_1 = Z", "H_1 = Z_2 + Z", "H_1 = Z_2 + Z^2",
266
# "H_1 = Z + Z_2 + Z"
267
if output.find('trivial') != -1:
268
if mod_p:
269
return {0: VectorSpace(base_ring, 0)}
270
else:
271
return {0: HomologyGroup(0)}
272
d = {}
273
h = re.compile("^H_([0-9]*) = (.*)$", re.M)
274
tors = re.compile("Z_([0-9]*)")
275
#
276
# homology groups
277
#
278
for m in h.finditer(output):
279
if verbose:
280
print m.groups()
281
# dim is the dimension of the homology group
282
dim = int(m.group(1))
283
# hom_str is the right side of the equation "H_n = Z^r + Z_k + ..."
284
hom_str = m.group(2)
285
# need to read off number of summands and their invariants
286
if hom_str.find("0") == 0:
287
if mod_p:
288
hom = VectorSpace(base_ring, 0)
289
else:
290
hom = HomologyGroup(0)
291
else:
292
rk = 0
293
if hom_str.find("^") != -1:
294
rk_srch = re.search(r'\^([0-9]*)\s?', hom_str)
295
rk = int(rk_srch.group(1))
296
rk += len(re.findall("(Z$)|(Z\s)", hom_str))
297
if mod_p:
298
rk = rk if rk != 0 else 1
299
if verbose:
300
print "dimension = %s, rank of homology = %s" % (dim, rk)
301
hom = VectorSpace(base_ring, rk)
302
else:
303
n = rk
304
invts = []
305
for t in tors.finditer(hom_str):
306
n += 1
307
invts.append(int(t.group(1)))
308
for i in range(rk):
309
invts.append(0)
310
if verbose:
311
print "dimension = %s, number of factors = %s, invariants = %s" %(dim, n, invts)
312
hom = HomologyGroup(n, invts)
313
314
#
315
# generators
316
#
317
if generators:
318
if cubical:
319
g = process_generators_cubical(gens, dim)
320
if verbose:
321
print "raw generators: %s" % g
322
if g:
323
module = CombinatorialFreeModule(base_ring, original_complex.n_cells(dim))
324
module._repr_option_bracket = True
325
module._prefix = ""
326
basis = module.basis()
327
output = []
328
for x in g:
329
v = module(0)
330
for term in x:
331
v += term[0] * basis[term[1]]
332
output.append(v)
333
g = output
334
elif simplicial:
335
g = process_generators_simplicial(gens, dim, complex)
336
if verbose:
337
print "raw generators: %s" % gens
338
if g:
339
module = CombinatorialFreeModule(base_ring, complex.n_cells(dim))
340
module._repr_option_bracket = False
341
module._prefix = ""
342
basis = module.basis()
343
output = []
344
for x in g:
345
v = module(0)
346
for term in x:
347
if complex._numeric:
348
v += term[0] * basis[term[1]]
349
else:
350
# if not numeric: have to
351
# translate vertices from numbers
352
# to their original form.
353
translate = dict([(b,a) for (a,b) in complex._numeric_translation])
354
simplex = Simplex([translate[a] for a in tuple(term[1])])
355
v += term[0] * basis[simplex]
356
output.append(v)
357
g = output
358
elif chain:
359
g = process_generators_chain(gens, dim, base_ring)
360
if verbose:
361
print "raw generators: %s" % gens
362
if g:
363
d[dim] = (hom, g)
364
else:
365
d[dim] = hom
366
else:
367
d[dim] = hom
368
369
if chain:
370
new_d = {}
371
bottom = min(complex.differential())
372
top = max(complex.differential())
373
for dim in d:
374
if complex._degree == -1: # chain complex
375
new_dim = bottom + dim
376
else: # cochain complex
377
new_dim = top - dim
378
if isinstance(d[dim], tuple):
379
# generators included.
380
group = d[dim][0]
381
gens = d[dim][1]
382
new_gens = []
383
dimension = complex.differential(new_dim).ncols()
384
# make sure that each vector is embedded in the
385
# correct ambient space: pad with a zero if
386
# necessary.
387
for v in gens:
388
v_dict = v.dict()
389
if dimension - 1 not in v.dict():
390
v_dict[dimension - 1] = 0
391
new_gens.append(vector(base_ring, v_dict))
392
else:
393
new_gens.append(v)
394
new_d[new_dim] = (group, new_gens)
395
else:
396
new_d[new_dim] = d[dim]
397
d = new_d
398
return d
399
400
def help(self, program):
401
"""
402
Print a help message for ``program``, a program from the CHomP suite.
403
404
:param program: which CHomP program to use
405
:type program: string
406
:return: nothing -- just print a message
407
408
EXAMPLES::
409
410
sage: from sage.interfaces.chomp import CHomP
411
sage: CHomP().help('homcubes') # optional - CHomP
412
HOMCUBES, ver. ... Copyright (C) ... by Pawel Pilarczyk...
413
"""
414
from subprocess import Popen, PIPE
415
print Popen([program, '-h'], stdout=PIPE).communicate()[0]
416
417
def homsimpl(complex=None, subcomplex=None, **kwds):
418
r"""
419
Compute the homology of a simplicial complex using the CHomP
420
program ``homsimpl``. If the argument ``subcomplex`` is present,
421
compute homology of ``complex`` relative to ``subcomplex``.
422
423
:param complex: a simplicial complex
424
:param subcomplex: a subcomplex of ``complex`` or None (the default)
425
:param base_ring: ring over which to perform computations -- must be `\ZZ` or `\GF{p}`.
426
:type base_ring: ring; optional, default `\ZZ`
427
:param generators: if True, also return list of generators
428
:type generators: boolean; optional, default False
429
:param verbose: if True, print helpful messages as the computation
430
progresses
431
:type verbose: boolean; optional, default False
432
:param help: if True, just print a help message and exit
433
:type help: boolean; optional, default False
434
:param extra_opts: options passed directly to ``program``
435
:type extra_opts: string
436
:return: homology groups as a dictionary indexed by dimension
437
438
EXAMPLES::
439
440
sage: from sage.interfaces.chomp import homsimpl
441
sage: T = simplicial_complexes.Torus()
442
sage: M8 = simplicial_complexes.MooreSpace(8)
443
sage: M4 = simplicial_complexes.MooreSpace(4)
444
sage: X = T.disjoint_union(T).disjoint_union(T).disjoint_union(M8).disjoint_union(M4)
445
sage: homsimpl(X)[1] # optional - CHomP
446
Z^6 x C4 x C8
447
448
Relative homology::
449
450
sage: S = simplicial_complexes.Simplex(3)
451
sage: bdry = S.n_skeleton(2)
452
sage: homsimpl(S, bdry)[3] # optional - CHomP
453
Z
454
455
Generators: these are given as a list after the homology group.
456
Each generator is specified as a linear combination of simplices::
457
458
sage: homsimpl(S, bdry, generators=True)[3] # optional - CHomP
459
(Z, [(0, 1, 2, 3)])
460
461
sage: homsimpl(simplicial_complexes.Sphere(1), generators=True) # optional - CHomP
462
{0: 0, 1: (Z, [(0, 1) - (0, 2) + (1, 2)])}
463
464
TESTS:
465
466
Generators for a simplicial complex whose vertices are not integers::
467
468
sage: S1 = simplicial_complexes.Sphere(1)
469
sage: homsimpl(S1.join(S1), generators=True, base_ring=GF(2))[3][1] # optional - CHomP
470
[('L0', 'L1', 'R0', 'R1') + ('L0', 'L1', 'R0', 'R2') + ('L0', 'L1', 'R1', 'R2') + ('L0', 'L2', 'R0', 'R1') + ('L0', 'L2', 'R0', 'R2') + ('L0', 'L2', 'R1', 'R2') + ('L1', 'L2', 'R0', 'R1') + ('L1', 'L2', 'R0', 'R2') + ('L1', 'L2', 'R1', 'R2')]
471
"""
472
from sage.homology.all import SimplicialComplex
473
help = kwds.get('help', False)
474
if help:
475
return CHomP().help('homsimpl')
476
# Check types here, because if you pass homsimpl a cubical
477
# complex, it tries to compute its homology as if it were a
478
# simplicial complex and gets terribly wrong answers.
479
if (isinstance(complex, SimplicialComplex)
480
and (subcomplex is None or isinstance(subcomplex, SimplicialComplex))):
481
return CHomP()('homsimpl', complex, subcomplex=subcomplex, **kwds)
482
else:
483
raise TypeError, "Complex and/or subcomplex are not simplicial complexes."
484
485
def homcubes(complex=None, subcomplex=None, **kwds):
486
r"""
487
Compute the homology of a cubical complex using the CHomP program
488
``homcubes``. If the argument ``subcomplex`` is present, compute
489
homology of ``complex`` relative to ``subcomplex``.
490
491
:param complex: a cubical complex
492
:param subcomplex: a subcomplex of ``complex`` or None (the default)
493
:param generators: if True, also return list of generators
494
:type generators: boolean; optional, default False
495
:param verbose: if True, print helpful messages as the computation progresses
496
:type verbose: boolean; optional, default False
497
:param help: if True, just print a help message and exit
498
:type help: boolean; optional, default False
499
:param extra_opts: options passed directly to ``homcubes``
500
:type extra_opts: string
501
:return: homology groups as a dictionary indexed by dimension
502
503
EXAMPLES::
504
505
sage: from sage.interfaces.chomp import homcubes
506
sage: S = cubical_complexes.Sphere(3)
507
sage: homcubes(S)[3] # optional - CHomP
508
Z
509
510
Relative homology::
511
512
sage: C3 = cubical_complexes.Cube(3)
513
sage: bdry = C3.n_skeleton(2)
514
sage: homcubes(C3, bdry) # optional - CHomP
515
{0: 0, 1: 0, 2: 0, 3: Z}
516
517
Generators: these are given as a list after the homology group.
518
Each generator is specified as a linear combination of cubes::
519
520
sage: homcubes(cubical_complexes.Sphere(1), generators=True, base_ring=GF(2))[1][1] # optional - CHomP
521
[[[1,1] x [0,1]] + [[0,1] x [1,1]] + [[0,1] x [0,0]] + [[0,0] x [0,1]]]
522
"""
523
from sage.homology.all import CubicalComplex
524
help = kwds.get('help', False)
525
if help:
526
return CHomP().help('homcubes')
527
# Type-checking is here for the same reason as in homsimpl above.
528
if (isinstance(complex, CubicalComplex)
529
and (subcomplex is None or isinstance(subcomplex, CubicalComplex))):
530
return CHomP()('homcubes', complex, subcomplex=subcomplex, **kwds)
531
else:
532
raise TypeError, "Complex and/or subcomplex are not cubical complexes."
533
534
def homchain(complex=None, **kwds):
535
r"""
536
Compute the homology of a chain complex using the CHomP program
537
``homchain``.
538
539
:param complex: a chain complex
540
:param generators: if True, also return list of generators
541
:type generators: boolean; optional, default False
542
:param verbose: if True, print helpful messages as the computation progresses
543
:type verbose: boolean; optional, default False
544
:param help: if True, just print a help message and exit
545
:type help: boolean; optional, default False
546
:param extra_opts: options passed directly to ``homchain``
547
:type extra_opts: string
548
:return: homology groups as a dictionary indexed by dimension
549
550
EXAMPLES::
551
552
sage: from sage.interfaces.chomp import homchain
553
sage: C = cubical_complexes.Sphere(3).chain_complex()
554
sage: homchain(C)[3] # optional - CHomP
555
Z
556
557
Generators: these are given as a list after the homology group.
558
Each generator is specified as a cycle, an element in the
559
appropriate free module over the base ring::
560
561
sage: C2 = delta_complexes.Sphere(2).chain_complex()
562
sage: homchain(C2, generators=True)[2] # optional - CHomP
563
(Z, [(1, -1)])
564
sage: homchain(C2, generators=True, base_ring=GF(2))[2] # optional - CHomP
565
(Vector space of dimension 1 over Finite Field of size 2, [(1, 1)])
566
567
TESTS:
568
569
Chain complexes concentrated in negative dimensions, cochain complexes, etc.::
570
571
sage: C = ChainComplex({-5: 4 * identity_matrix(ZZ, 2)}, degree=-1)
572
sage: homchain(C) # optional - CHomP
573
{-6: C4 x C4}
574
sage: C = ChainComplex({-5: 4 * identity_matrix(ZZ, 2)}, degree=1)
575
sage: homchain(C, generators=True) # optional - CHomP
576
{-4: (C4 x C4, [(1, 0), (0, 1)])}
577
"""
578
from sage.homology.all import ChainComplex
579
help = kwds.get('help', False)
580
if help:
581
return CHomP().help('homchain')
582
# Type-checking just in case.
583
if isinstance(complex, ChainComplex):
584
return CHomP()('homchain', complex, **kwds)
585
else:
586
raise TypeError, "Complex is not a chain complex."
587
588
def process_generators_cubical(gen_string, dim):
589
r"""
590
Process CHomP generator information for cubical complexes.
591
592
:param gen_string: generator output from CHomP
593
:type gen_string: string
594
:param dim: dimension in which to find generators
595
:type dim: integer
596
:return: list of generators in each dimension, as described below
597
598
``gen_string`` has the form ::
599
600
The 2 generators of H_1 follow:
601
generator 1
602
-1 * [(0,0,0,0,0)(0,0,0,0,1)]
603
1 * [(0,0,0,0,0)(0,0,1,0,0)]
604
...
605
generator 2
606
-1 * [(0,1,0,1,1)(1,1,0,1,1)]
607
-1 * [(0,1,0,0,1)(0,1,0,1,1)]
608
...
609
610
Each line consists of a coefficient multiplied by a cube; the cube
611
is specified by its "bottom left" and "upper right" corners.
612
613
For technical reasons, we remove the first coordinate of each
614
tuple, and using regular expressions, the remaining parts get
615
converted from a string to a pair ``(coefficient, Cube)``, with
616
the cube represented as a product of tuples. For example, the
617
first line in "generator 1" gets turned into ::
618
619
(-1, [0,0] x [0,0] x [0,0] x [0,1])
620
621
representing an element in the free abelian group with basis given
622
by cubes. Each generator is a list of such pairs, representing
623
the sum of such elements. These are reassembled in
624
:meth:`CHomP.__call__` to actual elements in the free module
625
generated by the cubes of the cubical complex in the appropriate
626
dimension.
627
628
Therefore the return value is a list of lists of pairs, one list
629
of pairs for each generator.
630
631
EXAMPLES::
632
633
sage: from sage.interfaces.chomp import process_generators_cubical
634
sage: s = "The 2 generators of H_1 follow:\ngenerator 1:\n-1 * [(0,0,0,0,0)(0,0,0,0,1)]\n1 * [(0,0,0,0,0)(0,0,1,0,0)]"
635
sage: process_generators_cubical(s, 1)
636
[[(-1, [0,0] x [0,0] x [0,0] x [0,1]), (1, [0,0] x [0,1] x [0,0] x [0,0])]]
637
sage: len(process_generators_cubical(s, 1)) # only one generator
638
1
639
"""
640
from sage.homology.cubical_complex import Cube
641
# each dim in gen_string starts with "The generator for
642
# H_3 follows:". So search for "T" to find the
643
# end of the current list of generators.
644
g_srch = re.compile(r'H_%s\sfollow[s]?:\n([^T]*)(?:T|$)' % dim)
645
g = g_srch.search(gen_string)
646
output = []
647
cubes_list = []
648
if g:
649
g = g.group(1)
650
if g:
651
# project g to one end of the cylinder [0,1] x complex:
652
#
653
# drop the first coordinate and eliminate duplicates, at least
654
# in positive dimensions, drop any line containing a
655
# degenerate cube
656
g = re.sub('\([01],', '(', g)
657
if dim > 0:
658
lines = g.splitlines()
659
newlines = []
660
for l in lines:
661
x = re.findall(r'(\([0-9,]*\))', l)
662
if x:
663
left, right = x
664
left = [int(a) for a in left.strip('()').split(',')]
665
right = [int(a) for a in right.strip('()').split(',')]
666
if sum([x-y for (x,y) in zip(right, left)]) == dim:
667
newlines.append(l)
668
else: # line like "generator 2"
669
newlines.append(l)
670
g = newlines
671
cubes_list = []
672
for l in g:
673
cubes = re.search(r'([+-]?)\s?([0-9]+)?\s?[*]?\s?\[\(([-0-9,]+)\)\(([-0-9,]+)\)\]', l)
674
if cubes:
675
if cubes.group(1) and re.search("-", cubes.group(1)):
676
sign = -1
677
else:
678
sign = 1
679
if cubes.group(2) and len(cubes.group(2)) > 0:
680
coeff = sign * int(cubes.group(2))
681
else:
682
coeff = sign * 1
683
cube1 = [int(a) for a in cubes.group(3).split(',')]
684
cube2 = [int(a) for a in cubes.group(4).split(',')]
685
cube = Cube(zip(cube1, cube2))
686
cubes_list.append((coeff, cube))
687
else:
688
if cubes_list:
689
output.append(cubes_list)
690
cubes_list = []
691
if cubes_list:
692
output.append(cubes_list)
693
return output
694
else:
695
return None
696
697
def process_generators_simplicial(gen_string, dim, complex):
698
r"""
699
Process CHomP generator information for simplicial complexes.
700
701
:param gen_string: generator output from CHomP
702
:type gen_string: string
703
:param dim: dimension in which to find generators
704
:type dim: integer
705
:param complex: simplicial complex under consideration
706
:return: list of generators in each dimension, as described below
707
708
``gen_string`` has the form ::
709
710
The 2 generators of H_1 follow:
711
generator 1
712
-1 * (1,6)
713
1 * (1,4)
714
...
715
generator 2
716
-1 * (1,6)
717
1 * (1,4)
718
...
719
720
where each line contains a coefficient and a simplex. Each line
721
is converted, using regular expressions, from a string to a pair
722
``(coefficient, Simplex)``, like ::
723
724
(-1, (1,6))
725
726
representing an element in the free abelian group with basis given
727
by simplices. Each generator is a list of such pairs,
728
representing the sum of such elements. These are reassembled in
729
:meth:`CHomP.__call__` to actual elements in the free module
730
generated by the simplices of the simplicial complex in the
731
appropriate dimension.
732
733
Therefore the return value is a list of lists of pairs, one list
734
of pairs for each generator.
735
736
EXAMPLES::
737
738
sage: from sage.interfaces.chomp import process_generators_simplicial
739
sage: s = "The 2 generators of H_1 follow:\ngenerator 1:\n-1 * (1,6)\n1 * (1,4)"
740
sage: process_generators_simplicial(s, 1, simplicial_complexes.Torus())
741
[[(-1, (1, 6)), (1, (1, 4))]]
742
"""
743
from sage.homology.all import Simplex
744
# each dim in gen_string starts with "The generator for H_3
745
# follows:". So search for "T" to find the end of the current
746
# list of generators.
747
g_srch = re.compile(r'H_%s\sfollow[s]?:\n([^T]*)(?:T|$)' % dim)
748
g = g_srch.search(gen_string)
749
output = []
750
simplex_list = []
751
if g:
752
g = g.group(1)
753
if g:
754
if not complex._numeric:
755
# need to rename the vertices
756
translate = dict([(b,a) for (a,b) in complex._numeric_translation])
757
lines = g.splitlines()
758
newlines = []
759
for l in lines:
760
simplex = re.search(r'([+-]?)\s?([0-9]+)?\s?[*]?\s?(\([0-9,]*\))', l)
761
if simplex:
762
if simplex.group(1) and re.search("-", simplex.group(1)):
763
sign = -1
764
else:
765
sign = 1
766
if simplex.group(2) and len(simplex.group(2)) > 0:
767
coeff = sign * int(simplex.group(2))
768
else:
769
coeff = sign * 1
770
simplex = Simplex([int(a) for a in simplex.group(3).strip('()').split(',')])
771
simplex_list.append((coeff, simplex))
772
else:
773
if simplex_list:
774
output.append(simplex_list)
775
simplex_list = []
776
if simplex_list:
777
output.append(simplex_list)
778
return output
779
else:
780
return None
781
782
def process_generators_chain(gen_string, dim, base_ring=None):
783
r"""
784
Process CHomP generator information for simplicial complexes.
785
786
:param gen_string: generator output from CHomP
787
:type gen_string: string
788
:param dim: dimension in which to find generators
789
:type dim: integer
790
:param base_ring: base ring over which to do the computations
791
:type base_ring: optional, default ZZ
792
:return: list of generators in each dimension, as described below
793
794
``gen_string`` has the form ::
795
796
[H_0]
797
a1
798
799
[H_1]
800
a2
801
a3
802
803
[H_2]
804
a1 - a2
805
806
For each homology group, each line lists a homology generator as a
807
linear combination of generators ``ai`` of the group of chains in
808
the appropriate dimension. The elements ``ai`` are indexed
809
starting with `i=1`. Each generator is converted, using regular
810
expressions, from a string to a vector (an element in the free
811
module over ``base_ring``), with ``ai`` representing the unit
812
vector in coordinate `i-1`. For example, the string ``a1 - a2``
813
gets converted to the vector ``(1, -1)``.
814
815
Therefore the return value is a list of vectors.
816
817
EXAMPLES::
818
819
sage: from sage.interfaces.chomp import process_generators_chain
820
sage: s = "[H_0]\na1\n\n[H_1]\na2\na3\n"
821
sage: process_generators_chain(s, 1)
822
[(0, 1), (0, 0, 1)]
823
sage: s = "[H_0]\na1\n\n[H_1]\n5 * a2 - a1\na3\n"
824
sage: process_generators_chain(s, 1, base_ring=ZZ)
825
[(-1, 5), (0, 0, 1)]
826
sage: process_generators_chain(s, 1, base_ring=GF(2))
827
[(1, 1), (0, 0, 1)]
828
"""
829
from sage.modules.all import vector
830
from sage.rings.all import ZZ
831
if base_ring is None:
832
base_ring = ZZ
833
# each dim in gens starts with a string like
834
# "[H_3]". So search for "[" to find the end of
835
# the current list of generators.
836
g_srch = re.compile(r'\[H_%s\]\n([^]]*)(?:\[|$)' % dim)
837
g = g_srch.search(gen_string)
838
if g:
839
g = g.group(1)
840
if g:
841
# each line in the string g is a linear
842
# combination of things like "a2", "a31", etc.
843
# indexing on the a's starts at 1.
844
lines = g.splitlines()
845
new_gens = []
846
for l in lines:
847
gen = re.compile(r"([+-]?)\s?([0-9]+)?\s?[*]?\s?a([0-9]*)(?:\s|$)")
848
v = {}
849
for term in gen.finditer(l):
850
if term.group(1) and re.search("-", term.group(1)):
851
sign = -1
852
else:
853
sign = 1
854
if term.group(2) and len(term.group(2)) > 0:
855
coeff = sign * int(term.group(2))
856
else:
857
coeff = sign * 1
858
idx = int(term.group(3))
859
v[idx-1] = coeff
860
if v:
861
new_gens.append(vector(base_ring, v))
862
g = new_gens
863
return g
864
865