Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
r""" Interface to CHomP
CHomP stands for "Computation Homology Program", and is good at computing homology of simplicial complexes, cubical complexes, and chain complexes. It can also compute homomorphisms induced on homology by maps. See the CHomP web page http://chomp.rutgers.edu/ for more information.
AUTHOR:
- John H. Palmieri """
""" Return True if this computer has ``program`` installed.
The first time it is run, this function caches its result in the variable ``_have_chomp`` -- a dictionary indexed by program name -- and any subsequent time, it just checks the value of the variable.
This program is used in the routine CHomP.__call__.
If this computer doesn't have CHomP installed, you may obtain it from http://chomp.rutgers.edu/.
EXAMPLES::
sage: from sage.interfaces.chomp import have_chomp sage: have_chomp() # random -- depends on whether CHomP is installed True sage: 'homsimpl' in sage.interfaces.chomp._have_chomp True sage: sage.interfaces.chomp._have_chomp['homsimpl'] == have_chomp() True """ global _have_chomp
""" Interface to the CHomP package.
:param program: which CHomP program to use :type program: string :param complex: a simplicial or cubical complex :param subcomplex: a subcomplex of ``complex`` or None (the default) :param base_ring: ring over which to perform computations -- must be `\ZZ` or `\GF{p}`. :type base_ring: ring; optional, default `\ZZ` :param generators: if True, also return list of generators :type generators: boolean; optional, default False :param verbose: if True, print helpful messages as the computation progresses :type verbose: boolean; optional, default False :param extra_opts: options passed directly to ``program`` :type extra_opts: string :return: homology groups as a dictionary indexed by dimension
The programs ``homsimpl``, ``homcubes``, and ``homchain`` are available through this interface. ``homsimpl`` computes the relative or absolute homology groups of simplicial complexes. ``homcubes`` computes the relative or absolute homology groups of cubical complexes. ``homchain`` computes the homology groups of chain complexes. For consistency with Sage's other homology computations, the answers produced by ``homsimpl`` and ``homcubes`` in the absolute case are converted to reduced homology.
Note also that CHomP can only compute over the integers or `\GF{p}`. CHomP is fast enough, though, that if you want rational information, you should consider using CHomP with integer coefficients, or with mod `p` coefficients for a sufficiently large `p`, rather than using Sage's built-in homology algorithms.
See also the documentation for the functions :func:`homchain`, :func:`homcubes`, and :func:`homsimpl` for more examples, including illustrations of some of the optional parameters.
EXAMPLES::
sage: from sage.interfaces.chomp import CHomP sage: T = cubical_complexes.Torus() sage: CHomP()('homcubes', T) # optional - CHomP {0: 0, 1: Z x Z, 2: Z}
Relative homology of a segment relative to its endpoints::
sage: edge = simplicial_complexes.Simplex(1) sage: ends = edge.n_skeleton(0) sage: CHomP()('homsimpl', edge) # optional - CHomP {0: 0} sage: CHomP()('homsimpl', edge, ends) # optional - CHomP {0: 0, 1: Z}
Homology of a chain complex::
sage: C = ChainComplex({3: 2 * identity_matrix(ZZ, 2)}, degree=-1) sage: CHomP()('homchain', C) # optional - CHomP {2: C2 x C2} """ """ String representation.
EXAMPLES::
sage: from sage.interfaces.chomp import CHomP sage: CHomP() # indirect doctest CHomP interface """
""" Call a CHomP program to compute the homology of a chain complex, simplicial complex, or cubical complex.
See :class:`CHomP` for full documentation.
EXAMPLES::
sage: from sage.interfaces.chomp import CHomP sage: T = cubical_complexes.Torus() sage: CHomP()('homcubes', T) # indirect doctest, optional - CHomP {0: 0, 1: Z x Z, 2: Z} """ from sage.misc.temporary_file import tmp_filename from sage.homology.all import CubicalComplex, cubical_complexes from sage.homology.all import SimplicialComplex, Simplex from sage.homology.chain_complex import HomologyGroup from subprocess import Popen, PIPE from sage.rings.all import QQ, ZZ from sage.modules.all import VectorSpace, vector from sage.combinat.free_module import CombinatorialFreeModule
if not have_chomp(program): raise OSError("Program %s not found" % program)
verbose = kwds.get('verbose', False) generators = kwds.get('generators', False) extra_opts = kwds.get('extra_opts', '') base_ring = kwds.get('base_ring', ZZ)
if extra_opts: extra_opts = extra_opts.split() else: extra_opts = []
# type of complex: cubical = False simplicial = False chain = False # CHomP seems to have problems with cubical complexes if the # first interval in the first cube defining the complex is # degenerate. So replace the complex X with [0,1] x X. if isinstance(complex, CubicalComplex): cubical = True edge = cubical_complexes.Cube(1) original_complex = complex complex = edge.product(complex) if verbose: print("Cubical complex") elif isinstance(complex, SimplicialComplex): simplicial = True if verbose: print("Simplicial complex") else: chain = True base_ring = kwds.get('base_ring', complex.base_ring()) if verbose: print("Chain complex over %s" % base_ring)
if base_ring == QQ: raise ValueError("CHomP doesn't compute over the rationals, only over Z or F_p.") if base_ring.is_prime_field(): p = base_ring.characteristic() extra_opts.append('-p%s' % p) mod_p = True else: mod_p = False
# # complex # try: data = complex._chomp_repr_() except AttributeError: raise AttributeError("Complex can not be converted to use with CHomP.")
datafile = tmp_filename() f = open(datafile, 'w') f.write(data) f.close()
# # subcomplex # if subcomplex is None: if cubical: subcomplex = CubicalComplex([complex.n_cells(0)[0]]) elif simplicial: m = re.search(r'\(([^,]*),', data) v = int(m.group(1)) subcomplex = SimplicialComplex([[v]]) else: # replace subcomplex with [0,1] x subcomplex. if cubical: subcomplex = edge.product(subcomplex) # # generators # if generators: genfile = tmp_filename() extra_opts.append('-g%s' % genfile)
# # call program # if subcomplex is not None: try: sub = subcomplex._chomp_repr_() except AttributeError: raise AttributeError("Subcomplex can not be converted to use with CHomP.") subfile = tmp_filename() f = open(subfile, 'w') f.write(sub) f.close() else: subfile = '' if verbose: print("Popen called with arguments", end="") print([program, datafile, subfile] + extra_opts) print("") print("CHomP output:") print("") # output = Popen([program, datafile, subfile, extra_opts], cmd = [program, datafile] if subfile: cmd.append(subfile) if extra_opts: cmd.extend(extra_opts) output = Popen(cmd, stdout=PIPE).communicate()[0] if verbose: print(output) print("End of CHomP output") print("") if generators: gens = open(genfile, 'r').read() if verbose: print("Generators:") print(gens) # # process output # if output.find('ERROR') != -1: raise RuntimeError('error inside CHomP') # output contains substrings of one of the forms # "H_1 = Z", "H_1 = Z_2 + Z", "H_1 = Z_2 + Z^2", # "H_1 = Z + Z_2 + Z" if output.find('trivial') != -1: if mod_p: return {0: VectorSpace(base_ring, 0)} else: return {0: HomologyGroup(0, ZZ)} d = {} h = re.compile("^H_([0-9]*) = (.*)$", re.M) tors = re.compile("Z_([0-9]*)") # # homology groups # for m in h.finditer(output): if verbose: print(m.groups()) # dim is the dimension of the homology group dim = int(m.group(1)) # hom_str is the right side of the equation "H_n = Z^r + Z_k + ..." hom_str = m.group(2) # need to read off number of summands and their invariants if hom_str.find("0") == 0: if mod_p: hom = VectorSpace(base_ring, 0) else: hom = HomologyGroup(0, ZZ) else: rk = 0 if hom_str.find("^") != -1: rk_srch = re.search(r'\^([0-9]*)\s?', hom_str) rk = int(rk_srch.group(1)) rk += len(re.findall("(Z$)|(Z\s)", hom_str)) if mod_p: rk = rk if rk != 0 else 1 if verbose: print("dimension = %s, rank of homology = %s" % (dim, rk)) hom = VectorSpace(base_ring, rk) else: n = rk invts = [] for t in tors.finditer(hom_str): n += 1 invts.append(int(t.group(1))) for i in range(rk): invts.append(0) if verbose: print("dimension = %s, number of factors = %s, invariants = %s" % (dim, n, invts)) hom = HomologyGroup(n, ZZ, invts)
# # generators # if generators: if cubical: g = process_generators_cubical(gens, dim) if verbose: print("raw generators: %s" % g) if g: module = CombinatorialFreeModule(base_ring, original_complex.n_cells(dim), prefix="", bracket=True) basis = module.basis() output = [] for x in g: v = module(0) for term in x: v += term[0] * basis[term[1]] output.append(v) g = output elif simplicial: g = process_generators_simplicial(gens, dim, complex) if verbose: print("raw generators: %s" % gens) if g: module = CombinatorialFreeModule(base_ring, complex.n_cells(dim), prefix="", bracket=False) basis = module.basis() output = [] for x in g: v = module(0) for term in x: if complex._is_numeric(): v += term[0] * basis[term[1]] else: translate = complex._translation_from_numeric() simplex = Simplex([translate[a] for a in term[1]]) v += term[0] * basis[simplex] output.append(v) g = output elif chain: g = process_generators_chain(gens, dim, base_ring) if verbose: print("raw generators: %s" % gens) if g: if not mod_p: # sort generators to match up with corresponding invariant g = [_[1] for _ in sorted(zip(invts, g), key=lambda x: x[0])] d[dim] = (hom, g) else: d[dim] = hom else: d[dim] = hom
if chain: new_d = {} diff = complex.differential() if len(diff) == 0: return {} bottom = min(diff) top = max(diff) for dim in d: if complex._degree_of_differential == -1: # chain complex new_dim = bottom + dim else: # cochain complex new_dim = top - dim if isinstance(d[dim], tuple): # generators included. group = d[dim][0] gens = d[dim][1] new_gens = [] dimension = complex.differential(new_dim).ncols() # make sure that each vector is embedded in the # correct ambient space: pad with a zero if # necessary. for v in gens: v_dict = v.dict() if dimension - 1 not in v.dict(): v_dict[dimension - 1] = 0 new_gens.append(vector(base_ring, v_dict)) else: new_gens.append(v) new_d[new_dim] = (group, new_gens) else: new_d[new_dim] = d[dim] d = new_d return d
""" Print a help message for ``program``, a program from the CHomP suite.
:param program: which CHomP program to use :type program: string :return: nothing -- just print a message
EXAMPLES::
sage: from sage.interfaces.chomp import CHomP sage: CHomP().help('homcubes') # optional - CHomP HOMCUBES, ver. ... Copyright (C) ... by Pawel Pilarczyk... """ from subprocess import Popen, PIPE print(Popen([program, '-h'], stdout=PIPE).communicate()[0])
r""" Compute the homology of a simplicial complex using the CHomP program ``homsimpl``. If the argument ``subcomplex`` is present, compute homology of ``complex`` relative to ``subcomplex``.
:param complex: a simplicial complex :param subcomplex: a subcomplex of ``complex`` or None (the default) :param base_ring: ring over which to perform computations -- must be `\ZZ` or `\GF{p}`. :type base_ring: ring; optional, default `\ZZ` :param generators: if True, also return list of generators :type generators: boolean; optional, default False :param verbose: if True, print helpful messages as the computation progresses :type verbose: boolean; optional, default False :param help: if True, just print a help message and exit :type help: boolean; optional, default False :param extra_opts: options passed directly to ``program`` :type extra_opts: string :return: homology groups as a dictionary indexed by dimension
EXAMPLES::
sage: from sage.interfaces.chomp import homsimpl sage: T = simplicial_complexes.Torus() sage: M8 = simplicial_complexes.MooreSpace(8) sage: M4 = simplicial_complexes.MooreSpace(4) sage: X = T.disjoint_union(T).disjoint_union(T).disjoint_union(M8).disjoint_union(M4) sage: homsimpl(X)[1] # optional - CHomP Z^6 x C4 x C8
Relative homology::
sage: S = simplicial_complexes.Simplex(3) sage: bdry = S.n_skeleton(2) sage: homsimpl(S, bdry)[3] # optional - CHomP Z
Generators: these are given as a list after the homology group. Each generator is specified as a linear combination of simplices::
sage: homsimpl(S, bdry, generators=True)[3] # optional - CHomP (Z, [(0, 1, 2, 3)])
sage: homsimpl(simplicial_complexes.Sphere(1), generators=True) # optional - CHomP {0: 0, 1: (Z, [(0, 1) - (0, 2) + (1, 2)])}
TESTS:
Generators for a simplicial complex whose vertices are not integers::
sage: S1 = simplicial_complexes.Sphere(1) sage: homsimpl(S1.join(S1), generators=True, base_ring=GF(2))[3][1] # optional - CHomP [('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')] """ from sage.homology.all import SimplicialComplex help = kwds.get('help', False) if help: return CHomP().help('homsimpl') # Check types here, because if you pass homsimpl a cubical # complex, it tries to compute its homology as if it were a # simplicial complex and gets terribly wrong answers. if (isinstance(complex, SimplicialComplex) and (subcomplex is None or isinstance(subcomplex, SimplicialComplex))): return CHomP()('homsimpl', complex, subcomplex=subcomplex, **kwds) else: raise TypeError("Complex and/or subcomplex are not simplicial complexes.")
r""" Compute the homology of a cubical complex using the CHomP program ``homcubes``. If the argument ``subcomplex`` is present, compute homology of ``complex`` relative to ``subcomplex``.
:param complex: a cubical complex :param subcomplex: a subcomplex of ``complex`` or None (the default) :param base_ring: ring over which to perform computations -- must be `\ZZ` or `\GF{p}`. :type base_ring: ring; optional, default `\ZZ` :param generators: if True, also return list of generators :type generators: boolean; optional, default False :param verbose: if True, print helpful messages as the computation progresses :type verbose: boolean; optional, default False :param help: if True, just print a help message and exit :type help: boolean; optional, default False :param extra_opts: options passed directly to ``homcubes`` :type extra_opts: string :return: homology groups as a dictionary indexed by dimension
EXAMPLES::
sage: from sage.interfaces.chomp import homcubes sage: S = cubical_complexes.Sphere(3) sage: homcubes(S)[3] # optional - CHomP Z
Relative homology::
sage: C3 = cubical_complexes.Cube(3) sage: bdry = C3.n_skeleton(2) sage: homcubes(C3, bdry) # optional - CHomP {0: 0, 1: 0, 2: 0, 3: Z}
Generators: these are given as a list after the homology group. Each generator is specified as a linear combination of cubes::
sage: homcubes(cubical_complexes.Sphere(1), generators=True, base_ring=GF(2))[1][1] # optional - CHomP [[[1,1] x [0,1]] + [[0,1] x [1,1]] + [[0,1] x [0,0]] + [[0,0] x [0,1]]] """ from sage.homology.all import CubicalComplex help = kwds.get('help', False) if help: return CHomP().help('homcubes') # Type-checking is here for the same reason as in homsimpl above. if (isinstance(complex, CubicalComplex) and (subcomplex is None or isinstance(subcomplex, CubicalComplex))): return CHomP()('homcubes', complex, subcomplex=subcomplex, **kwds) else: raise TypeError("Complex and/or subcomplex are not cubical complexes.")
r""" Compute the homology of a chain complex using the CHomP program ``homchain``.
:param complex: a chain complex :param generators: if True, also return list of generators :type generators: boolean; optional, default False :param verbose: if True, print helpful messages as the computation progresses :type verbose: boolean; optional, default False :param help: if True, just print a help message and exit :type help: boolean; optional, default False :param extra_opts: options passed directly to ``homchain`` :type extra_opts: string :return: homology groups as a dictionary indexed by dimension
EXAMPLES::
sage: from sage.interfaces.chomp import homchain sage: C = cubical_complexes.Sphere(3).chain_complex() sage: homchain(C)[3] # optional - CHomP Z
Generators: these are given as a list after the homology group. Each generator is specified as a cycle, an element in the appropriate free module over the base ring::
sage: C2 = delta_complexes.Sphere(2).chain_complex() sage: homchain(C2, generators=True)[2] # optional - CHomP (Z, [(1, -1)]) sage: homchain(C2, generators=True, base_ring=GF(2))[2] # optional - CHomP (Vector space of dimension 1 over Finite Field of size 2, [(1, 1)])
TESTS:
Chain complexes concentrated in negative dimensions, cochain complexes, etc.::
sage: C = ChainComplex({-5: 4 * identity_matrix(ZZ, 2)}, degree=-1) sage: homchain(C) # optional - CHomP {-6: C4 x C4} sage: C = ChainComplex({-5: 4 * identity_matrix(ZZ, 2)}, degree=1) sage: homchain(C, generators=True) # optional - CHomP {-4: (C4 x C4, [(1, 0), (0, 1)])} """ from sage.homology.chain_complex import ChainComplex_class help = kwds.get('help', False) if help: return CHomP().help('homchain') # Type-checking just in case. if isinstance(complex, ChainComplex_class): return CHomP()('homchain', complex, **kwds) else: raise TypeError("Complex is not a chain complex.")
r""" Process CHomP generator information for cubical complexes.
:param gen_string: generator output from CHomP :type gen_string: string :param dim: dimension in which to find generators :type dim: integer :return: list of generators in each dimension, as described below
``gen_string`` has the form ::
The 2 generators of H_1 follow: generator 1 -1 * [(0,0,0,0,0)(0,0,0,0,1)] 1 * [(0,0,0,0,0)(0,0,1,0,0)] ... generator 2 -1 * [(0,1,0,1,1)(1,1,0,1,1)] -1 * [(0,1,0,0,1)(0,1,0,1,1)] ...
Each line consists of a coefficient multiplied by a cube; the cube is specified by its "bottom left" and "upper right" corners.
For technical reasons, we remove the first coordinate of each tuple, and using regular expressions, the remaining parts get converted from a string to a pair ``(coefficient, Cube)``, with the cube represented as a product of tuples. For example, the first line in "generator 1" gets turned into ::
(-1, [0,0] x [0,0] x [0,0] x [0,1])
representing an element in the free abelian group with basis given by cubes. Each generator is a list of such pairs, representing the sum of such elements. These are reassembled in :meth:`CHomP.__call__` to actual elements in the free module generated by the cubes of the cubical complex in the appropriate dimension.
Therefore the return value is a list of lists of pairs, one list of pairs for each generator.
EXAMPLES::
sage: from sage.interfaces.chomp import process_generators_cubical 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)]" sage: process_generators_cubical(s, 1) [[(-1, [0,0] x [0,0] x [0,0] x [0,1]), (1, [0,0] x [0,1] x [0,0] x [0,0])]] sage: len(process_generators_cubical(s, 1)) # only one generator 1 """ # each dim in gen_string starts with "The generator for # H_3 follows:". So search for "T" to find the # end of the current list of generators. # project g to one end of the cylinder [0,1] x complex: # # drop the first coordinate and eliminate duplicates, at least # in positive dimensions, drop any line containing a # degenerate cube else: # line like "generator 2" else: else: coeff = sign * 1 else: output.append(cubes_list) cubes_list = [] else: return None
r""" Process CHomP generator information for simplicial complexes.
:param gen_string: generator output from CHomP :type gen_string: string :param dim: dimension in which to find generators :type dim: integer :param complex: simplicial complex under consideration :return: list of generators in each dimension, as described below
``gen_string`` has the form ::
The 2 generators of H_1 follow: generator 1 -1 * (1,6) 1 * (1,4) ... generator 2 -1 * (1,6) 1 * (1,4) ...
where each line contains a coefficient and a simplex. Each line is converted, using regular expressions, from a string to a pair ``(coefficient, Simplex)``, like ::
(-1, (1,6))
representing an element in the free abelian group with basis given by simplices. Each generator is a list of such pairs, representing the sum of such elements. These are reassembled in :meth:`CHomP.__call__` to actual elements in the free module generated by the simplices of the simplicial complex in the appropriate dimension.
Therefore the return value is a list of lists of pairs, one list of pairs for each generator.
EXAMPLES::
sage: from sage.interfaces.chomp import process_generators_simplicial sage: s = "The 2 generators of H_1 follow:\ngenerator 1:\n-1 * (1,6)\n1 * (1,4)" sage: process_generators_simplicial(s, 1, simplicial_complexes.Torus()) [[(-1, (1, 6)), (1, (1, 4))]] """ # each dim in gen_string starts with "The generator for H_3 # follows:". So search for "T" to find the end of the current # list of generators. else: else: coeff = sign * 1 else: output.append(simplex_list) simplex_list = [] else: return None
r""" Process CHomP generator information for simplicial complexes.
:param gen_string: generator output from CHomP :type gen_string: string :param dim: dimension in which to find generators :type dim: integer :param base_ring: base ring over which to do the computations :type base_ring: optional, default ZZ :return: list of generators in each dimension, as described below
``gen_string`` has the form ::
[H_0] a1
[H_1] a2 a3
[H_2] a1 - a2
For each homology group, each line lists a homology generator as a linear combination of generators ``ai`` of the group of chains in the appropriate dimension. The elements ``ai`` are indexed starting with `i=1`. Each generator is converted, using regular expressions, from a string to a vector (an element in the free module over ``base_ring``), with ``ai`` representing the unit vector in coordinate `i-1`. For example, the string ``a1 - a2`` gets converted to the vector ``(1, -1)``.
Therefore the return value is a list of vectors.
EXAMPLES::
sage: from sage.interfaces.chomp import process_generators_chain sage: s = "[H_0]\na1\n\n[H_1]\na2\na3\n" sage: process_generators_chain(s, 1) [(0, 1), (0, 0, 1)] sage: s = "[H_0]\na1\n\n[H_1]\n5 * a2 - a1\na3\n" sage: process_generators_chain(s, 1, base_ring=ZZ) [(-1, 5), (0, 0, 1)] sage: process_generators_chain(s, 1, base_ring=GF(2)) [(1, 1), (0, 0, 1)] """ # each dim in gens starts with a string like # "[H_3]". So search for "[" to find the end of # the current list of generators. # each line in the string g is a linear # combination of things like "a2", "a31", etc. # indexing on the a's starts at 1. else: else: |