diff --git a/cutgeneratingfunctionology/igp/parametric.sage b/cutgeneratingfunctionology/igp/parametric.sage index 73b49bda0..06b7b71df 100644 --- a/cutgeneratingfunctionology/igp/parametric.sage +++ b/cutgeneratingfunctionology/igp/parametric.sage @@ -2250,7 +2250,7 @@ class SemialgebraicComplex(SageObject): else: self.add_new_component(var_value, bddbsa=self.bddbsa, flip_ineq_step=0, goto_lower_dim=False) - def plot(self, alpha=0.5, plot_points=300, slice_value=None, goto_lower_dim=False, restart=True, **kwds): + def plot(self, alpha=0.5, plot_points=300, slice_value=None, goto_lower_dim=False, restart=True, animated=False, **kwds): r""" Plot the complex and store the graph. @@ -2283,6 +2283,10 @@ class SemialgebraicComplex(SageObject): Plot the slice in (x,z)-space with y=x/2:: sage: complex.plot(slice_value=[u, u/2, v]) # not tested + + Animate the plot:: + + sage: complex.plot(slice_value=[None, None,4], animated=True) # not tested """ if restart: self.graph = Graphics() @@ -2308,49 +2312,84 @@ class SemialgebraicComplex(SageObject): ymax = kwds['ymax'] else: ymax = self.default_var_bound[1] - # # FIXME: zorder is broken in region_plot/ContourPlot. - # for c in self.components[self.num_plotted_components::]: - # num_eq = len(list(c.bsa.eq_poly())) - # gc = c.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, default_var_bound=self.default_var_bound, goto_lower_dim=goto_lower_dim, zorder=num_eq, **kwds) + # # zorder is broken in Region_Plot/ContourPlot. see issue #35992 in sagemath/sage. contorplot has regions with zorder=1 and lines with zorder=2. + # # code for when zorder is fixed. + # for c in self.components: + # num_eq = len(list((c.bsa.eq_poly()))) + # gc = c.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, default_var_bound=self.default_var_bound, goto_lower_dim=goto_lower_dim, zorder=num_eq, **kwds) # if gc: # need this because (empty g + empty gc) forgets about xmin xmax ymin ymax. # self.graph += gc - # Workaround. - # FIXME: bug example plot(goto_lower_dim=True) in doctest of SemialgebraicComplex - components = self.components[self.num_plotted_components::] - for c in components: - if not list(c.bsa.eq_poly()): - if not goto_lower_dim: - self.graph += c.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, default_var_bound=self.default_var_bound, goto_lower_dim=False, zorder=0, **kwds) - else: - color = find_region_color(c.region_type) - self.graph += (c.bsa).plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, color='white', fill_color=color, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) - if goto_lower_dim: - for c in components: - if not list(c.bsa.eq_poly()): - color = find_region_color(c.region_type) - for l in c.bsa.le_poly(): - new_bsa = copy(c.bsa) - new_bsa.add_polynomial_constraint(l, operator.eq) - self.graph += new_bsa.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, color=color, fill_color=color, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) - for c in components: - if len(list(c.bsa.eq_poly()))==1: - self.graph += c.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, default_var_bound=self.default_var_bound, goto_lower_dim=False, zorder=0, **kwds) - if goto_lower_dim: - for c in components: - if len(list(c.bsa.eq_poly()))==1: - color = find_region_color(c.region_type) - for l in c.bsa.lt_poly(): - new_bsa = BasicSemialgebraicSet_eq_lt_le_sets(eq=list(c.bsa.eq_poly())+[l], lt=[ll for ll in c.bsa.lt_poly() if ll != l], le=list(c.bsa.le_poly())) - self.graph += new_bsa.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, color='white', fill_color='white', xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) - for l in c.bsa.le_poly(): - new_bsa = copy(c.bsa) - new_bsa.add_polynomial_constraint(l, operator.eq) - self.graph += new_bsa.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, color=color, fill_color=color, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) - for c in components: - if len(list(c.bsa.eq_poly()))==2: - ptcolor = find_region_color(c.region_type) - self.graph += point(c.var_value, color = ptcolor, zorder=10) - self.num_plotted_components = len(self.components) + # # Work around is done by having semialgebraic complex plot things in order. + cells_by_dim = sorted(self.components, key=lambda cell: len(list((cell.bsa.eq_poly())))) ## Since Z order is defined for the complex by no. equations, sort by no. equations + grouped_cells = [] #group cells by no. eqns + for k, g in itertools.groupby(cells_by_dim, key=lambda cell: len(list((cell.bsa.eq_poly())))): + grouped_cells.append(list(g)) + if animated: + artists = [] + if not goto_lower_dim: #goto_lower_dim=True, plotting by no. equations mimics zorder request with lower no. equations first , then higest number last. The cell plotting method works fine for this. + for c in cells_by_dim: + num_eq = len(list((c.bsa.eq_poly()))) + gc = c.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, default_var_bound=self.default_var_bound, goto_lower_dim=goto_lower_dim, zorder=num_eq, **kwds) + if gc: # need this because (empty g + empty gc) forgets about xmin xmax ymin ymax. + self.graph += gc + if animated: + artists.append(self.graph) + + if animated: + return animate(artists) #use sage's build in animation from a list of graphics objects + else: #in the goto_lower_dim case=False, semialgebraic complex handels all plotting to mimic zorder done by the proof cells. + for cell_group in grouped_cells: #cell_groups are all cells with the same number of equations ordered least to greatest by number of equations + for c in cell_group: #plot components which should have zorder 3*no_eqns + color = kwds.pop('color', find_region_color(c.region_type)) + zorder = len(list((c.bsa.eq_poly()))) + if not list(c.bsa.eq_poly()): + g = (c.bsa).plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, color='white', fill_color=color, zorder=3*zorder, **kwds) + if g: + self.graph += g + if animated: + artists.append(self.graph) + else: + g = (c.bsa).plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, color=color, fill_color=color, zorder=3*zorder, **kwds) + if g: + self.graph += g + if animated: + artists.append(self.graph) + for c in cell_group: #plot components which should have zorder 3*no_eqns + 1 + zorder = len(list((c.bsa.eq_poly()))) + color = kwds.pop('color', find_region_color(c.region_type)) + if not list(c.bsa.eq_poly()): + pass + else: + for l in c.bsa.lt_poly(): # + new_bsa = BasicSemialgebraicSet_eq_lt_le_sets(eq=list(c.bsa.eq_poly())+[l], lt=[ll for ll in c.bsa.lt_poly() if ll != l], le=list(c.bsa.le_poly())) + g = new_bsa.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, color='white', fill_color='white', zorder=3*zorder+1, **kwds) + if g: + self.graph += g + if animated: + artists.append(self.graph) + for c in cell_group: #plot components which should have zorder 3*no_eqns + 2 + zorder = len(list((c.bsa.eq_poly()))) + color = kwds.pop('color', find_region_color(c.region_type)) + if not list(c.bsa.eq_poly()): + for l in c.bsa.le_poly(): + new_bsa = copy(c.bsa) + new_bsa.add_polynomial_constraint(l, operator.eq) + g = new_bsa.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, color=color, fill_color=color, zorder=3*zorder+2, **kwds) + if g: + self.graph += g + if animated: + artists.append(self.graph) + else: + for l in c.bsa.le_poly(): # fill in boundaries if possible it makes sense. + new_bsa = copy(c.bsa) + new_bsa.add_polynomial_constraint(l, operator.eq) + g = new_bsa.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, color=color, fill_color=color, zorder=3*zorder+2, **kwds) + if g: + self.graph += g + if animated: + artists.append(self.graph) + if animated: + return animate(artists, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) return self.graph def plot_bfs_tree(self, **kwds): diff --git a/cutgeneratingfunctionology/spam/basic_semialgebraic.py b/cutgeneratingfunctionology/spam/basic_semialgebraic.py index 9ba0ab697..5d6e26412 100644 --- a/cutgeneratingfunctionology/spam/basic_semialgebraic.py +++ b/cutgeneratingfunctionology/spam/basic_semialgebraic.py @@ -749,13 +749,13 @@ def plot(self, alpha=0.5, plot_points=300, slice_value=None, color='blue', fill_ section = self.section(slice_value) else: raise NotImplementedError("Plotting with dimension not equal to 2 is not implemented.") - return section.plot(alpha=alpha, plot_points=plot_points, slice_value=None, + return BasicSemialgebraicSet_veronese.from_bsa(section).plot(alpha=alpha, plot_points=plot_points, slice_value=None, color=color, fill_color=fill_color, constraints_color=constraints_color, constraints_fill_color=constraints_fill_color, eq_constraints_kwds=eq_constraints_kwds, le_constraints_kwds=le_constraints_kwds, - lt_constraints_kwds=lt_constraints_kwds, **kwds) + lt_constraints_kwds=lt_constraints_kwds, **kwds) # use a BSA that can prove emptiness and reduce inequalities g = Graphics() for bv in 'xmin', 'xmax', 'ymin', 'ymax': b = kwds.get(bv, None) @@ -783,33 +783,35 @@ def plot(self, alpha=0.5, plot_points=300, slice_value=None, color='blue', fill_ x, y = SR.var('x, y') x_range = (x, xmin-0.01, xmax+0.01) y_range = (y, ymin-0.01, ymax+0.01) - eq_constraints = [ l(x, y) == 0 for l in self.eq_poly() ] + eq_constraints = [ l(x, y) == 0 for l in self.eq_poly() ] lt_constraints = [ l(x, y) < 0 for l in self.lt_poly() ] le_constraints = [ l(x, y) <= 0 for l in self.le_poly() ] + # note: region_plot does not plot inconsistent inequalites + # TO DO: Dealing with plotting R^2 as represetned as empty list in eq, lt, le constraints. + # TO DO: Special case for dealing when the 2d constraints that reduce to a point. For now region_plot can handle it with a warning. + # note region_plot has issues with zorder keyword see issue #35992 in sagemath/sage if constraints_fill_color: for c in chain(lt_constraints, le_constraints): if not isinstance(c, bool): g += region_plot([c], x_range, y_range, - alpha=0.1, plot_points=plot_points, - incol=constraints_fill_color, bordercol=None, **kwds) + alpha=0.1, plot_points=plot_points, + incol=constraints_fill_color, bordercol=None, **kwds) if color or fill_color: - all_constraints = eq_constraints + lt_constraints + le_constraints - non_trivial_constraints = [ c for c in all_constraints if not isinstance(c, bool) ] - if not any(c is False for c in all_constraints): - g += region_plot(non_trivial_constraints, x_range, y_range, - alpha=alpha, plot_points=plot_points, - incol=fill_color, bordercol=color, **kwds) + all_constraints = eq_constraints + lt_constraints + le_constraints + g += region_plot(all_constraints, x_range, y_range, + alpha=alpha, plot_points=plot_points, + incol=fill_color, bordercol=color, **kwds) if constraints_color: for polys, plot_kwds in [(self.eq_poly(), eq_constraints_kwds), - (self.lt_poly(), lt_constraints_kwds), - (self.le_poly(), le_constraints_kwds)]: + (self.lt_poly(), lt_constraints_kwds), + (self.le_poly(), le_constraints_kwds)]: for l in polys: c = l(x, y) == 0 if not isinstance(c, bool): effective_kwds = copy(kwds) effective_kwds['color'] = constraints_color effective_kwds.update(plot_kwds) - g += implicit_plot(c, x_range, y_range, **effective_kwds) + g += implicit_plot(c, x_range, y_range, **effective_kwds) # implicit_plot has no zorder information. return g class BasicSemialgebraicSet_polyhedral(BasicSemialgebraicSet_base): @@ -1015,13 +1017,14 @@ def eq_poly(self): Together, ``eq_poly``, ``lt_poly``, and ``le_poly`` describe ``self``. """ + ## add tests for c in self._polyhedron.minimized_constraints(): if c.is_equality(): coeff = c.coefficients() # observe: coeffients in a constraint of NNC_Polyhedron could have gcd != 1. # take care of this. gcd_c = gcd(gcd(coeff), c.inhomogeneous_term()) - t = sum(QQ(x)/gcd_c*y for x, y in zip(coeff, self.poly_ring().gens())) + QQ(c.inhomogeneous_term())/gcd_c - yield t + t = sum(QQ(x)/gcd_c*y for x, y in zip(coeff, self.poly_ring().gens())) + QQ(c.inhomogeneous_term())/gcd_c # not type stable, make it type stable + yield self.poly_ring()(t) def lt_poly(self): r""" @@ -1036,7 +1039,7 @@ def lt_poly(self): gcd_c = gcd(gcd(coeff), c.inhomogeneous_term()) # constraint is written with '>', while lt_poly records '<' relation t = sum(-QQ(x)/gcd_c*y for x, y in zip(coeff, self.poly_ring().gens())) - QQ(c.inhomogeneous_term())/gcd_c - yield t + yield self.poly_ring()(t) def le_poly(self): r""" @@ -1051,7 +1054,7 @@ def le_poly(self): gcd_c = gcd(gcd(coeff), c.inhomogeneous_term()) # constraint is written with '>=', while lt_poly records '<=' relation t = sum(-QQ(x)/gcd_c*y for x, y in zip(coeff, self.poly_ring().gens())) - QQ(c.inhomogeneous_term())/gcd_c - yield t + yield self.poly_ring()(t) # override the default implementation def __contains__(self, point): @@ -1733,7 +1736,7 @@ def __init__(self, upstairs_bsa, polynomial_map, poly_ring=None, ambient_dim=Non super(BasicSemialgebraicSet_section, self).__init__(base_ring=base_ring, ambient_dim=ambient_dim, names=names, poly_ring=poly_ring) self._upstairs_bsa = upstairs_bsa self._polynomial_map = polynomial_map - + def __copy__(self): r""" Make a copy of ``self``. @@ -1744,15 +1747,15 @@ def __copy__(self): def eq_poly(self): for p in self._upstairs_bsa.eq_poly(): - yield p(self._polynomial_map) + yield self.poly_ring()(p(self._polynomial_map)) # these need to be type stable def le_poly(self): for p in self._upstairs_bsa.le_poly(): - yield p(self._polynomial_map) + yield self.poly_ring()(p(self._polynomial_map)) def lt_poly(self): for p in self._upstairs_bsa.lt_poly(): - yield p(self._polynomial_map) + yield self.poly_ring()(p(self._polynomial_map)) def _repr_(self): return 'BasicSemialgebraicSet_section({}, polynomial_map={})'.format(self._upstairs_bsa, self._polynomial_map) @@ -1939,7 +1942,7 @@ def __init__(self, upstairs_bsa=None, polynomial_map=None, v_dict=None, self._v_dict = v_dict @classmethod - def from_bsa(cls, bsa, poly_ring=None, **init_kwds): + def from_bsa(cls, bsa, poly_ring=None, **init_kwds): #add example if poly_ring is None: poly_ring = bsa.poly_ring() return super(BasicSemialgebraicSet_veronese, cls).from_bsa(bsa, poly_ring=poly_ring, **init_kwds)