Ticket #9188: trac_9188_fix_facet_normal_reviewer.patch

File trac_9188_fix_facet_normal_reviewer.patch, 6.5 KB (added by novoselt, 12 years ago)

Apply on top of the original patch.

  • sage/geometry/lattice_polytope.py

    # HG changeset patch
    # User Andrey Novoseltsev <novoselt@gmail.com>
    # Date 1276029901 21600
    # Node ID 125efc24ad6b1db4e2e95c1ae54458f8c9307dfc
    # Parent  09dbbfb5c3f65c62405778b8d5c4b15e5da91771
    Trac 9188: Bug in facet_normal of lattice polytopes - reviewer changes.
    
    Added abs(...) to dual scaling and prettified doctests.
    
    diff -r 09dbbfb5c3f6 -r 125efc24ad6b sage/geometry/lattice_polytope.py
    a b  
    521521            self._shift_vector = p0
    522522            if compute_vertices:
    523523                self._vertices = self._embed(H_polytope.vertices())
    524 
     524            # In order to use facet normals obtained from subpolytopes, we
     525            # need the following (see Trac #9188).
    525526            M = self._embedding_matrix
    526             basis = matrix( M.columns() + M.integer_kernel().basis() ).transpose()
    527             self._dual_embedding_scale = basis.det()
     527            # Basis for the ambient space with spanned subspace in front
     528            basis = M.columns() + M.integer_kernel().basis()
     529            # Let's represent it as columns of a matrix
     530            basis = matrix(basis).transpose()
     531            # Absolute value helps to keep normals "inner"
     532            self._dual_embedding_scale = abs(basis.det())
    528533            dualbasis = matrix(ZZ, self._dual_embedding_scale * basis.inverse())
    529534            self._dual_embedding_matrix = dualbasis.submatrix(0,0,M.ncols())
    530535       
     
    15131518            sage: p.facet_constant(0)
    15141519            1       
    15151520
    1516         A lattice polytope of strictly smaller dimension (=2) than the
    1517         ambient dimension(=4)::
    1518 
    1519             sage: LatticePolytope(matrix([(1, -1, 0), (-1, -1, 0), (1, 1, 0), (3, 3, 0)]))
     1521        This is a 2-dimensional lattice polytope in a 4-dimensional space::
     1522
     1523            sage: m = matrix([(1,-1,1,3), (-1,-1,1,3), (0,0,0,0)])
     1524            sage: p = LatticePolytope(m.transpose())
     1525            sage: p
    15201526            A lattice polytope: 2-dimensional, 3 vertices.
    1521             sage: lp = LatticePolytope(matrix([(1, -1, 0), (-1, -1, 0), (1, 1, 0), (3, 3, 0)]))
    1522             sage: lp.vertices()
     1527            sage: p.vertices()
    15231528            [ 1 -1  0]
    15241529            [-1 -1  0]
    15251530            [ 1  1  0]
    15261531            [ 3  3  0]
    1527             sage: matrix([[ lp.facet_normal(i)*lp.vertex(j) + lp.facet_constant(i) for i in range(0,3)] for j in range(0,3)])
    1528             [-22   0   0]
    1529             [  0 -22   0]
    1530             [  0   0 -11]
     1532            sage: fns = [p.facet_normal(i) for i in range(p.nfacets())]
     1533            sage: fns
     1534            [(11, -1, 1, 3), (-11, -1, 1, 3), (0, 1, -1, -3)]
     1535            sage: fcs = [p.facet_constant(i) for i in range(p.nfacets())]
     1536            sage: fcs
     1537            [0, 0, 11]
     1538           
     1539        Now we manually compute the distance matrix of this polytope. Since it
     1540        is a triangle, each line (corresponding to a facet) should have two
     1541        zeros (vertices of the corresponding facet) and one positive number
     1542        (since our normals are inner)::
     1543       
     1544            sage: matrix([[fns[i] * p.vertex(j) + fcs[i]
     1545            ...            for j in range(p.nvertices())]
     1546            ...           for i in range(p.nfacets())])
     1547            [22  0  0]
     1548            [ 0 22  0]
     1549            [ 0  0 11]
    15311550        """
    15321551        if self.is_reflexive():
    15331552            return 1
    15341553        elif self.ambient_dim() == self.dim():
    15351554            return self._facet_constants[i]
    15361555        else:
    1537             return (self._sublattice_polytope.facet_constant(i) * self._dual_embedding_scale
    1538                     - self.facet_normal(i) * self._shift_vector )
     1556            return (self._sublattice_polytope.facet_constant(i)
     1557                    * self._dual_embedding_scale
     1558                    - self.facet_normal(i) * self._shift_vector)
    15391559       
    15401560    def facet_normal(self, i):
    15411561        r"""
    15421562        Return the inner normal to the ``i``-th facet of this polytope.
    15431563       
    15441564        If this polytope is not full-dimensional, facet normals will be
    1545         parallel to the affine subspace spanned by this polytope.
     1565        orthogonal to the integer kernel of the affine subspace spanned by
     1566        this polytope.
    15461567       
    15471568        INPUT:
    15481569       
     
    15881609            sage: p.facet_constant(0)
    15891610            1       
    15901611
    1591         Here is an example of a 3-dimensional polytope in 4d::
    1592 
    1593             sage: lp = LatticePolytope(matrix([[1,1,-1,0],[1,-1,-1,0],[1,1,1,0],[3,3,3,0]]))
    1594             sage: lp.vertices()
    1595             [ 1  1 -1  0]
    1596             [ 1 -1 -1  0]
    1597             [ 1  1  1  0]
    1598             [ 3  3  3  0]
    1599             sage: ker = lp.vertices().integer_kernel().matrix()
     1612        Here is an example of a 3-dimensional polytope in a 4-dimensional
     1613        space::
     1614
     1615            sage: m = matrix([(0,0,0,0), (1,1,1,3), (1,-1,1,3), (-1,-1,1,3)])
     1616            sage: p = LatticePolytope(m.transpose())
     1617            sage: p.vertices()
     1618            [ 0  1  1 -1]
     1619            [ 0  1 -1 -1]
     1620            [ 0  1  1  1]
     1621            [ 0  3  3  3]
     1622            sage: ker = p.vertices().integer_kernel().matrix()
    16001623            sage: ker
    16011624            [ 0  0  3 -1]
    1602             sage: [ ker * lp.facet_normal(i) for i in range(0,4) ]
     1625            sage: [ker * p.facet_normal(i) for i in range(p.nfacets())]
    16031626            [(0), (0), (0), (0)]
    1604             sage: matrix([lp.facet_normal(i) for i in range(0,4)]) * lp.vertices()
    1605             [  0   0 -20   0]
    1606             [  0 -20   0   0]
    1607             [ 10  10  10   0]
    1608             [-20   0   0   0]
    1609             sage: matrix([[ lp.facet_normal(i)*lp.vertex(j) + lp.facet_constant(i) for i in range(0,4)] for j in range(0,4)])
    1610             [  0   0   0 -20]
    1611             [  0 -20   0   0]
    1612             [-20   0   0   0]
    1613             [  0   0 -10   0]
     1627           
     1628        Now we manually compute the distance matrix of this polytope. Since it
     1629        is a simplex, each line (corresponding to a facet) should consist of
     1630        zeros (indicating generating vertices of the corresponding facet) and
     1631        a single positive number (since our normals are inner)::
     1632       
     1633            sage: matrix([[p.facet_normal(i) * p.vertex(j)
     1634            ...            + p.facet_constant(i)
     1635            ...            for j in range(p.nvertices())]
     1636            ...           for i in range(p.nfacets())])
     1637            [ 0  0  0 20]
     1638            [ 0  0 20  0]
     1639            [ 0 20  0  0]
     1640            [10  0  0  0]
    16141641        """
    16151642        if self.is_reflexive():
    16161643            return self.polar().vertex(i)