Ticket #14116: trac_14116-polymake-upgrade-rebased.patch

File trac_14116-polymake-upgrade-rebased.patch, 41.0 KB (added by kcrisman, 9 years ago)
  • module_list.py

    # HG changeset patch
    # User Timo Kluck <tkluck@infty.nl>
    # Date 1356283340 0
    # Node ID fe02b73f7bc5351589ae2bf506decdd48588ddeb
    # Parent  c73851d094ac59f9793ad6006ea4678286accbfc
    Add wrappers for new version of polymake
    
    diff --git a/module_list.py b/module_list.py
    a b  
    19861986                  libraries = ["csage", "stdc++", "gurobi"])
    19871987        )
    19881988
     1989if is_package_installed('polymake'):
     1990    ext_modules.append(
     1991        Extension('sage.geometry.polymake',
     1992                  sources = ["sage/geometry/polymake/polymake.pyx"],
     1993                  include_dirs = [SAGE_LOCAL + '/include'],
     1994                  libraries=["polymake","gmp","xml2"],
     1995                  library_dirs=[SAGE_LOCAL + "/lib"],
     1996                  extra_compile_args=["-DPOLYMAKE_DEBUG=0"],
     1997                  extra_link_args=["-Wl,-rpath,%s/lib" % SAGE_LOCAL],
     1998                  language="c++",
     1999                  depends = ["sage/geometry/polymake/polymake_wrappers.h"])
     2000        )
    19892001
    19902002if is_package_installed('coxeter3'):
    19912003    ext_modules.append(
  • 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 *
    46 from sage.env import SAGE_LOCAL
     32            on the command line""")
     33        return getattr(pm.polymake, attr)
     34           
    4735
    48 path = os.path.join(SAGE_LOCAL,'polymake','bin')
    49 polymake_command = path + 'polymake'
    50 
    51 if os.path.exists(path):
    52     os.environ['PATH'] = '%s:'%path + os.environ['PATH']
    53 
    54 tmp_file = os.path.join(SAGE_TMP, 'tmp.poly')
    55 
    56 class Polytope(SageObject):
    57     """
    58     Create a polytope.
    59 
    60     EXAMPLES::
    61 
    62         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
    63 
    64     .. note::
    65 
    66        If you have trouble with this module do::
    67 
    68           sage: !polymake --reconfigure   # not tested
    69 
    70        at the command line.
    71     """
    72     def __init__(self, datafile, desc):
    73         self.__data = datafile
    74         self.__desc = desc
    75 
    76     def _repr_(self):
    77         s = self.__desc
    78         # vertices, facets, points, inequalities
    79         try:
    80             s += '\nVertices:\n%s'%self.__vertices
    81         except AttributeError:
    82             pass
    83         try:
    84             s += '\nFacets:\n%s'%self.__facets
    85         except AttributeError:
    86             pass
    87         return s
    88 
    89 
    90 
    91     def __add__(self, other):
    92         if not isinstance(other, Polytope):
    93             raise TypeError, "other (=%s) must be a polytope"%other
    94         output_file = tmp_filename()
    95         infile1 = tmp_filename()
    96         open(infile1,'w').write(self.__data)
    97         infile2 = tmp_filename()
    98         open(infile2,'w').write(other.__data)
    99         cmd = "minkowski_sum %s 1 %s 1 %s"%(output_file, infile1,
    100                                             infile2)
    101         stdin, stdout, stderr = os.popen3(cmd)
    102         stdin.close()
    103         err = stderr.read()
    104         if len(err) > 0:
    105             raise RuntimeError, err
    106         print stdout.read(), err
    107         S = polymake.from_data(open(output_file).read())
    108         os.unlink(infile1)
    109         os.unlink(infile2)
    110         os.unlink(output_file)
    111         return S
    112 
    113 
    114     def data(self):
    115         return self.__data
    116 
    117     def write(self, filename):
    118         open(filename,'w').write(self.__data)
    119 
    120     def cmd(self, cmd):
    121         cmd = cmd.upper()
    122         # First check if the value of the command
    123         # is already known.
    124         D = self.data()
    125         cmd2='\n%s\n'%cmd
    126         i = D.find(cmd2)
    127         if i != -1:
    128             j = D[i:].find('\n\n')
    129             if j == -1:
    130                 j = len(D)
    131             else:
    132                 j += i
    133             return D[i+len(cmd2)-1:j]
    134 
    135         F = tmp_file
    136         open(F,'w').write(self.__data)
    137         c = '%s %s %s'%(polymake_command, F, cmd)
    138         polymake_processes = Popen([polymake_command, F, cmd],stdout=PIPE,stderr=PIPE)
    139         ans, err = polymake_processes.communicate()
    140         if len(err) > 0:
    141             raise RuntimeError, err
    142         if len(ans) == 0:
    143             raise ValueError, "%s\nError executing polymake command %s"%(
    144                 err,cmd)
    145         self.__data = open(F).read()
    146         return ans
    147 
    148 
    149     def facets(self):
    150         """
    151         EXAMPLES::
    152 
    153             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
    154             sage: P.facets()                            # optional - polymake
    155             [(0, 0, 0, 1), (0, 1, 0, 0), (0, 0, 1, 0), (1, 0, 0, -1), (1, 0, -1, 0), (1, -1, 0, 0)]
    156         """
    157         try:
    158             return self.__facets
    159         except AttributeError:
    160             pass
    161         s = self.cmd('FACETS')
    162         s = s.rstrip().split('\n')[1:]
    163         if len(s) == 0:
    164             ans = Sequence([], immutable=True)
    165         else:
    166             n = len(s[0].split())
    167             V = VectorSpace(QQ, n)
    168             ans = Sequence((V(x.split()) for x in s), immutable=True)
    169         self.__facets = ans
    170         return ans
    171 
    172     def vertices(self):
    173         """
    174         EXAMPLES::
    175 
    176             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
    177             sage: P.vertices()                            # optional - polymake
    178             [(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)]
    179         """
    180         try:
    181             return self.__vertices
    182         except AttributeError:
    183             pass
    184         s = self.cmd('VERTICES')
    185         s = s.rstrip().split('\n')[1:]
    186         if len(s) == 0:
    187             ans = Sequence([], immutable=True)
    188         else:
    189             n = len(s[0].split())
    190             V = VectorSpace(QQ, n)
    191             ans = Sequence((V(x.split()) for x in s), immutable=True)
    192         self.__vertices = ans
    193         return ans
    194 
    195     def visual(self):
    196         try:
    197             self.cmd('visual')
    198         except ValueError:
    199             pass
    200 
    201     def graph(self):
    202         try:
    203             return self.__graph
    204         except AttributeError:
    205             pass
    206         g = self.cmd('GRAPH')
    207         return g
    208 
    209     def is_simple(self):
    210         r"""
    211         Return True if this polytope is simple.
    212 
    213         A polytope is *simple* if the degree of each vertex equals the
    214         dimension of the polytope.
    215 
    216         EXAMPLES::
    217 
    218             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
    219             sage: P.is_simple()                              # optional - polymake
    220             True
    221 
    222         AUTHORS:
    223 
    224         - Edwin O'Shea (2006-05-02): Definition of simple.
    225         """
    226         try:
    227             return self.__is_simple
    228         except AttributeError:
    229             pass
    230         s = self.cmd('SIMPLE')
    231         i = s.find('\n')
    232         self.__is_simple = bool(int(s[i:]))
    233         return self.__is_simple
    234 
    235 
    236 
    237     def n_facets(self):
    238         """
    239         EXAMPLES::
    240 
    241             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
    242             sage: P.n_facets()                            # optional - polymake
    243             6
    244         """
    245         try:
    246             return self.__n_facets
    247         except AttributeError:
    248             pass
    249         s = self.cmd('N_FACETS')
    250         i = s.find('\n')
    251         self.__n_facets = Integer(s[i:])
    252         return self.__n_facets
    253 
    254 class Polymake:
    255     def __repr__(self):
    256         return "Object that makes polytopes."
    257 
    258     def __make(self, cmd, name):
    259         os.system(cmd)
    260         try:
    261             d = open(tmp_file).read()
    262         except IOError:
    263             raise RuntimeError, "You may need to install the polymake package"
    264         return Polytope(d, name)
    265 
    266     def reconfigure(self):
    267         """
    268         Reconfigure polymake.
    269 
    270         Remember to run polymake.reconfigure() as soon as you have changed
    271         the customization file and/or installed missing software!
    272         """
    273         os.system("polymake --reconfigure")
    274 
    275     def associahedron(self, dimension):
    276         return self.__make('associahedron %s %s'%(tmp_file, dimension),
    277                            '%s-dimensional associahedron'%dimension)
    278 
    279     def birkhoff(self, n):
    280         return self.__make('birkhoff %s %s'%(tmp_file, n),
    281                            'Birkhoff %s'%n)
    282 
    283 
    284     def cell24(self):
    285         """
    286         EXAMPLES::
    287 
    288             sage: polymake.cell24()            # optional - polymake
    289             The 24-cell
    290         """
    291         return self.__make('24-cell %s'%tmp_file,
    292                            'The 24-cell')
    293 
    294     def convex_hull(self, points=[]):
    295         r"""
    296         EXAMPLES::
    297 
    298             sage: R.<x,y,z> = PolynomialRing(QQ,3)
    299             sage: f = x^3 + y^3 + z^3 + x*y*z
    300             sage: e = f.exponents()
    301             sage: a = [[1] + list(v) for v in e]
    302             sage: a
    303             [[1, 3, 0, 0], [1, 0, 3, 0], [1, 1, 1, 1], [1, 0, 0, 3]]
    304             sage: n = polymake.convex_hull(a)       # optional - polymake
    305             sage: n                                 # optional - polymake
    306             Convex hull of points [[1, 0, 0, 3], [1, 0, 3, 0], [1, 1, 1, 1], [1, 3, 0, 0]]
    307             sage: n.facets()                        # optional - polymake
    308             [(0, 1, 0, 0), (3, -1, -1, 0), (0, 0, 1, 0)]
    309             sage: n.is_simple()                     # optional - polymake
    310             True
    311             sage: n.graph()                         # optional - polymake
    312             'GRAPH\n{1 2}\n{0 2}\n{0 1}\n\n'
    313         """
    314         f = 'POINTS\n'
    315         points.sort()
    316         for p in points:
    317             f += ' '.join(str(x) for x in p) + '\n'
    318         f += '\n'
    319         return Polytope(f, 'Convex hull of points %s'%points)
    320 
    321     def cube(self, dimension, scale=0):
    322         return self.__make('cube %s %s %s'%(tmp_file, dimension, scale),
    323                            'Cube of dimension %s (scale %s)'%(dimension, scale))
    324 
    325     def from_data(self, data):
    326         return Polytope(data, 'A Polytope')
    327 
    328     def rand01(self, d, n, seed=None):
    329         cmd = 'rand01 %s %s %s'%(tmp_file, d, n)
    330         if not seed is None:
    331             cmd += ' -seed %s'%seed
    332         return self.__make(cmd,
    333               '%s-dimensional 0/1-polytope with %s random vertices (uniform distribution)'%(d, n))
    334 
    335 
    336 
    337 polymake = Polymake()
     36polymake = PolymakeProxy()