Ticket #14116: trac_14116_polymake_upgrade.patch

File trac_14116_polymake_upgrade.patch, 41.0 KB (added by tkluck, 9 years ago)
  • module_list.py

    # HG changeset patch
    # User Timo Kluck <tkluck@infty.nl>
    # Date 1356283340 0
    # Node ID fd215fcf41cc1b6f27057f066c989fc6c6124e0e
    # Parent  db92b42520772efe29a41110c8aa0f850950f519
    Add wrappers for new version of polymake
    
    diff --git a/module_list.py b/module_list.py
    a b  
    19641964                  depends = [SAGE_LOCAL + "/include/lrcalc/symfcn.h"]), # should include all .h
    19651965        )
    19661966
     1967if is_package_installed('polymake'):
     1968    ext_modules.append(
     1969        Extension('sage.geometry.polymake',
     1970                  sources = ["sage/geometry/polymake/polymake.pyx"],
     1971                  include_dirs = [SAGE_LOCAL + '/include'],
     1972                  libraries=["polymake","gmp","xml2"],
     1973                  library_dirs=[SAGE_LOCAL + "/lib"],
     1974                  extra_compile_args=["-DPOLYMAKE_DEBUG=0"],
     1975                  extra_link_args=["-Wl,-rpath,%s/lib" % SAGE_LOCAL],
     1976                  language="c++",
     1977                  depends = ["sage/geometry/polymake/polymake_wrappers.h"])
     1978        )
  • new file sage/geometry/polymake/polymake.pyx

    diff --git a/sage/geometry/polymake/polymake.pyx b/sage/geometry/polymake/polymake.pyx
    new file mode 100644
    - +  
     1"""
     2This module provides access to polymake, which 'has been developed
     3since 1997 in the Discrete Geometry group at the Institute of
     4Mathematics of Technische Universitat Berlin. Since 2004 the
     5development is shared with Fachbereich Mathematik, Technische
     6Universitat Darmstadt. The system offers access to a wide variety
     7of algorithms and packages within a common framework. polymake is
     8flexible and continuously expanding. The software supplies C++ and
     9Perl interfaces which make it highly adaptable to individual
     10needs.'
     11
     12AUTHORS:
     13
     14- Ewgenij Gawrilow, Michael Joswig: main authors of polymake
     15
     16- Timo Kluck (2012-11-29): total rewrite to take advantage of libary interface,
     17  keeping backwards compatibility
     18
     19- William Stein: Sage interface
     20
     21EXAMPLES::
     22
     23Use :func:`polytope` to construct cones. The easiest way it to define it by a
     24set of affine equations::
     25
     26    sage: import sage.geometry.polymake as pm
     27    sage: p1 = pm.polytope([x < 3, x > 0], homogeneous_coordinates=(1,x))
     28    sage: p2 = pm.polytope([x < 2, x > -3], homogeneous_coordinates=(1,x))
     29    sage: p3 = p1.intersection(p2)
     30    sage: p3.get_defining_equations()
     31    [x > 0, -x + 2 > 0]
     32    sage: p3.volume()
     33    2
     34
     35Another option is to use a set of homogeneous equations::
     36
     37    sage: x,y,z,w = var('x,y,z,w')
     38    sage: p1 = pm.polytope([x > 0, y > 0, z > 0], homogeneous_coordinates=(x,y,z))
     39    sage: p1.get_defining_equations()
     40    [x > 0, y > 0, z > 0]
     41    sage: p1.is_simple()
     42    True
     43    sage: p1.is_feasible()
     44    True
     45    sage: p1.volume()
     46    1
     47
     48This polytope is empty because of the extra constraint::
     49
     50    sage: p1 = pm.polytope([x+y+z==1, x> 1/2, y > 1/2, z > 1/2], homogeneous_coordinates = (1,x,y,z))
     51    sage: #p1.get_defining_equations()
     52    [x + y + z - 1 == 0, x - 1/10 > 0, y - 1/10 > 0, z - 1/10 > 0]
     53    sage: p1.is_feasible()
     54    False
     55    sage: p1.volume()
     56    0
     57
     58Another call allows you to specify affine equations, which are then interpreted
     59on the affine space you spacify::
     60
     61    sage: p1 = pm.polytope([x > 1/10, y > 1/10, z > 1/10], homogeneous_coordinates=(x,y,z), affine_plane_equation=(x+y+z==1))
     62
     63
     64So these two calls are equivalent::
     65
     66    sage: p1 = pm.polytope([x > 1/2, y > 1/2, z > 1/2], homogeneous_coordinates = (1,x,y,z))
     67    sage: p1 = pm.polytope([x > 1/2, y > 1/2, z > 1/2], homogeneous_coordinates = (w,x,y,z), affine_plane_equation=(w==1))
     68
     69You can obtain any polymake property, even if it hasn't been wrapped, by
     70using the underlying Perl object::
     71
     72    sage: p1._polymake_obj.VERTICES_IN_FACETS
     73    '{0 1 3}\n{0 1 2}\n{0 2 3}\n{1 2 3}\n'
     74
     75
     76Doctests from the previous version of this interface::
     77
     78    sage: P = polymake.convex_hull([[1,0,0,0], [1,0,0,1], [1,0,1,0], [1,0,1,1],  [1,1,0,0], [1,1,0,1], [1,1,1,0], [1,1,1,1]])   # optional - polymake
     79    sage: P.facets()                            # optional - polymake
     80    [(0, 0, 0, 1), (0, 1, 0, 0), (0, 0, 1, 0), (1, 0, 0, -1), (1, 0, -1, 0), (1, -1, 0, 0)]
     81    sage: P.vertices()                            # optional - polymake
     82    [(1, 0, 0, 0), (1, 0, 0, 1), (1, 0, 1, 0), (1, 0, 1, 1), (1, 1, 0, 0), (1, 1, 0, 1), (1, 1, 1, 0), (1, 1, 1, 1)]
     83    sage: P.is_simple()                              # optional - polymake
     84    True
     85    sage: P.n_facets()                            # optional - polymake
     86    6
     87    sage: polymake.cell24()            # optional - polymake
     88    The 24-cell
     89    sage: R.<x,y,z> = PolynomialRing(QQ,3)
     90    sage: f = x^3 + y^3 + z^3 + x*y*z
     91    sage: e = f.exponents()
     92    sage: a = [[1] + list(v) for v in e]
     93    sage: a
     94    [[1, 3, 0, 0], [1, 0, 3, 0], [1, 1, 1, 1], [1, 0, 0, 3]]
     95    sage: n = polymake.convex_hull(a)       # optional - polymake
     96    sage: n                                 # optional - polymake
     97    Convex hull of points [[1, 0, 0, 3], [1, 0, 3, 0], [1, 1, 1, 1], [1, 3, 0, 0]]
     98    sage: n.facets()                        # optional - polymake
     99    [(0, 1, 0, 0), (3, -1, -1, 0), (0, 0, 1, 0)]
     100    sage: n.is_simple()                     # optional - polymake
     101    True
     102    sage: n.graph()                         # optional - polymake
     103    'GRAPH\n{1 2}\n{0 2}\n{0 1}\n\n'
     104   
     105
     106
     107NOTE:
     108
     109The part that wraps the C++ interface to polymake is engineered to work
     110with Cython 0.17pre. It appears that during the preparation of the actual 0.17
     111release, lots of things have changed (improved) in the Cython C++ implementation,
     112specifically in the stdlib wrappers. I've made a note where I know that a Cython
     113upgrade will need a change in this file. In general, it should be possible
     114to get rid of the explicit Python object -> std::string casts, that I'm now doing
     115by
     116   
     117    std::string cpp_string = <string><char*>python_string
     118
     119because that cast is builtin for higher versions of Cython.
     120
     121CITE:
     122
     123Ewgenij Gawrilow and Michael Joswig. polymake: a framework for analyzing
     124convex polytopes. Polytopes—combinatorics and computation (Oberwolfach,
     1251997), 43–73, DMV Sem., 29, Birkhäuser, Basel, 2000. MR1785292 (2001f:52033).
     126
     127"""
     128#*****************************************************************************
     129#       Copyright (C) 2012 Timo Kluck <tkluck@infty.nl>
     130#       Copyright (C) 2012 William Stein <wstein@gmail.com>
     131#
     132#  Distributed under the terms of the GNU General Public License (GPL)
     133#  as published by the Free Software Foundation; either version 2 of
     134#  the License, or (at your option) any later version.
     135#                  http://www.gnu.org/licenses/
     136#*****************************************************************************
     137
     138from libcpp.string cimport string
     139from libcpp.vector cimport vector
     140import operator
     141from sage.matrix.constructor import matrix
     142from sage.modules.all import vector as sage_vector
     143from sage.rings.rational_field import QQ
     144from sage.symbolic.expression import Expression
     145from sage.symbolic.ring import SR
     146from sage.structure.sage_object import SageObject
     147import tempfile
     148import subprocess
     149
     150
     151cdef extern from "polymake/Matrix.h" namespace "pm":
     152    cdef cppclass Matrix[T]:
     153        Matrix(int, int, vector   #.iterator  # needed in 17.2, not in 17pre
     154                            ) except +
     155        T operator[](int, int) except +
     156
     157cdef extern from "polymake/Rational.h" namespace "pm":
     158    cdef cppclass Rational:
     159        Rational()
     160        Rational(int, int)
     161
     162cdef extern from "polymake/Rational.h" namespace "":
     163    int numerator(Rational)
     164    int denominator(Rational)
     165
     166cdef extern from "polymake/Main.h" namespace "polymake":
     167    cdef cppclass Main:
     168        Main(string)
     169        void set_application(string) except +
     170 
     171cdef extern from "polymake/client.h" namespace "pm::perl":
     172    cdef cppclass PropertyValue:
     173        void operator>>(string) except +
     174
     175    cdef cppclass PropertyOut:
     176        void operator<<(string) except +
     177        void operator<<(Matrix[Rational]) except +
     178
     179    cdef cppclass Object:
     180        Object()
     181        Object(string) except +
     182        #PropertyValue give(string) #except +  #/// FIXME: cython wants to instantiate PropertyValue if we
     183                                              # add exception handling, but it has a private constructor
     184        #PropertyOut take(string)  #except +   #/// idem
     185                                              # the workaround is to use the wrappers below
     186        string description() except +
     187        void save(string) except +
     188
     189cdef extern from "polymake/client.h" namespace "pm::perl::Object":
     190    Object load(string) except +
     191
     192cdef extern from "polymake_wrappers.h":
     193    #These classes wrap the preprocessor macro CallPolymakeFunction. This
     194    #is easier that figuring out all the type casts and template instantiations
     195    #involved in the macro expanded C++ code
     196    void VoidCallPolymakeFunction_wrapper(string) except +
     197    void VoidCallPolymakeFunction_wrapper(string, Object) except +
     198    void VoidCallPolymakeFunction_wrapper(string, Object, Object) except +
     199    void VoidCallPolymakeFunction_wrapper(string, Object, Object, Object) except +
     200    Object *CallPolymakeFunction_wrapper(string) except +
     201    Object *CallPolymakeFunction_wrapper(string, int) except +
     202    Object *CallPolymakeFunction_wrapper(string, int, bint) except +
     203    Object *CallPolymakeFunction_wrapper(string, int, Rational) except +
     204    Object *CallPolymakeFunction_wrapper(string, Object) except +
     205    Object *CallPolymakeFunction_wrapper(string, Object, Object) except +
     206    Object *CallPolymakeFunction_wrapper(string, Object, Object, Object) except +
     207    bint boolCallPolymakeFunction_wrapper(string, Object, Object) except +
     208    Rational rationalCallPolymakeFunction_wrapper(string, Object, Object) except +
     209    # these wrappers cast an expression that may throw an exception.
     210    # Cython cannot compile that expression because it wants to instantiate
     211    # a temporary copy of a class that has a private constructor
     212    void set_property_wrapper(Object, string, Matrix[Rational]) except +
     213    void set_property_wrapper(Object, string, string) except +
     214    string get_property_wrapper(Object, string) except +
     215
     216cdef class _Main:
     217    #The class _Main initializes the polymake interface. It is
     218    #instantiated by this module and has no user-serviceable parts
     219    cdef Main *thisptr
     220    def __cinit__(self, user_settings):
     221        self.thisptr = new Main(<string><char*>user_settings)
     222    def __dealloc__(self):
     223        del self.thisptr
     224    def set_application(self, application):
     225        self.thisptr.set_application(<string><char*>application)
     226
     227cdef class PerlObject:
     228    #This class is a Python wrapper for the Object cython class. The latter
     229    #one is a C++ proxy for Perl objects
     230    cdef Object *thisptr
     231    def __cinit__(self, typename=None):
     232        if typename is not None:
     233            self.thisptr = new Object(<string><char*>typename)
     234        else:
     235            self.thisptr = new Object()
     236    def __dealloc__(self):
     237        del self.thisptr
     238
     239    def load(self, filename):
     240        self.thisptr[0] = load(<string><char*>filename)
     241
     242    def save(self, filename):
     243        self.thisptr.save(<string><char*>filename)
     244
     245    def __getattr__(self, propertyname):
     246        propertyname = propertyname.upper()
     247        return <str>get_property_wrapper(self.thisptr[0], <string><char*>propertyname).c_str()
     248
     249    def __setattr__(self, propertyname, value):
     250        propertyname = propertyname.upper()
     251        cdef int rows
     252        cdef int cols
     253        cdef vector[Rational] entries
     254        cdef Matrix[Rational] *matrix
     255        cdef string val
     256        if isinstance(value, list):
     257            # assume that we've got a matrix
     258            rows = len(value)
     259            cols = len(value[0]) if rows > 0 else 0
     260            flat_list = sum([list(row) for row in value], [])
     261            for i in flat_list:
     262                if hasattr(i, 'numerator') and hasattr(i, 'denominator'):
     263                    p = i.numerator() if callable(i.numerator) else i.numerator
     264                    q = i.denominator() if callable(i.denominator) else i.denominator
     265                    entries.push_back(Rational(p,q))
     266                else:
     267                    entries.push_back(<Rational><int>i)
     268            matrix = new Matrix[Rational](rows, cols, entries.begin())
     269            set_property_wrapper(self.thisptr[0], <string><char*>propertyname, matrix[0])
     270            del matrix
     271        else:
     272            val = <string><char*>value
     273            set_property_wrapper(self.thisptr[0], <string><char*>propertyname, val)
     274
     275
     276def call_function(name, *args, **kwds):
     277    return _call_function(name, args, kwds)
     278
     279cdef _are_instances(values, classes):
     280    return len(values) == len(classes) and all([isinstance(value, klass) for value, klass in zip(values, classes)])
     281
     282cdef _call_function(name, args, kwds):
     283
     284    cdef string name_str = <string><char*>name
     285    cdef int int_arg
     286    cdef bint bool_arg
     287    cdef Rational rational_arg
     288    cdef Rational rational_ret
     289    cdef PerlObject ob1
     290    cdef PerlObject ob2
     291
     292    retval = PerlObject()
     293    if _are_instances(args, (PerlObject, PerlObject)):
     294        ob1 = args[0]
     295        ob2 = args[1]
     296        if kwds.has_key('return_type'):
     297            if kwds['return_type'] == bool:
     298                return boolCallPolymakeFunction_wrapper(name_str, ob1.thisptr[0], ob2.thisptr[0])
     299            elif kwds['return_type'] == QQ:
     300                raise NotImplementedError
     301                #rational_ret = rationalCallPolymakeFunction_wrapper(name_str, ob1.thisptr[0], ob2.thisptr[0])
     302                #return QQ(numerator(rational_ret), denominator(rational_ret))
     303        else:
     304            del retval.thisptr;
     305            retval.thisptr = CallPolymakeFunction_wrapper(name_str, ob1.thisptr[0], ob2.thisptr[0])
     306    elif _are_instances(args, tuple()):
     307        del retval.thisptr;
     308        retval.thisptr = CallPolymakeFunction_wrapper(name_str)
     309    elif _are_instances(args, (int,)):
     310        int_arg = args[0]
     311        del retval.thisptr;
     312        retval.thisptr =  CallPolymakeFunction_wrapper(name_str, int_arg)
     313    elif _are_instances(args, (int,bool)):
     314        int_arg = args[0]
     315        bool_arg = args[1]
     316        del retval.thisptr;
     317        retval.thisptr =  CallPolymakeFunction_wrapper(name_str, int_arg, bool_arg)
     318    elif len(args) == 2 and isinstance(args[0], int) and args[1] in QQ:
     319        int_arg = args[0]
     320        rational_arg = Rational(args[1].numerator(), args[1].denominator())
     321        del retval.thisptr;
     322        retval.thisptr =  CallPolymakeFunction_wrapper(name_str, int_arg, rational_arg)
     323    return retval
     324
     325class Polytope(SageObject):
     326    """
     327    A class representing a polytope
     328    """
     329
     330    def __init__(self, data=None, homogeneous_coordinates=None, affine_plane_equation=None, name=None):
     331        """
     332        return a new polytope object
     333
     334        """
     335        self._homogeneous_coordinates = homogeneous_coordinates
     336        self._affine_plane_equation = affine_plane_equation
     337        self._name = "A Polytope" if name is None else name
     338        if isinstance(data, PerlObject):
     339            self._polymake_obj = data
     340        elif isinstance(data, str):
     341            # backward compatibity, assume that we got a
     342            # (filename, description) pair
     343            self.load(data)
     344            if isinstance(homogeneous_coordinates, str):
     345                self._name = homogeneous_coordinates
     346                self._homogeneous_coordinates = None
     347        else:
     348            self._polymake_obj = PerlObject("Polytope<Rational>")
     349            if data:
     350                if isinstance(data, list):
     351                    if isinstance(data[0], (Expression, bool)):
     352                        # Polytope([x < 2, x > 1], (1,x))
     353                        self.set_defining_equations(data, homogeneous_coordinates, affine_plane_equation)
     354                    else:
     355                        # assume that it is numeric and that they represent inequalities
     356                        self._set_inequalities_coefficients(data)
     357                   
     358
     359    def _repr_(self):
     360        return self._name
     361
     362    def load(self, filename):
     363        self._polymake_obj = PerlObject()
     364        self._polymake_obj.load(filename)
     365
     366    def write(self, filename):
     367        self._polymake_obj.save(filename)
     368
     369    def _parse_matrix(self, string):
     370        row_strings = string.strip().split("\n")
     371        return [[QQ(number_string) for number_string in row_string.strip().split(" ")] for row_string in row_strings]
     372
     373    def _parse_vector(self, string):
     374        numbers_string = string.strip()
     375        return [QQ(number_string) for number_string in numbers_string.split(" ")]
     376
     377    def _parse_bool(self, string):
     378        return '1' in string
     379
     380    def _parse_rational(self, string):
     381        return QQ(string)
     382
     383    def _parse_int(self, string):
     384        return int(string)
     385
     386    def _encode_vector(self, entries):
     387        return " ".join(str(entry) for entry in entries)
     388
     389    def _encode_bool(self, value):
     390        return str(int(value))
     391
     392    def _encode_rational(self, value):
     393        return str(value)
     394
     395    def _encode_int(self, value):
     396        return str(value)
     397       
     398    def _get_bool_prop(self, prop):
     399        return self._parse_bool(getattr(self._polymake_obj, prop))
     400
     401    def _get_rational_prop(self, prop):
     402        return self._parse_rational(getattr(self._polymake_obj, prop))
     403
     404    def _get_int_prop(self, prop):
     405        return self._parse_int(getattr(self._polymake_obj, prop))
     406
     407    def _get_matrix_prop(self, prop):
     408        return matrix(self._parse_matrix(getattr(self._polymake_obj, prop)))
     409
     410    def _get_vector_list_prop(self, prop):
     411        return [sage_vector(row) for row in self._parse_matrix(getattr(self._polymake_obj, prop))]
     412
     413    def _set_bool_prop(self, prop, value):
     414        return setattr(self._polymake_obj, prop, self._encode_bool(value))
     415
     416    def _set_rational_prop(self, prop, value):
     417        return setattr(self._polymake_obj, prop, self._encode_rational(value))
     418
     419    def _set_int_prop(self, prop, value):
     420        return setattr(self._polymake_obj, prop, self._encode_int(value))
     421
     422    def _set_matrix_prop(self, prop, value):
     423        return setattr(self._polymake_obj, prop, value)
     424
     425    def affine_hull(self):
     426        return self._get_matrix_prop("AFFINE_HULL")
     427
     428    def is_bounded(self):
     429        return self._get_bool_prop("BOUNDED")
     430
     431    def is_centered(self):
     432        return self._get_bool_prop("CENTERED")
     433
     434    def is_feasible(self):
     435        return self._get_bool_prop("FEASIBLE")
     436
     437    def is_simple(self):
     438        return self._get_bool_prop("SIMPLE")
     439
     440    def gale_transform(self):
     441        return self._get_matrix_prop("GALE_TRANSFORM")
     442 
     443    def n_points(self):
     444        return self._get_int_prop("N_POINTS")
     445
     446    def _points(self):
     447        return self._get_vector_list_prop("POINTS")
     448
     449    def vertices(self):
     450        return self._get_vector_list_prop("VERTICES")
     451
     452    def facets(self):
     453        return self._get_vector_list_prop("FACETS")
     454
     455    def n_facets(self):
     456        return self._get_int_prop("N_FACETS")
     457
     458    def visual(self):
     459        tmp_file = tempfile.NamedTemporaryFile()
     460        # FIXME: polymake doesn't use our temporary file, but a file
     461        # with suffix .poly instead. That's probably nonexistent and
     462        # writeable so we should be fine
     463        self.write(tmp_file.name)
     464        result_filename = tmp_file.name + ".poly"
     465        subprocess.Popen('polymake "%s" VISUAL; rm "%s"' % (result_filename, result_filename), shell=True)
     466
     467    def equations(self):
     468        try:
     469            return self._get_vector_list_prop("EQUATIONS")
     470        except:
     471            return []
     472   
     473    def _set_inequalities_coefficients(self, inequalities_coefficients):
     474        return self._set_matrix_prop("INEQUALITIES", inequalities_coefficients)
     475
     476    def _set_equality_coefficients(self, equality_coefficients):
     477        return self._set_matrix_prop("EQUATIONS", equality_coefficients)
     478
     479    def set_defining_equations(self, equations, homogeneous_coordinates, affine_plane_equation):
     480        self._homogeneous_coordinates = homogeneous_coordinates
     481        self._affine_plane_equation = affine_plane_equation
     482        symbol_coordinates = [c for c in homogeneous_coordinates if hasattr(c, 'is_symbol') and c.is_symbol()]
     483        if len(symbol_coordinates) < len(homogeneous_coordinates) - 1:
     484            raise ValueError("Too many non-symbol coordinates")
     485        if len(symbol_coordinates) == len(homogeneous_coordinates) - 1:
     486            if affine_plane_equation != None:
     487                raise ValueError("Cannot specify both affine coordinates and an affine plane equation")
     488            d = SR.symbol()
     489            homogeneous_coordinates = [c if hasattr(c, 'is_symbol') and c.is_symbol() else d for c in homogeneous_coordinates]
     490            affine_plane_equation = (d == 1)
     491        # solve the affine plane equation to find a homogeneous representation of 1
     492        if affine_plane_equation != None:
     493            zero_expr = affine_plane_equation.lhs() - affine_plane_equation.rhs()
     494            const_coeff = QQ(zero_expr.subs({var:0 for var in homogeneous_coordinates}))
     495            homogeneous_one = (const_coeff - zero_expr) / const_coeff
     496            if homogeneous_one.subs({var:0 for var in homogeneous_coordinates}) != 0:
     497                raise ValueError("Cannot represent affine_plane equation homogeneously")
     498           
     499        inequalities_coefficients = []
     500        equalities_coefficients = []
     501        for expr in equations:
     502            if expr is True:
     503                continue
     504            elif expr is False:
     505                inequalities_coefficients.append([-1] + [0 for _ in homogeneous_coordinates][1:])
     506                continue
     507            elif expr.operator() in (operator.le, operator.lt):
     508                zero_expr = expr.rhs() - expr.lhs()
     509            elif expr.operator() in (operator.ge, operator.gt):
     510                zero_expr = expr.lhs() - expr.rhs()
     511            elif expr.operator() == operator.eq:
     512                zero_expr = expr.lhs() - expr.rhs()
     513            else:
     514                raise ValueError("equations should have operator >, <, >=, <= or ==")
     515
     516            if affine_plane_equation != None:
     517                const_coeff = QQ(zero_expr.subs({var:0 for var in homogeneous_coordinates}))
     518                zero_expr = zero_expr - const_coeff + const_coeff * homogeneous_one
     519
     520            coeffs = [QQ(zero_expr.coefficient(var)) for var in homogeneous_coordinates]
     521            if sum(coeff*var for coeff, var in zip(coeffs, homogeneous_coordinates)) != zero_expr:
     522                raise ValueError("equation is not homogeneous or affine")
     523
     524            if expr.operator() == operator.eq:
     525                equalities_coefficients.append(coeffs)
     526            else:
     527                inequalities_coefficients.append(coeffs)
     528        self._set_inequalities_coefficients(inequalities_coefficients)
     529        self._set_equality_coefficients(equalities_coefficients)
     530
     531    def get_defining_equations(self):
     532        ineq_coeffs = self.facets()
     533        if self._homogeneous_coordinates is None:
     534            self._homogeneous_coordinates = [SR.symbol() for _ in ineq_coeffs[0]]
     535        ineqs = [sum(c*var for c,var in zip(coeff, self._homogeneous_coordinates)) > 0
     536                    for coeff in ineq_coeffs]
     537        eq_coeffs = self.equations()
     538        eqs = [sum(c*var for c,var in zip(coeff, self._homogeneous_coordinates)) == 0
     539                    for coeff in eq_coeffs]
     540        return eqs + ineqs
     541
     542    def graph(self):
     543        return getattr(self._polymake_obj, "GRAPH")
     544
     545    def volume(self):
     546        return self._get_rational_prop("VOLUME")
     547
     548    def _binary_operation(self, other, name):
     549        retval = Polytope(call_function(name, self._polymake_obj, other._polymake_obj))
     550        retval._homogeneous_coordinates = self._homogeneous_coordinates
     551        retval._affine_plane_equation = self._affine_plane_equation
     552        return retval
     553
     554    def intersection(self, other):
     555        return self._binary_operation(other, "intersection")
     556
     557    def join(self, other):
     558        return self._binary_operation(other, "join_polytopes")
     559
     560    def __eq__(self, other):
     561        return call_function("equal_polyhedra", self._polymake_obj, other._polymake_obj, return_type=bool)
     562
     563    def congruent(self, other):
     564        # FIXME: we'd like to return the Rational that gives the congruence scale factor,
     565        # but I haven't found a way to Cythonize a call to numerator(Rational) etc.
     566        return call_function("equal_polyhedra", self._polymake_obj, other._polymake_obj, return_type=bool)
     567
     568
     569#constructor
     570polytope = Polytope
     571
     572class Polymake:
     573    """
     574    This class contains members to generate predefined polytopes
     575    """
     576    def associahedron(self, dimension):
     577        return Polytope(call_function("associahedron", int(dimension)), name="%s-dimensional associahedron" % dimension)
     578
     579    def birkhoff(self, n, even=True):
     580        # FIXME: should even default to True or False for backward compat?
     581        return Polytope(call_function("birkhoff", int(n), bool(even)), name="%s Birkhoff %s" % ("even" if even else "odd", n))
     582
     583    def cell24(self):
     584        return Polytope(call_function("create_24_cell"), name="The 24-cell")
     585
     586    def convex_hull(self, points=[]):
     587        points.sort()
     588        retval = Polytope(name="Convex hull of points %s" % points)
     589        retval._set_matrix_prop("POINTS", points)
     590        return retval
     591
     592    def cube(self, dimension, scale=0):
     593        return Polytope(call_function("cube", int(dimension), QQ(scale)), name="Cube of dimension %s (scale %s)" % (dimension, scale))
     594
     595    def from_data(self, data):
     596        raise NotImplementedError
     597
     598    def rand01(self, d, n, seed=None):
     599        raise NotImplementedError
     600
     601polymake = Polymake()
     602
     603todo = """MINIMAL_VERTEX_ANGLE: common::Float
     604 
     605ONE_VERTEX: common::Vector<Scalar>
     606 
     607STEINER_POINTS: common::Matrix<Scalar, NonSymmetric>
     608 
     609VALID_POINT: common::Vector<Scalar>
     610 
     611VERTEX_BARYCENTER: common::Vector<Scalar>
     612 
     613VERTEX_LABELS: common::Array<String>
     614 
     615VERTEX_NORMALS: common::Matrix<Scalar, NonSymmetric>
     616 
     617ZONOTOPE_INPUT_VECTORS: common::Matrix<Scalar, NonSymmetric>
     618    def _
     619    """
     620
     621# "none" signifies that we don't want polymake to save user configuration data
     622main = _Main("none")
     623main.set_application("polytope")
     624
  • new file sage/geometry/polymake/polymake_wrappers.h

    diff --git a/sage/geometry/polymake/polymake_wrappers.h b/sage/geometry/polymake/polymake_wrappers.h
    new file mode 100644
    - +  
     1#include "polymake/client.h"
     2
     3using namespace pm;
     4
     5inline void VoidCallPolymakeFunction_wrapper(std::string name) {
     6    VoidCallPolymakeFunction(name);
     7}
     8inline void VoidCallPolymakeFunction_wrapper(std::string name, perl::Object ob1) {
     9    VoidCallPolymakeFunction(name, ob1);
     10}
     11inline void VoidCallPolymakeFunction_wrapper(std::string name, perl::Object ob1, perl::Object ob2) {
     12    VoidCallPolymakeFunction(name, ob1, ob2);
     13}
     14inline void VoidCallPolymakeFunction_wrapper(std::string name, perl::Object ob1, perl::Object ob2, perl::Object ob3) {
     15    VoidCallPolymakeFunction(name, ob1, ob2, ob3);
     16}
     17
     18inline perl::Object *CallPolyMakeFunction_wrapper(std::string name) {
     19    perl::Object *retval = new perl::Object();
     20    CallPolymakeFunction(name) >> *retval;
     21    return retval;
     22}
     23inline perl::Object *CallPolymakeFunction_wrapper(std::string name, perl::Object ob1) {
     24    perl::Object *retval = new perl::Object();
     25    CallPolymakeFunction(name, ob1) >> *retval;
     26    return retval;
     27
     28}
     29inline perl::Object *CallPolymakeFunction_wrapper(std::string name, perl::Object ob1, perl::Object ob2) {
     30    perl::Object *retval = new perl::Object();
     31    CallPolymakeFunction(name, ob1, ob2) >> *retval;
     32    return retval;
     33
     34}
     35
     36inline perl::Object *CallPolymakeFunction_wrapper(std::string name, perl::Object ob1, perl::Object ob2, perl::Object ob3) {
     37    perl::Object *retval = new perl::Object();
     38    CallPolymakeFunction(name, ob1, ob2, ob3) >> *retval;
     39    return retval;
     40
     41}
     42
     43inline perl::Object *CallPolymakeFunction_wrapper(std::string name) {
     44    perl::Object *retval = new perl::Object();
     45    CallPolymakeFunction(name) >> *retval;
     46    return retval;
     47}
     48
     49inline perl::Object *CallPolymakeFunction_wrapper(std::string name, int value) {
     50    perl::Object *retval = new perl::Object();
     51    CallPolymakeFunction(name, value) >> *retval;
     52    return retval;
     53}
     54
     55inline perl::Object *CallPolymakeFunction_wrapper(std::string name, int value1, bool value2) {
     56    perl::Object *retval = new perl::Object();
     57    CallPolymakeFunction(name, value1, value2) >> *retval;
     58    return retval;
     59}
     60
     61inline perl::Object *CallPolymakeFunction_wrapper(std::string name, int value1, Rational value2) {
     62    perl::Object *retval = new perl::Object();
     63    CallPolymakeFunction(name, value1, value2) >> *retval;
     64    return retval;
     65}
     66
     67inline bool boolCallPolymakeFunction_wrapper(std::string name, perl::Object ob1, perl::Object ob2) {
     68    bool retval;
     69    CallPolymakeFunction(name, ob1, ob2) >> retval;
     70    return retval;
     71
     72}
     73
     74inline Rational rationalCallPolymakeFunction_wrapper(std::string name, perl::Object ob1, perl::Object ob2) {
     75    Rational retval;
     76    CallPolymakeFunction(name, ob1, ob2) >> retval;
     77    return retval;
     78
     79}
     80
     81inline void set_property_wrapper(perl::Object object, std::string propertyname, Matrix<Rational> &matrix) {
     82    object.take(propertyname) << matrix;
     83}
     84
     85inline void set_property_wrapper(perl::Object object, std::string propertyname, std::string value) {
     86    object.take(propertyname) << value;
     87}
     88
     89inline std::string get_property_wrapper(perl::Object object, std::string propertyname) {
     90    std::string retval;
     91    object.give(propertyname) >> retval;
     92    return retval;
     93}
  • sage/geometry/polytope.py

    diff --git a/sage/geometry/polytope.py b/sage/geometry/polytope.py
    a b  
    11"""
    22Polytopes
    33
    4 This module provides access to polymake, which 'has been developed
    5 since 1997 in the Discrete Geometry group at the Institute of
    6 Mathematics of Technische Universitat Berlin. Since 2004 the
    7 development is shared with Fachbereich Mathematik, Technische
    8 Universitat Darmstadt. The system offers access to a wide variety
    9 of algorithms and packages within a common framework. polymake is
    10 flexible and continuously expanding. The software supplies C++ and
    11 Perl interfaces which make it highly adaptable to individual
    12 needs.'
     4This file is a wrapper for the polymake module, which cannot be
     5imported by default because it is linked to an optional shared
     6library
    137
    14 .. note::
    15 
    16    If you have trouble with this module do::
    17 
    18        sage: !polymake --reconfigure   # not tested
    19 
    20    at the command line.
    21 
    22 AUTHORS:
    23 
    24 - Ewgenij Gawrilow, Michael Joswig: main authors of polymake
    25 
    26 - William Stein: Sage interface
    278"""
    289
    2910########################################################################
    30 #       Copyright (C) 2006 William Stein <wstein@gmail.com>
     11#       Copyright (C) 2012 Timo Kluck <tkluck@infty.nl>
     12#       Copyright (C) 2012 William Stein <wstein@gmail.com>
    3113#
    3214#  Distributed under the terms of the GNU General Public License (GPL)
    3315#
    3416#                  http://www.gnu.org/licenses/
    3517########################################################################
    3618
     19class PolymakeProxy(object):
     20    # TODO: use some wrapper decorator to make sure that this proxy has
     21    # the right tab completion and docstrings
     22    def __getattr__(self, attr):
     23        try:
     24            import sage.geometry.polymake as pm
     25        except:
     26            raise RuntimeError(
     27            """Polymake requires the optional polymake package. It can be installed
     28            by typing
    3729
    38 from sage.misc.all import SAGE_TMP, tmp_filename
    39 from sage.rings.all import Integer, QQ
    40 from sage.structure.all import Sequence
    41 from sage.modules.all import VectorSpace
    42 from sage.structure.sage_object import SageObject
     30                 sage -i polymake
    4331
    44 import os
    45 from subprocess import *
     32            on the command line""")
     33        return getattr(pm.polymake, attr)
     34           
    4635
    47 path = '%s/local/polymake/bin/'%os.environ['SAGE_ROOT']
    48 polymake_command = path + 'polymake'
    49 
    50 if os.path.exists(path):
    51     os.environ['PATH'] = '%s:'%path + os.environ['PATH']
    52 
    53 tmp_file = os.path.join(SAGE_TMP, 'tmp.poly')
    54 
    55 class Polytope(SageObject):
    56     """
    57     Create a polytope.
    58 
    59     EXAMPLES::
    60 
    61         sage: P = polymake.convex_hull([[1,0,0,0], [1,0,0,1], [1,0,1,0], [1,0,1,1],  [1,1,0,0], [1,1,0,1], [1,1,1,0], [1,1,1,1]])   # optional - polymake
    62 
    63     .. note::
    64 
    65        If you have trouble with this module do::
    66 
    67           sage: !polymake --reconfigure   # not tested
    68 
    69        at the command line.
    70     """
    71     def __init__(self, datafile, desc):
    72         self.__data = datafile
    73         self.__desc = desc
    74 
    75     def _repr_(self):
    76         s = self.__desc
    77         # vertices, facets, points, inequalities
    78         try:
    79             s += '\nVertices:\n%s'%self.__vertices
    80         except AttributeError:
    81             pass
    82         try:
    83             s += '\nFacets:\n%s'%self.__facets
    84         except AttributeError:
    85             pass
    86         return s
    87 
    88 
    89 
    90     def __add__(self, other):
    91         if not isinstance(other, Polytope):
    92             raise TypeError, "other (=%s) must be a polytope"%other
    93         output_file = tmp_filename()
    94         infile1 = tmp_filename()
    95         open(infile1,'w').write(self.__data)
    96         infile2 = tmp_filename()
    97         open(infile2,'w').write(other.__data)
    98         cmd = "minkowski_sum %s 1 %s 1 %s"%(output_file, infile1,
    99                                             infile2)
    100         stdin, stdout, stderr = os.popen3(cmd)
    101         stdin.close()
    102         err = stderr.read()
    103         if len(err) > 0:
    104             raise RuntimeError, err
    105         print stdout.read(), err
    106         S = polymake.from_data(open(output_file).read())
    107         os.unlink(infile1)
    108         os.unlink(infile2)
    109         os.unlink(output_file)
    110         return S
    111 
    112 
    113     def data(self):
    114         return self.__data
    115 
    116     def write(self, filename):
    117         open(filename,'w').write(self.__data)
    118 
    119     def cmd(self, cmd):
    120         cmd = cmd.upper()
    121         # First check if the value of the command
    122         # is already known.
    123         D = self.data()
    124         cmd2='\n%s\n'%cmd
    125         i = D.find(cmd2)
    126         if i != -1:
    127             j = D[i:].find('\n\n')
    128             if j == -1:
    129                 j = len(D)
    130             else:
    131                 j += i
    132             return D[i+len(cmd2)-1:j]
    133 
    134         F = tmp_file
    135         open(F,'w').write(self.__data)
    136         c = '%s %s %s'%(polymake_command, F, cmd)
    137         polymake_processes = Popen([polymake_command, F, cmd],stdout=PIPE,stderr=PIPE)
    138         ans, err = polymake_processes.communicate()
    139         if len(err) > 0:
    140             raise RuntimeError, err
    141         if len(ans) == 0:
    142             raise ValueError, "%s\nError executing polymake command %s"%(
    143                 err,cmd)
    144         self.__data = open(F).read()
    145         return ans
    146 
    147 
    148     def facets(self):
    149         """
    150         EXAMPLES::
    151 
    152             sage: P = polymake.convex_hull([[1,0,0,0], [1,0,0,1], [1,0,1,0], [1,0,1,1],  [1,1,0,0], [1,1,0,1], [1,1,1,0], [1,1,1,1]])   # optional - polymake
    153             sage: P.facets()                            # optional - polymake
    154             [(0, 0, 0, 1), (0, 1, 0, 0), (0, 0, 1, 0), (1, 0, 0, -1), (1, 0, -1, 0), (1, -1, 0, 0)]
    155         """
    156         try:
    157             return self.__facets
    158         except AttributeError:
    159             pass
    160         s = self.cmd('FACETS')
    161         s = s.rstrip().split('\n')[1:]
    162         if len(s) == 0:
    163             ans = Sequence([], immutable=True)
    164         else:
    165             n = len(s[0].split())
    166             V = VectorSpace(QQ, n)
    167             ans = Sequence((V(x.split()) for x in s), immutable=True)
    168         self.__facets = ans
    169         return ans
    170 
    171     def vertices(self):
    172         """
    173         EXAMPLES::
    174 
    175             sage: P = polymake.convex_hull([[1,0,0,0], [1,0,0,1], [1,0,1,0], [1,0,1,1],  [1,1,0,0], [1,1,0,1], [1,1,1,0], [1,1,1,1]])     # optional - polymake
    176             sage: P.vertices()                            # optional - polymake
    177             [(1, 0, 0, 0), (1, 0, 0, 1), (1, 0, 1, 0), (1, 0, 1, 1), (1, 1, 0, 0), (1, 1, 0, 1), (1, 1, 1, 0), (1, 1, 1, 1)]
    178         """
    179         try:
    180             return self.__vertices
    181         except AttributeError:
    182             pass
    183         s = self.cmd('VERTICES')
    184         s = s.rstrip().split('\n')[1:]
    185         if len(s) == 0:
    186             ans = Sequence([], immutable=True)
    187         else:
    188             n = len(s[0].split())
    189             V = VectorSpace(QQ, n)
    190             ans = Sequence((V(x.split()) for x in s), immutable=True)
    191         self.__vertices = ans
    192         return ans
    193 
    194     def visual(self):
    195         try:
    196             self.cmd('visual')
    197         except ValueError:
    198             pass
    199 
    200     def graph(self):
    201         try:
    202             return self.__graph
    203         except AttributeError:
    204             pass
    205         g = self.cmd('GRAPH')
    206         return g
    207 
    208     def is_simple(self):
    209         r"""
    210         Return True if this polytope is simple.
    211 
    212         A polytope is *simple* if the degree of each vertex equals the
    213         dimension of the polytope.
    214 
    215         EXAMPLES::
    216 
    217             sage: P = polymake.convex_hull([[1,0,0,0], [1,0,0,1], [1,0,1,0], [1,0,1,1],  [1,1,0,0], [1,1,0,1], [1,1,1,0], [1,1,1,1]])        # optional - polymake
    218             sage: P.is_simple()                              # optional - polymake
    219             True
    220 
    221         AUTHORS:
    222 
    223         - Edwin O'Shea (2006-05-02): Definition of simple.
    224         """
    225         try:
    226             return self.__is_simple
    227         except AttributeError:
    228             pass
    229         s = self.cmd('SIMPLE')
    230         i = s.find('\n')
    231         self.__is_simple = bool(int(s[i:]))
    232         return self.__is_simple
    233 
    234 
    235 
    236     def n_facets(self):
    237         """
    238         EXAMPLES::
    239 
    240             sage: P = polymake.convex_hull([[1,0,0,0], [1,0,0,1], [1,0,1,0], [1,0,1,1],  [1,1,0,0], [1,1,0,1], [1,1,1,0], [1,1,1,1]])     # optional - polymake
    241             sage: P.n_facets()                            # optional - polymake
    242             6
    243         """
    244         try:
    245             return self.__n_facets
    246         except AttributeError:
    247             pass
    248         s = self.cmd('N_FACETS')
    249         i = s.find('\n')
    250         self.__n_facets = Integer(s[i:])
    251         return self.__n_facets
    252 
    253 class Polymake:
    254     def __repr__(self):
    255         return "Object that makes polytopes."
    256 
    257     def __make(self, cmd, name):
    258         os.system(cmd)
    259         try:
    260             d = open(tmp_file).read()
    261         except IOError:
    262             raise RuntimeError, "You may need to install the polymake package"
    263         return Polytope(d, name)
    264 
    265     def reconfigure(self):
    266         """
    267         Reconfigure polymake.
    268 
    269         Remember to run polymake.reconfigure() as soon as you have changed
    270         the customization file and/or installed missing software!
    271         """
    272         os.system("polymake --reconfigure")
    273 
    274     def associahedron(self, dimension):
    275         return self.__make('associahedron %s %s'%(tmp_file, dimension),
    276                            '%s-dimensional associahedron'%dimension)
    277 
    278     def birkhoff(self, n):
    279         return self.__make('birkhoff %s %s'%(tmp_file, n),
    280                            'Birkhoff %s'%n)
    281 
    282 
    283     def cell24(self):
    284         """
    285         EXAMPLES::
    286 
    287             sage: polymake.cell24()            # optional - polymake
    288             The 24-cell
    289         """
    290         return self.__make('24-cell %s'%tmp_file,
    291                            'The 24-cell')
    292 
    293     def convex_hull(self, points=[]):
    294         r"""
    295         EXAMPLES::
    296 
    297             sage: R.<x,y,z> = PolynomialRing(QQ,3)
    298             sage: f = x^3 + y^3 + z^3 + x*y*z
    299             sage: e = f.exponents()
    300             sage: a = [[1] + list(v) for v in e]
    301             sage: a
    302             [[1, 3, 0, 0], [1, 0, 3, 0], [1, 1, 1, 1], [1, 0, 0, 3]]
    303             sage: n = polymake.convex_hull(a)       # optional - polymake
    304             sage: n                                 # optional - polymake
    305             Convex hull of points [[1, 0, 0, 3], [1, 0, 3, 0], [1, 1, 1, 1], [1, 3, 0, 0]]
    306             sage: n.facets()                        # optional - polymake
    307             [(0, 1, 0, 0), (3, -1, -1, 0), (0, 0, 1, 0)]
    308             sage: n.is_simple()                     # optional - polymake
    309             True
    310             sage: n.graph()                         # optional - polymake
    311             'GRAPH\n{1 2}\n{0 2}\n{0 1}\n\n'
    312         """
    313         f = 'POINTS\n'
    314         points.sort()
    315         for p in points:
    316             f += ' '.join(str(x) for x in p) + '\n'
    317         f += '\n'
    318         return Polytope(f, 'Convex hull of points %s'%points)
    319 
    320     def cube(self, dimension, scale=0):
    321         return self.__make('cube %s %s %s'%(tmp_file, dimension, scale),
    322                            'Cube of dimension %s (scale %s)'%(dimension, scale))
    323 
    324     def from_data(self, data):
    325         return Polytope(data, 'A Polytope')
    326 
    327     def rand01(self, d, n, seed=None):
    328         cmd = 'rand01 %s %s %s'%(tmp_file, d, n)
    329         if not seed is None:
    330             cmd += ' -seed %s'%seed
    331         return self.__make(cmd,
    332               '%s-dimensional 0/1-polytope with %s random vertices (uniform distribution)'%(d, n))
    333 
    334 
    335 
    336 polymake = Polymake()
     36polymake = PolymakeProxy()