Ticket #11564: trac_11564_preliminary_effort.patch

File trac_11564_preliminary_effort.patch, 12.1 KB (added by mhampton, 8 years ago)

Just for comparison with other effort, don't review

  • sage/geometry/polyhedra.py

    # HG changeset patch
    # User Marshall Hampton <hamptonio@gmail.com>
    # Date 1311268446 18000
    # Node ID 871b7c32e4561aab4e801a41513786645264d5d7
    # Parent  f69e26a8724f5e1ee34bcdbdaca6feaad72040b1
    For 11564, this is just a preliminary effort for comparison.
    
    diff -r f69e26a8724f -r 871b7c32e456 sage/geometry/polyhedra.py
    a b  
    15351535    projection_3d = Projection(polyhedron).schlegel(projection_direction)
    15361536    return render_3d(projection_3d, **kwds)
    15371537
    1538 
    15391538#########################################################################
    15401539class Polyhedron(SageObject):
    15411540    """
     
    17311730            cmdline_arg = '--reps'
    17321731        self._init_from_cdd_input(s, cmdline_arg, verbose)
    17331732
     1733        if self.dim() == 3:
     1734
     1735            def face_to_2d(self,a,b,a2,b2,face_index, dict_2d = {}, e2sign = 1):
     1736                '''
     1737
     1738                INPUT:
     1739
     1740                - ``p`` -- a polyhedron
     1741                - ``a`` -- index of a vertex
     1742                - ``b`` -- index of a second vertex
     1743                - ``a2`` -- 2d coordinates of a
     1744                - ``b2`` -- 2d coordinates of b
     1745                -  ``face_index`` -- index of the face we are laying down in 2d
     1746
     1747                OUTPUT:
     1748
     1749                - ``output`` - this is a modified-in-place version of dict_2d
     1750                '''
     1751                output = dict_2d
     1752                output[face_index] = {a:a2,b:b2}
     1753                vertices3d = self.vertices()
     1754                A = vector(RDF, vertices3d[a]) # 3d coordinates of vertex a
     1755                B = vector(RDF, vertices3d[b]) # 3d coordinates of vertex b
     1756                E1 = A - B
     1757                E1 = E1/norm(E1)
     1758                e1 = vector(a2) - vector(b2)
     1759                e1 = e1/norm(e1)
     1760                e2 = e2sign*vector([e1[1],-e1[0]])
     1761                face_vertex_indices = self.facial_incidences()[face_index][1][:]
     1762                face_vertex_indices.remove(a)
     1763                face_vertex_indices.remove(b)
     1764                for i in face_vertex_indices:
     1765                    # add in new 2d points, skipping a and b
     1766                    C = vector(vertices3d[i]) # 3d coordinates of vertex i
     1767                    dac = norm(A-C)
     1768                    cos_th = E1.inner_product(A-C)/dac
     1769                    output[face_index][i] = vector(a2) - e1*dac*cos_th + e2*dac*sqrt(1.0-cos_th**2)
     1770                return output
     1771
     1772            def find_e2sign(p,a,b,face_list,dict_2d):
     1773                '''
     1774
     1775                INPUT:
     1776
     1777                - ``p`` -- a polyhedron
     1778                - ``a`` -- 2d coordinate of unfolded vertex
     1779                - ``b`` -- 2d coordinate of unfolded vertex
     1780                - ``face_list`` -- list of adjacent faces
     1781                - ``dict_2d`` -- dictionary of face vertex indices to 2d coordinates
     1782
     1783                OUTPUT:
     1784
     1785                - Either 1 or -1 to indicate a direction relative to an ordered edge
     1786                '''
     1787                # a,b are coords of vertices
     1788                # face_list is the set of adjacent faces
     1789                oldface = face_list[0]
     1790                center = sum([vector(dict_2d[oldface][i]) for i in p.facial_incidences()[oldface][1]])
     1791                center = center/len(p.facial_incidences()[oldface][1])
     1792                e1 = vector(a) - vector(b)
     1793                e1 = e1/norm(e1)
     1794                e2 = vector([e1[1],-e1[0]])
     1795                if norm(e2+vector(a)-vector(center)) > norm(-e2+vector(a)-vector(center)):
     1796                    return 1
     1797                return -1
     1798
     1799            def facial_order(pin,algorithm='Kruskal'):
     1800                '''
     1801                Create a dictionary of adjacent faces in a minimal spanning tree.
     1802                '''
     1803                output = {}
     1804                pg = Graph(dict(pin.facial_adjacencies())).min_spanning_tree(algorithm=algorithm)
     1805                j = 0
     1806                for i in pg:
     1807                    output[j] = [i[0],i[1]]
     1808                    j = j+1
     1809                return output
     1810
     1811            def shared_vertices(p, a, b):
     1812                '''
     1813                Returns the vertices shared by an ridge (two faces with
     1814                indices a and b).
     1815                '''
     1816                fv = p.facial_incidences()[a][1]
     1817                sv = p.facial_incidences()[b][1]
     1818                fverts = [x for x in fv if x in sv]
     1819                return fverts
     1820
     1821            def shared_faces(p, a, b):
     1822                '''
     1823                Returns the faces shared by an edge (two vertices with
     1824                indices a and b).
     1825                '''
     1826                fv = p.vertex_incidences()[a][1]
     1827                sv = p.vertex_incidences()[b][1]
     1828                fverts = [x for x in fv if x in sv]
     1829                return fverts
     1830
     1831            def check_order(facialorder):
     1832                '''
     1833                Given a dictionary of face indices (keys are simply an ordering)
     1834                this returns them in an order that works for an unfolding.
     1835                '''
     1836                facelist = facialorder[0][:]
     1837                checkedorder = {0:facialorder[0][:]}
     1838                faceorder_todo = facialorder.keys()
     1839                pop = faceorder_todo.pop(0)
     1840                i = 1
     1841                j = 1
     1842                while faceorder_todo != []:
     1843                    for i in faceorder_todo:
     1844                        f1 = facialorder[i][0]
     1845                        f2 = facialorder[i][1]
     1846                        if f1 in facelist:
     1847                            facelist.append(f2)
     1848                            checkedorder[j] = facialorder[i]
     1849                            faceorder_todo.remove(i)
     1850                            j += 1
     1851                        elif f2 in facelist:
     1852                            f1,f2 = f2,f1
     1853                            facelist.append(f2)
     1854                            checkedorder[j] = [f1,f2]
     1855                            faceorder_todo.remove(i)
     1856                            j += 1
     1857                return checkedorder
     1858
     1859            def draw_net(pin = self, tree_algorithm = 'Kruskal', verbose = False, faceorder = ''):
     1860                '''
     1861                Returns a dictionary of 2-d unfolded faces given a 3-d polyhedron.
     1862                '''
     1863                if faceorder=='':
     1864                    faceorder = facial_order(pin, algorithm = tree_algorithm)
     1865                if verbose: print 'initial faceorder:', faceorder
     1866                faceorder = check_order(faceorder)
     1867                if verbose: print 'checked faceorder:',  faceorder
     1868                facial_order_graph = Graph(faceorder.values())
     1869                if verbose: facial_order_graph.show(graph_border=True)
     1870
     1871                v1,v2 = shared_vertices(pin, faceorder[0][0], faceorder[0][1])
     1872
     1873                #length of first edge:
     1874                edge1 = RDF(norm(vector(pin.vertices()[v1])-vector(pin.vertices()[v2])))
     1875
     1876                # enter the first face
     1877                dict2d = {}
     1878                dict2d = face_to_2d(pin, v1, v2, (0,0), (edge1, 0), 0, dict_2d = dict2d)
     1879
     1880                # laying down the rest of the faces
     1881                for i in range(0,max(faceorder.keys())+1):
     1882                    f1,f2 = faceorder[i]
     1883                    try:
     1884                        v1,v2 = shared_vertices(pin,faceorder[i][0], faceorder[i][1])
     1885                    except:
     1886                        print 'Problem at vertices,faces: ',v1,v2,i,faceorder[i]
     1887                        break
     1888                    es = find_e2sign(pin, dict2d[f1][v1], dict2d[f1][v2], [f1,f2], dict2d)
     1889                    dict2d = face_to_2d(pin, v1, v2, dict2d[f1][v1], dict2d[f1][v2], f2, dict_2d=dict2d, e2sign=es)
     1890
     1891                return dict2d
     1892
     1893            def steepestedge(p):
     1894                c = [random() for i in range(3)]
     1895                c = vector(c)
     1896                vs = p.vertices()
     1897                vs = [vector(x) for x in vs]
     1898                heights = [[c.inner_product(vs[i]),i] for i in range(len(vs))]
     1899                heights.sort()
     1900                vminus = heights[0][1]
     1901                vplus = heights[-1][1]
     1902                cutedges = []
     1903                for j,v in enumerate(vs):
     1904                    adjs = p.vertex_adjacencies()[j][1]
     1905                    steeps = [[(c.inner_product(v-vs[i])/(norm(v-vs[i]))),i] for i in adjs]
     1906                    cutedges.append([j,max(steeps)[1]])
     1907                return cutedges
     1908
     1909            def stedgeorder(pin):
     1910                '''
     1911                Returns a dictionary whose values are face indices, computed
     1912                from a steepest-edge cut.
     1913                '''
     1914                facialpairs = []
     1915                verticesset = []
     1916                cutedges = steepestedge(pin)
     1917
     1918                for i in range(len(pin.facial_adjacencies())):
     1919                    for j in pin.facial_adjacencies()[i][1]:
     1920                        if i<j:
     1921                            facialpairs.append([i,j])
     1922                for i in range(0, len(facialpairs)):
     1923                    verticesset.append(shared_vertices(pin, facialpairs[i][0],facialpairs[i][1]))
     1924                for x in cutedges:
     1925                    x.sort()
     1926                for x in verticesset:
     1927                    x.sort()
     1928
     1929                facialorder = {}
     1930                k=0
     1931                for i in verticesset:
     1932                    if i not in(cutedges):
     1933                        fp = shared_faces(pin, i[0],i[1])
     1934                        facialorder[k] = fp
     1935                        k+=1
     1936                facialorder = check_order(facialorder)       
     1937
     1938                return facialorder
     1939
     1940            def draw_edges(alpha = 1, thickness = 1, rgbcolor = 'black'):
     1941                '''
     1942                Returns the edges of the polyhedron unfolded by edge-cuts
     1943                (a polyhedral net) as a Graphics object.
     1944
     1945                EXAMPLES::
     1946
     1947                    sage: cube = polytopes.n_cube(3)
     1948                    sage: show(cube.draw_unfolding(rgbcolor='red'), aspect_ratio=1, axes = False)
     1949                '''
     1950                faces2di = draw_net()
     1951                edges = []
     1952                for v in self.vertex_adjacencies():
     1953                    v1 = v[0]
     1954                    for v2 in v[1]:
     1955                        if v2 > v1: edges.append([v1,v2])
     1956
     1957                edges2d = 0
     1958
     1959                for i in faces2di.keys():         
     1960                        edges2d += sum([line2d([faces2di[i][x],faces2di[i][y]], alpha=alpha, thickness=thickness, rgbcolor=rgbcolor) for x in faces2di[i].keys() for y in faces2di[i].keys() if [x,y] in edges])
     1961                return edges2d
     1962
     1963            self.draw_unfolding = draw_edges
    17341964
    17351965    def __getstate__(self):
    17361966        r"""
     
    62486478        verts = verts + permutations([0 for i in range(dim_n-1)] + [-1])
    62496479        return Polyhedron(vertices = verts)
    62506480
    6251    
    6252     def parallelotope(self, generators):
    6253         r"""
    6254         Return the parallelotope spanned by the generators.
    6255 
    6256         INPUT:
    6257        
    6258         - ``generators`` -- an iterable of anything convertible to vector
    6259           (for example, a list of vectors) such that the vectors all
    6260           have the same dimension.
    6261 
    6262         OUTPUT:
    6263 
    6264         The parallelotope. This is the multi-dimensional
    6265         generalization of a parallelogram (2 generators) and a
    6266         parallelepiped (3 generators).
    6267 
    6268         EXAMPLES::
    6269        
    6270             sage: polytopes.parallelotope([ (1,0), (0,1) ])
    6271             A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 4 vertices.
    6272             sage: polytopes.parallelotope([[1,2,3,4],[0,1,0,7],[3,1,0,2],[0,0,1,0]])
    6273             A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 16 vertices.
    6274         """
    6275         try:
    6276             generators = [ vector(QQ,v) for v in generators ]
    6277             field = QQ
    6278         except TypeError:
    6279             generators = [ vector(RDF,v) for v in generators ]
    6280             field = RDF
    6281            
    6282         from sage.combinat.combination import Combinations
    6283         par =  [ 0*generators[0] ]
    6284         par += [ sum(c) for c in Combinations(generators) if c!=[] ]
    6285         return Polyhedron(vertices=par, field=field)
    6286                            
    6287 
    62886481
    62896482polytopes = Polytopes()
    62906483