Ticket #10336: trac_10336_plot_non_strictly_convex_cones.patch

File trac_10336_plot_non_strictly_convex_cones.patch, 19.9 KB (added by novoselt, 12 years ago)

  • sage/geometry/cone.py 

    # HG changeset patch
# User Andrey Novoseltsev <novoselt@gmail.com>
# Date 1290975403 25200
# Node ID 2572cc23715e5f6f229844c4ee398a16777b43ce
# Parent  3e3e79e371f99a8ba6d9d3417c71584cf0560d47
Trac 10336: Allow plotting of non-strictly convex cones.

This patch also adds computation of the face lattice for non-strictly convex
cones and improves arrow scaling for plots of cones and fans.

diff -r 3e3e79e371f9 -r 2572cc23715e sage/geometry/cone.py
    a b  
    17601760        """
    17611761        if "_face_lattice" not in self.__dict__:
    17621762            if self._ambient is self:
    1763                 # We need to compute face lattice on our own
    1764                 if not self.is_strictly_convex():
    1765                     raise NotImplementedError("face lattice is currently "
    1766                                 "implemented only for strictly convex cones!")
    1767                 def ConeFace(rays, facets):
     1763                # We need to compute face lattice on our own. To accommodate
     1764                # non-strictly convex cones we split rays (or rather their
     1765                # indicies) into those in the linear subspace and others, which
     1766                # we refer to as atoms.
     1767                S = self.linear_subspace()
     1768                subspace_rays = []
     1769                atom_to_ray = []
     1770                atom_to_facets = []
     1771                normals = self.facet_normals()
     1772                facet_to_atoms = [[] for normal in normals]
     1773                for i, ray in enumerate(self):
     1774                    if ray in S:
     1775                        subspace_rays.append(i)
     1776                    else:
     1777                        facets = [j for j, normal in enumerate(normals)
     1778                                    if ray * normal == 0]
     1779                        atom_to_facets.append(facets)
     1780                        atom = len(atom_to_ray)
     1781                        for j in facets:
     1782                            facet_to_atoms[j].append(atom)
     1783                        atom_to_ray.append(i)
     1784                           
     1785                def ConeFace(atoms, facets):
    17681786                    if facets:
     1787                        rays = sorted([atom_to_ray[a] for a in atoms]
     1788                                      + subspace_rays)
    17691789                        face = ConvexRationalPolyhedralCone(
    17701790                                    ambient=self, ambient_ray_indices=rays)
    17711791                        # It may be nice if this functionality is exposed,
    … …  
    17751795                        return face
    17761796                    else:
    17771797                        return self
    1778                 ray_to_facets = []
    1779                 facet_to_rays = []
    1780                 for ray in self:
    1781                     ray_to_facets.append([])
    1782                 for normal in self.facet_normals():
    1783                     facet_to_rays.append([])
    1784                 for i, ray in enumerate(self):
    1785                     for j, normal in enumerate(self.facet_normals()):
    1786                         if ray * normal == 0:
    1787                             ray_to_facets[i].append(j)
    1788                             facet_to_rays[j].append(i)
     1798                       
    17891799                self._face_lattice = Hasse_diagram_from_incidences(
    1790                                     ray_to_facets, facet_to_rays, ConeFace)
     1800                                    atom_to_facets, facet_to_atoms, ConeFace)
    17911801            else:
    17921802                # Get face lattice as a sublattice of the ambient one
    17931803                allowed_indices = frozenset(self._ambient_ray_indices)
    … …  
    18521862          :class:`tuple` of :class:`cones <ConvexRationalPolyhedralCone>`;
    18531863
    18541864        - if neither ``dim`` nor ``codim`` is given, the output will be the
    1855           :class:`tuple` of tuples as above, giving faces of all dimension. If
    1856           you care about inclusion relations between faces, consider using
    1857           :meth:`face_lattice` or :meth:`adjacent`, :meth:`facet_of`, and
    1858           :meth:`facets`.
    1859 
     1865          :class:`tuple` of tuples as above, giving faces of all existing
     1866          dimensions. If you care about inclusion relations between faces,
     1867          consider using :meth:`face_lattice` or :meth:`adjacent`,
     1868          :meth:`facet_of`, and :meth:`facets`.
     1869         
    18601870        EXAMPLES:
    18611871
    18621872        Let's take a look at the faces of the first quadrant::
    … …  
    18851895            sage: quadrant.ray(0)
    18861896            N(1, 0)
    18871897           
     1898        Note that it is OK to ask for faces of too small or high dimension::
     1899       
     1900            sage: quadrant.faces(-1)
     1901            ()
     1902            sage: quadrant.faces(3)
     1903            ()
     1904           
     1905        In the case of non-strictly convex cones even faces of small
     1906        non-negative dimension may be missing::
     1907       
     1908            sage: halfplane = Cone([(1,0), (0,1), (-1,0)])
     1909            sage: halfplane.faces(0)
     1910            ()
     1911            sage: halfplane.faces()
     1912            ((1-d face of 2-d cone in 2-d lattice N,),
     1913             (2-d cone in 2-d lattice N,))
     1914            sage: plane = Cone([(1,0), (0,1), (-1,-1)])
     1915            sage: plane.faces(1)
     1916            ()
     1917            sage: plane.faces()
     1918            ((2-d cone in 2-d lattice N,),)
     1919           
    18881920        TESTS:
    18891921
    18901922        Now we check that "general" cones whose dimension is smaller than the
    … …  
    19181950            raise ValueError(
    19191951                    "dimension and codimension cannot be specified together!")
    19201952        dim = self.dim() - codim if codim is not None else dim
    1921         if dim is not None and (dim < 0 or dim > self.dim()):
    1922             raise ValueError("%s does not have faces of dimension %d!"
    1923                              % (self, dim))
    19241953        if "_faces" not in self.__dict__:
    19251954            self._faces = tuple(self._sort_faces(e.element for e in level)
    19261955                                for level in self.face_lattice().level_sets())
    19271956            # To avoid duplication and ensure order consistency
    1928             if self.dim() > 0:
     1957            if len(self._faces) > 1:
    19291958                self._facets = self._faces[-2]
    19301959        if dim is None:
    19311960            return self._faces
    19321961        else:
    1933             return self._faces[dim]
     1962            lsd = self.linear_subspace().dimension()
     1963            return self._faces[dim - lsd] if lsd <= dim <= self.dim() else ()
    19341964
    19351965    def facet_normals(self):
    19361966        r"""
    … …  
    22062236        if self._ambient is cone._ambient:
    22072237            # Cones are always faces of their ambient structure, so
    22082238            return self.ray_set().issubset(cone.ray_set())
    2209         # Cases of the origin and the whole cone (we don't have empty cones)
    2210         if self.nrays() == 0 or self.is_equivalent(cone):
     2239        if self.is_equivalent(cone):
    22112240            return True
    22122241        # Obviously False cases
    22132242        if (self.dim() >= cone.dim() # if == and face, we return True above
    … …  
    25182547        # Modify ray labels to match the ambient cone or fan.
    25192548        tp.ray_label = label_list(tp.ray_label, self.nrays(), deg <= 2,
    25202549                                   self.ambient_ray_indices())
    2521         result = tp.plot_lattice() + tp.plot_rays() + tp.plot_generators()
    2522         if self.dim() >= 2:
    2523             # Modify wall labels to match the ambient cone or fan too.
    2524             walls = self.faces(2)
    2525             try:
    2526                 ambient_walls = self.ambient().faces(2)
    2527             except AttributeError:
    2528                 ambient_walls = self.ambient().cones(2)
    2529             tp.wall_label = label_list(tp.wall_label, len(walls), deg <= 2,
    2530                                 [ambient_walls.index(wall) for wall in walls])
    2531             tp.set_rays(self.ambient().rays())
    2532             result += tp.plot_walls(walls)
     2550        result = tp.plot_lattice() + tp.plot_generators()
     2551        # To deal with non-strictly convex cones we separate rays and labels.
     2552        result += tp.plot_ray_labels()
     2553        tp.ray_label = None
     2554        lsd = self.linear_subspace().dimension()
     2555        if lsd == 1:
     2556            # Plot only rays of the line
     2557            v = self.lines()[0]
     2558            tp.set_rays([v, -v])
     2559        if lsd <= 1:
     2560            result += tp.plot_rays()       
     2561        # Modify wall labels to match the ambient cone or fan too.
     2562        walls = self.faces(2)
     2563        try:
     2564            ambient_walls = self.ambient().faces(2)
     2565        except AttributeError:
     2566            ambient_walls = self.ambient().cones(2)
     2567        tp.wall_label = label_list(tp.wall_label, len(walls), deg <= 2,
     2568                            [ambient_walls.index(wall) for wall in walls])
     2569        tp.set_rays(self.ambient().rays())
     2570        result += tp.plot_walls(walls)
    25332571        return result
    25342572
    25352573    def polyhedron(self):

  • sage/geometry/fan.py 

    diff -r 3e3e79e371f9 -r 2572cc23715e sage/geometry/fan.py
    a b  
    16301630
    16311631        - :class:`tuple` of cones of ``self`` of the specified (co)dimension,
    16321632          if either ``dim`` or ``codim`` is given. Otherwise :class:`tuple` of
    1633           such tuples for all dimensions.
     1633          such tuples for all existing dimensions.
    16341634
    16351635        EXAMPLES::
    16361636
    … …  
    16471647            N(-1, 0)
    16481648            sage: fan.cones(2)
    16491649            (2-d cone of Rational polyhedral fan in 2-d lattice N,)
     1650           
     1651        You cannot specify both dimension and codimension, even if they
     1652        "agree"::
     1653       
    16501654            sage: fan(dim=1, codim=1)
    16511655            Traceback (most recent call last):
    16521656            ...
    16531657            ValueError: dimension and codimension
    16541658            cannot be specified together!
     1659           
     1660        But it is OK to ask for cones of too high or low (co)dimension::
     1661       
     1662            sage: fan(-1)
     1663            ()
     1664            sage: fan(3)
     1665            ()
     1666            sage: fan(codim=4)
     1667            ()
    16551668        """
    16561669        if "_cones" not in self.__dict__:
    16571670            levels = [(e.element for e in level) # Generators
    … …  
    16741687                rays = list(levels[1])
    16751688                rays.sort(key=lambda cone: cone.ambient_ray_indices()[0])
    16761689                levels[1] = rays
    1677             self._cones = tuple(tuple(level) for level in levels)
    1678         if dim is None and codim is None:
    1679             return self._cones
    1680         elif dim is None:
    1681             return self._cones[self.dim() - codim]
    1682         elif codim is None:
    1683             return self._cones[dim]
    1684         raise ValueError(
     1690            self._cones = tuple(tuple(level) for level in levels)           
     1691        if dim is None:
     1692            if codim is None:
     1693                return self._cones
     1694            dim = self.dim() - codim
     1695        elif codim is not None:
     1696            raise ValueError(
    16851697                    "dimension and codimension cannot be specified together!")
     1698        return self._cones[dim] if 0 <= dim < len(self._cones) else ()
    16861699
    16871700    def contains(self, cone):
    16881701        r"""

  • sage/geometry/polyhedra.py 

    diff -r 3e3e79e371f9 -r 2572cc23715e sage/geometry/polyhedra.py
    a b  
    65706570    coatom_to_atoms = [frozenset(cta) for cta in coatom_to_atoms]
    65716571    C = frozenset(range(len(coatom_to_atoms)))  # All coatoms
    65726572    # Comments with numbers correspond to steps in Section 2.5 of the article
    6573     L = DiGraph()       # 3: initialize L
     6573    L = DiGraph(1)       # 3: initialize L
    65746574    faces = dict()
    65756575    atoms = frozenset()
    65766576    coatoms = C

  • sage/geometry/toric_plotter.py 

    diff -r 3e3e79e371f9 -r 2572cc23715e sage/geometry/toric_plotter.py
    a b  
    5656from sage.plot.all import (Color, Graphics,
    5757                           arrow, disk, line, point, polygon, rainbow, text)
    5858from sage.plot.plot3d.all import text3d
    59 from sage.rings.all import RDF, RR
     59from sage.rings.all import RDF
    6060from sage.structure.sage_object import SageObject
    6161
    6262
    … …  
    213213            raise ValueError("toric objects can be plotted only for "
    214214                             "dimensions 1, 2, and 3, not %s!" % dimension)
    215215        self.dimension = dimension
    216         self.origin = vector(RR, max(dimension, 2)) # 1-d is plotted in 2-d
     216        self.origin = vector(RDF, max(dimension, 2)) # 1-d is plotted in 2-d
    217217        if self.mode not in ["box", "generators", "round"]:
    218218            raise ValueError("unrecognized plotting mode: %s!" % mode)
     219        # If radius was explicitly set by the user, it sets other bounds too.
     220        # If we don't take it into account here, they will be replaced by
     221        # automatically computed values.
     222        if sd["radius"] is not None:
     223            for key in ["xmin", "ymin", "zmin"]:
     224                if sd[key] is None:
     225                    sd[key] = - sd["radius"]
     226            for key in ["xmax", "ymax", "zmax"]:
     227                if sd[key] is None:
     228                    sd[key] = sd["radius"]
    219229        # We also set some of the "extra_options" if they were not given.
    220230        if "axes" not in extra_options:
    221231            extra_options["axes"] = False
    … …  
    288298        sd = self.__dict__
    289299        bounds = ["radius", "xmin", "xmax", "ymin", "ymax", "zmin", "zmax"]
    290300        bounds = [abs(sd[bound]) for bound in bounds if sd[bound] is not None]
    291         r = RR(max(bounds + [0.5]) if bounds else 2.5)
     301        r = RDF(max(bounds + [0.5]) if bounds else 2.5)
    292302        self.radius = r
    293303        round = self.mode == "round"
    294304        for key in ["xmin", "ymin", "zmin"]:
    … …  
    296306                sd[key] = - r
    297307            if sd[key] > - 0.5:
    298308                sd[key] = - 0.5
    299             sd[key] = RR(sd[key])
     309            sd[key] = RDF(sd[key])
    300310        for key in ["xmax", "ymax", "zmax"]:
    301311            if round or sd[key] is None:
    302312                sd[key] = r
    303313            if sd[key] < 0.5:
    304314                sd[key] = 0.5
    305             sd[key] = RR(sd[key])
     315            sd[key] = RDF(sd[key])
    306316        if self.show_lattice is None:
    307317            self.show_lattice = (r <= 5) if d <= 2 else r <= 3
    308318       
    … …  
    340350        """
    341351        if not points:
    342352            return
    343         points = [vector(RR, pt) for pt in points]
     353        points = [vector(RDF, pt) for pt in points]
    344354        sd = self.__dict__
    345355       
    346356        def update(bound, new_value, points):
    … …  
    402412                if d <= 2:
    403413                    result += arrow(origin, generator,
    404414                                    color=color, width=thickness,
     415                                    arrowsize=thickness + 1,
    405416                                    zorder=zorder, **extra_options)
    406417                else:
    407418                    result += line([origin, generator], arrow_head=True,
    … …  
    510521        return point(points, color=self.point_color, size=self.point_size,
    511522                     zorder=self.point_zorder, **self.extra_options)
    512523       
     524    def plot_ray_labels(self):
     525        r"""
     526        Plot ray labels.
     527       
     528        Usually ray labels are plotted together with rays, but in some cases it
     529        is desirable to output them separately.
     530       
     531        Ray generators must be specified during construction or using
     532        :meth:`set_rays` before calling this method.
     533                 
     534        OUTPUT:
     535       
     536        - a plot.       
     537       
     538        EXAMPLES::
     539       
     540            sage: from sage.geometry.toric_plotter import ToricPlotter
     541            sage: tp = ToricPlotter(dict(), 2, [(3,4)])
     542            sage: print tp.plot_ray_labels()
     543            Graphics object consisting of 1 graphics primitive
     544        """
     545        return self.plot_labels(self.ray_label,
     546                                [1.1 * ray for ray in self.rays])
     547       
    513548    def plot_rays(self):
    514549        r"""
    515         Plot ray generators and their labels.
     550        Plot rays and their labels.
    516551       
    517552        Ray generators must be specified during construction or using
    518553        :meth:`set_rays` before calling this method.
    … …  
    541576            result += line([origin, end],
    542577                           color=color, thickness=thickness,
    543578                           zorder=zorder, **extra_options)
    544         result += self.plot_labels(self.ray_label, [1.1 * ray for ray in rays])
     579        result += self.plot_ray_labels()
    545580        return result
    546581       
    547582    def plot_walls(self, walls):
    … …  
    594629        elif mode == "generators":
    595630            origin = self.origin
    596631            for wall, color in zip(walls, colors):
     632                vertices = [rays[i] for i in wall.ambient_ray_indices()]
     633                vertices.append(origin)
     634                result += Polyhedron(vertices=vertices, field=RDF).render_solid(
     635                    alpha=alpha, color=color, zorder=zorder, **extra_options)
     636        label_sectors = []
     637        round = mode == "round"
     638        for wall, color in zip(walls, colors):
     639            S = wall.linear_subspace()
     640            lsd = S.dimension()
     641            if lsd == 0:    # Strictly convex wall
    597642                r1, r2 = (rays[i] for i in wall.ambient_ray_indices())
    598                 result += polygon([origin, r1, r2],
    599                     alpha=alpha, color=color, zorder=zorder, **extra_options)
    600         elif mode == "round":
    601             for wall, color in zip(walls, colors):
    602                 r1, r2 = (rays[i] for i in wall.ambient_ray_indices())
     643            elif lsd == 1:  # wall is a half-plane
     644                for i, ray in zip(wall.ambient_ray_indices(), wall.rays()):
     645                    if ray in S:
     646                        r1 = rays[i]
     647                    else:
     648                        r2 = rays[i]
     649                if round:
     650                    # Plot one "extra" sector
     651                    result += sector(- r1, r2,
     652                      alpha=alpha, color=color, zorder=zorder, **extra_options)
     653            else:           # wall is a plane
     654                r1, r2 = S.basis()
     655                r1 = vector(RDF, r1)
     656                r1 = r1 / r1.norm() * self.radius
     657                r2 = vector(RDF, r2)
     658                r2 = r2 / r2.norm() * self.radius
     659                if round:
     660                    # Plot three "extra" sectors
     661                    result += sector(r1, - r2,
     662                      alpha=alpha, color=color, zorder=zorder, **extra_options)
     663                    result += sector(- r1, r2,
     664                      alpha=alpha, color=color, zorder=zorder, **extra_options)
     665                    result += sector(- r1, - r2,
     666                      alpha=alpha, color=color, zorder=zorder, **extra_options)
     667            label_sectors.append([r1, r2])
     668            if round:
    603669                result += sector(r1, r2,
    604670                    alpha=alpha, color=color, zorder=zorder, **extra_options)
    605671        result += self.plot_labels(self.wall_label,
    606                           [sum(rays[i] for i in wall.ambient_ray_indices()) / 3
    607                            for wall in walls])
     672                    [sum(label_sector) / 3 for label_sector in label_sectors])
    608673        return result
    609674
    610675    def set_rays(self, generators):
    … …  
    637702        """
    638703        d = self.dimension
    639704        if d == 1:
    640             generators = [vector(RR, 2, (gen[0], 0)) for gen in generators]
     705            generators = [vector(RDF, 2, (gen[0], 0)) for gen in generators]
    641706        else:
    642             generators = [vector(RR, d, gen) for gen in generators]
     707            generators = [vector(RDF, d, gen) for gen in generators]
    643708        self.generators = generators
    644709        if self.mode == "box":
    645710            rays = []
    … …  
    10261091        sage: print sector((3,2,1), (1,2,3))
    10271092        Graphics3d Object
    10281093    """
    1029     ray1 = vector(RR, ray1)
    1030     ray2 = vector(RR, ray2)
     1094    ray1 = vector(RDF, ray1)
     1095    ray2 = vector(RDF, ray2)
    10311096    r = ray1.norm()
    10321097    if len(ray1) == 2:
    10331098        # Plot an honest sector
    … …  
    10451110        dr = (ray2 - ray1) / n
    10461111        points = (ray1 + i * dr for i in range(n + 1))
    10471112        points = [r / pt.norm() * pt for pt in points]
    1048         points.append(vector(RR, 3))
     1113        points.append(vector(RDF, 3))
    10491114        return polygon(points, **extra_options)