Ticket #12553: trac_12553_palp_database.patch

File trac_12553_palp_database.patch, 16.5 KB (added by vbraun, 10 years ago)

Updated patch

  • new file sage/geometry/polyhedron/palp_database.py

    # HG changeset patch
    # User Volker Braun <vbraun@stp.dias.ie>
    # Date 1341159870 -3600
    # Node ID 5947def5642fa35471081616b4d812f9464a035e
    # Parent  6c2ab10b686cb990b61ff95c152bc5a696f00f4c
    Trac #12553: Add interface for PALP polytope databases
    
    Implements a new PALPreader class to parse the PALP listing of polytopes
    
    diff --git a/sage/geometry/polyhedron/palp_database.py b/sage/geometry/polyhedron/palp_database.py
    new file mode 100644
    - +  
     1"""
     2Access the PALP database(s) of reflexive lattice polytopes
     3
     4EXAMPLES::
     5
     6    sage: from sage.geometry.polyhedron.palp_database import PALPreader
     7    sage: for lp in PALPreader(2):
     8    ...       cone = Cone([(1,r[0],r[1]) for r in lp.vertices()])
     9    ...       fan = Fan([cone])
     10    ...       X = ToricVariety(fan)
     11    ...       ideal = X.affine_algebraic_patch(cone).defining_ideal()
     12    ...       print lp.n_vertices(), ideal.hilbert_series()
     13    3 (-t^2 - 7*t - 1)/(t^3 - 3*t^2 + 3*t - 1)
     14    3 (-t^2 - t - 1)/(t^3 - 3*t^2 + 3*t - 1)
     15    3 (t^2 + 6*t + 1)/(-t^3 + 3*t^2 - 3*t + 1)
     16    3 (t^2 + 2*t + 1)/(-t^3 + 3*t^2 - 3*t + 1)
     17    3 (t^2 + 4*t + 1)/(-t^3 + 3*t^2 - 3*t + 1)
     18    4 (-t^2 - 5*t - 1)/(t^3 - 3*t^2 + 3*t - 1)
     19    4 (-t^2 - 3*t - 1)/(t^3 - 3*t^2 + 3*t - 1)
     20    4 (t^2 + 2*t + 1)/(-t^3 + 3*t^2 - 3*t + 1)
     21    4 (t^2 + 6*t + 1)/(-t^3 + 3*t^2 - 3*t + 1)
     22    4 (t^2 + 6*t + 1)/(-t^3 + 3*t^2 - 3*t + 1)
     23    4 (t^2 + 2*t + 1)/(-t^3 + 3*t^2 - 3*t + 1)
     24    4 (t^2 + 4*t + 1)/(-t^3 + 3*t^2 - 3*t + 1)
     25    5 (-t^2 - 3*t - 1)/(t^3 - 3*t^2 + 3*t - 1)
     26    5 (-t^2 - 5*t - 1)/(t^3 - 3*t^2 + 3*t - 1)
     27    5 (t^2 + 4*t + 1)/(-t^3 + 3*t^2 - 3*t + 1)
     28    6 (t^2 + 4*t + 1)/(-t^3 + 3*t^2 - 3*t + 1)
     29"""
     30
     31from subprocess import Popen, PIPE
     32
     33from sage.structure.sage_object import SageObject
     34from sage.matrix.all import matrix
     35from sage.rings.all import Integer, ZZ
     36
     37from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
     38from sage.geometry.polyhedron.constructor import Polyhedron
     39
     40
     41
     42#########################################################################
     43class PALPreader(SageObject):
     44    """
     45    Read PALP database of polytopes.
     46
     47
     48    INPUT:
     49   
     50    - ``dim`` -- integer. The dimension of the poylhedra
     51       
     52    - ``data_basename`` -- string or ``None`` (default). The directory
     53      and database base filename (PALP usually uses ``'zzdb'``) name
     54      containing the PALP database to read. Defaults to the built-in
     55      database location.
     56       
     57    - ``output`` -- string. How to return the reflexive polyhedron
     58      data. Allowed values = ``'list'``, ``'Polyhedron'`` (default),
     59      ``'pointcollection'``, and ``'PPL'``. Case is ignored.
     60
     61    EXAMPLES::
     62
     63        sage: from sage.geometry.polyhedron.palp_database import PALPreader
     64        sage: polygons = PALPreader(2)
     65        sage: [ (p.n_Vrepresentation(), len(p.integral_points())) for p in polygons ]
     66        [(3, 4), (3, 10), (3, 5), (3, 9), (3, 7), (4, 6), (4, 8), (4, 9),
     67         (4, 5), (4, 5), (4, 9), (4, 7), (5, 8), (5, 6), (5, 7), (6, 7)]
     68
     69        sage: iter(PALPreader(2, output='list')).next()
     70        [[1, 0], [0, 1], [-1, -1]]
     71        sage: type(_)
     72        <type 'list'>
     73
     74        sage: iter(PALPreader(2, output='Polyhedron')).next()
     75        A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 3 vertices
     76        sage: type(_)
     77        <class 'sage.geometry.polyhedron.backend_ppl.Polyhedra_ZZ_ppl_with_category.element_class'>
     78
     79        sage: iter(PALPreader(2, output='PPL')).next()
     80        A 2-dimensional lattice polytope in ZZ^2 with 3 vertices
     81        sage: type(_)
     82        <class 'sage.geometry.polyhedron.ppl_lattice_polytope.LatticePolytope_PPL_class'>
     83
     84        sage: iter(PALPreader(2, output='PointCollection')).next()
     85        [ 1,  0],
     86        [ 0,  1],
     87        [-1, -1]
     88        in Ambient free module of rank 2 over the principal ideal domain Integer Ring
     89        sage: type(_)
     90        <type 'sage.geometry.point_collection.PointCollection'>
     91    """
     92
     93    def __init__(self, dim, data_basename=None, output='Polyhedron'):
     94        """
     95        The Python constructor
     96       
     97        See :class:`PALPreader` for documentation.
     98       
     99        TESTS::
     100       
     101            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     102            sage: polygons = PALPreader(2)
     103        """
     104        self._dim = dim
     105        if data_basename is not None:
     106            self._data_basename = data_basename
     107        else:
     108            import os
     109            SAGE_DATA = os.getenv('SAGE_DATA')
     110            self._data_basename = os.path.join(SAGE_DATA, 'reflexive_polytopes',
     111                                               'Full'+str(dim)+'d', 'zzdb')
     112            info = self._data_basename + '.info'
     113            if not os.path.exists(info):
     114                raise ValueError('Cannot find PALP database: '+info)
     115        from sage.geometry.polyhedron.parent import Polyhedra
     116        self._polyhedron_parent = Polyhedra(ZZ, dim)
     117        self._output = output.lower()
     118       
     119    def _palp_Popen(self):
     120        """
     121        Open PALP.
     122       
     123        OUTPUT:
     124
     125        A PALP subprocess.
     126
     127        EXAMPLES::
     128
     129            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     130            sage: polygons = PALPreader(2)
     131            sage: polygons._palp_Popen()
     132            <subprocess.Popen object at 0x...>
     133        """
     134        return Popen(["class.x", "-b2a", "-di", self._data_basename], stdout=PIPE)
     135
     136    def _read_vertices(self, stdout, rows, cols):
     137        """
     138        Read vertex data from the PALP output pipe.
     139
     140        OUTPUT:
     141       
     142        A list of lists.
     143
     144        EXAMPLES::
     145
     146            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     147            sage: polygons = PALPreader(2)
     148            sage: palp = polygons._palp_Popen()
     149            sage: palp.stdout.readline()
     150            '2 3  \n'
     151            sage: polygons._read_vertices(palp.stdout, 2, 3)
     152            [[1, 0], [0, 1], [-1, -1]]
     153        """
     154        m = [ [] for col in range(0,cols) ]
     155        for row in range(0,rows):
     156            for col,x in enumerate(stdout.readline().split()):
     157                m[col].append(ZZ(x))
     158        return m
     159   
     160    def _read_vertices_transposed(self, stdout, rows, cols):
     161        """
     162        Read vertex data from the PALP output pipe.
     163
     164        OUTPUT:
     165       
     166        A list of lists.
     167
     168        EXAMPLES::
     169
     170            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     171            sage: polygons = PALPreader(2)
     172            sage: palp = polygons._palp_Popen()
     173            sage: palp.stdout.readline()
     174            '2 3  \n'
     175            sage: polygons._read_vertices_transposed(palp.stdout, 2, 3)
     176            [[1, 0, -1], [0, 1, -1]]
     177        """
     178        m = []
     179        for row in range(0,rows):
     180            m.append( [ ZZ(x) for x in stdout.readline().split() ] )
     181        return m
     182   
     183    def _iterate_list(self, start, stop, step):
     184        """
     185        Iterate over the reflexive polytopes.
     186
     187        INPUT:
     188
     189        - ``start``, ``stop``, ``step`` -- integers specifying the
     190          range to iterate over.
     191
     192        OUTPUT:
     193
     194        A generator for vertex data as a list of lists.
     195
     196        EXAMPLES::
     197
     198            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     199            sage: polygons = PALPreader(2)
     200            sage: iter = polygons._iterate_list(0,4,2)
     201            sage: iter.next()
     202            [[1, 0], [0, 1], [-1, -1]]
     203        """
     204        if start is None:
     205            start = 0
     206        if step is None:
     207            step = 1
     208        palp = self._palp_Popen()
     209        try:
     210            palp_out = palp.stdout
     211            i = 0
     212            while True:
     213                l = palp_out.readline().strip()
     214                if l=='' or l.startswith('#'):
     215                    return  # EOF
     216                l=l.split()
     217                dim = ZZ(l[0]);  # dimension
     218                n = ZZ(l[1]);    # number of vertices
     219                if i>=start and (i-start) % step == 0:
     220                    if dim == self._dim:
     221                        vertices = self._read_vertices(palp_out, dim, n)
     222                    elif n == self._dim:  # PALP sometimes returns transposed data #@!#@
     223                        vertices = self._read_vertices_transposed(palp_out, dim, n)
     224                    else:
     225                        raise ValueError('PALP output dimension mismatch.')
     226                    yield vertices
     227                else:
     228                    for row in range(0,dim):
     229                        palp_out.readline()
     230                i += 1
     231                if stop is not None and i>=stop:
     232                    return
     233        finally:
     234            palp.poll()
     235            if palp.returncode is None:
     236                palp.terminate()
     237            palp.poll()
     238            if palp.returncode is None:
     239                palp.kill()
     240   
     241
     242    def _iterate_Polyhedron(self, start, stop, step):
     243        """
     244        Iterate over the reflexive polytopes.
     245
     246        INPUT:
     247
     248        - ``start``, ``stop``, ``step`` -- integers specifying the
     249          range to iterate over.
     250
     251        OUTPUT:
     252
     253        A generator for lattice polyhedra.
     254
     255        EXAMPLES::
     256
     257            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     258            sage: polygons = PALPreader(2)
     259            sage: iter = polygons._iterate_Polyhedron(0,4,2)
     260            sage: iter.next()
     261            A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 3 vertices
     262        """
     263        parent = self._polyhedron_parent
     264        for vertices in self._iterate_list(start, stop, step):
     265            p = parent.element_class(parent, [vertices,[],[]], None)
     266            yield p
     267
     268    def _iterate_PPL(self, start, stop, step):
     269        """
     270        Iterate over the reflexive polytopes.
     271
     272        INPUT:
     273
     274        - ``start``, ``stop``, ``step`` -- integers specifying the
     275          range to iterate over.
     276
     277        OUTPUT:
     278
     279        A generator for PPL-based lattice polyhedra.
     280
     281        EXAMPLES::
     282
     283            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     284            sage: polygons = PALPreader(2)
     285            sage: iter = polygons._iterate_PPL(0,4,2)
     286            sage: iter.next()
     287            A 2-dimensional lattice polytope in ZZ^2 with 3 vertices
     288        """
     289        for vertices in self._iterate_list(start, stop, step):
     290            yield LatticePolytope_PPL(*vertices)
     291
     292    def _iterate_PointCollection(self, start, stop, step):
     293        """
     294        Iterate over the reflexive polytopes.
     295
     296        INPUT:
     297
     298        - ``start``, ``stop``, ``step`` -- integers specifying the
     299          range to iterate over.
     300
     301        OUTPUT:
     302
     303        A generator for PPL-based lattice polyhedra.
     304
     305        EXAMPLES::
     306
     307            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     308            sage: polygons = PALPreader(2)
     309            sage: iter = polygons._iterate_PointCollection(0,4,2)
     310            sage: iter.next()
     311            [ 1,  0],
     312            [ 0,  1],
     313            [-1, -1]
     314            in Ambient free module of rank 2 over the principal ideal domain Integer Ring
     315        """
     316        from sage.modules.free_module import FreeModule
     317        N = FreeModule(ZZ, self._dim)
     318        from sage.geometry.point_collection import PointCollection
     319        for vertices in self._iterate_list(start, stop, step):
     320            yield PointCollection(vertices, module=N)
     321
     322    def _iterate(self, output=None):
     323        """
     324        Iterate over the reflexive polytopes.
     325
     326        INPUT:
     327
     328        - ``output`` -- as in the :class:`PALPreader` constructor.
     329
     330        OUTPUT:
     331
     332        A function generating lattice polytopes in the specified output format.
     333
     334        EAMPLES::
     335
     336            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     337            sage: polygons = PALPreader(2)
     338            sage: func = polygons._iterate(output='list')
     339            sage: func
     340            <bound method PALPreader._iterate_list of <class 'sage.geometry.polyhedron.palp_database.PALPreader'>>
     341            sage: iter = func(0,1,1)
     342            sage: iter.next()
     343            [[1, 0], [0, 1], [-1, -1]]
     344        """
     345        if output is None:
     346            output = self._output
     347        if output == 'polyhedron':
     348            return self._iterate_Polyhedron
     349        elif output == 'ppl':
     350            return self._iterate_PPL
     351        elif output == 'pointcollection':
     352            return self._iterate_PointCollection
     353        elif output == 'list':
     354            return self._iterate_list
     355        else:
     356            raise TypeError('Unknown output format (='+str(self._output)+').')
     357
     358    def __iter__(self):
     359        """
     360        Iterate over all polytopes.
     361       
     362        OUTPUT:
     363
     364        An iterator for all polytopes.
     365
     366        TESTS::
     367       
     368            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     369            sage: polygons = PALPreader(2)
     370            sage: polygons.__iter__()
     371            <generator object _iterate_Polyhedron at 0x...>
     372        """
     373        iterator = self._iterate()
     374        return iterator(None, None, None)
     375
     376    def __getitem__(self, item):
     377        """
     378        Return the polytopes(s) indexed by ``item``.
     379
     380        EXAMPLES::
     381
     382            sage: from sage.geometry.polyhedron.palp_database import PALPreader
     383            sage: palp = PALPreader(3)
     384            sage: list(palp[10:30:10])
     385            [A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
     386             A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 4 vertices]
     387        """
     388        iterator = self._iterate()
     389        if isinstance(item, slice):
     390            return iterator(item.start, item.stop, item.step)
     391        else:
     392            try:
     393                return iterator(item, item+1, 1).next()
     394            except StopIteration:
     395                raise IndexError('Index out of range.')
     396
     397
     398
     399#########################################################################
     400class Reflexive4dHodge(PALPreader):
     401    """
     402    Read the PALP database for Hodge numbers of 4d polytopes.
     403
     404    The database is very large and not installed by default. You can
     405    install it with the command ``install_package('polytopes_db_4d')``.
     406   
     407    INPUT:
     408
     409    - ``h11``, ``h21`` -- Integers. The Hodge numbers of the reflexive
     410      polytopes to list.
     411
     412    Any additional keyword arguments are passed to
     413    :class:`PALPreader`.
     414
     415    EXAMPLES:
     416
     417        sage: from sage.geometry.polyhedron.palp_database import Reflexive4dHodge
     418        sage: ref = Reflexive4dHodge(1,101)             # optional - polytopes_db_4d
     419        sage: iter(ref).next().Vrepresentation()        # optional - polytopes_db_4d
     420        (A vertex at (-1, -1, -1, -1), A vertex at (0, 0, 0, 1),
     421         A vertex at (0, 0, 1, 0), A vertex at (0, 1, 0, 0), A vertex at (1, 0, 0, 0))
     422    """
     423    def __init__(self, h11, h21, data_basename=None, **kwds):
     424        """
     425        The Python constructor.
     426       
     427        See :class:`Reflexive4dHodge` for documentation.
     428
     429        TESTS::
     430
     431        sage: from sage.geometry.polyhedron.palp_database import Reflexive4dHodge
     432        sage: Reflexive4dHodge(1,101)  # optional - polytopes_db_4d
     433        <class 'sage.geometry.polyhedron.palp_database.Reflexive4dHodge'>
     434        """
     435        dim = 4
     436        if data_basename is None:
     437            import os
     438            SAGE_DATA = os.getenv('SAGE_DATA')
     439            data_basename = os.path.join(SAGE_DATA, 'reflexive_polytopes',
     440                                         'Hodge4d', 'all')
     441            info = data_basename + '.vinfo'
     442            if not os.path.exists(info):
     443                raise ValueError('Cannot find PALP database: '+info+
     444                                 '. Did you install the polytopes_db_4d optional spkg?')
     445        PALPreader.__init__(self, dim, data_basename=data_basename, **kwds)
     446        self._h11 = h11
     447        self._h21 = h21
     448       
     449    def _palp_Popen(self):
     450        """
     451        Open PALP.
     452       
     453        OUTPUT:
     454
     455        A PALP subprocess.
     456
     457        EXAMPLES::
     458
     459            sage: from sage.geometry.polyhedron.palp_database import Reflexive4dHodge
     460            sage: polygons = Reflexive4dHodge(1, 101)   # optional - polytopes_db_4d
     461            sage: polygons._palp_Popen()                # optional - polytopes_db_4d
     462            <subprocess.Popen object at 0x...>
     463        """
     464        return Popen(['class-4d.x', '-He',
     465                      'H'+str(self._h21)+':'+str(self._h11)+'L100000000',
     466                      '-di', self._data_basename], stdout=PIPE)
     467       
     468
     469