From 14d044f61682f34790e9ff32eda809126dba732c Mon Sep 17 00:00:00 2001 From: Acadia Larsen <102884863+ComboProblem@users.noreply.github.com> Date: Wed, 6 Sep 2023 08:09:16 -0700 Subject: [PATCH 1/2] Update parametric.sage fix plotting method in semialgebraic complex with concern to broken z order of parametric plot and region plot. added animated feature to semialgebraic complex plot --- .../igp/parametric.sage | 123 ++++++++++++------ 1 file changed, 81 insertions(+), 42 deletions(-) diff --git a/cutgeneratingfunctionology/igp/parametric.sage b/cutgeneratingfunctionology/igp/parametric.sage index 73b49bda0..d193a62c5 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) return self.graph def plot_bfs_tree(self, **kwds): From 1bdb42573a95a36163e5f7128e30081a2ad10f47 Mon Sep 17 00:00:00 2001 From: Acadia Larsen <102884863+ComboProblem@users.noreply.github.com> Date: Wed, 6 Sep 2023 08:13:22 -0700 Subject: [PATCH 2/2] Update basic_semialgebraic.py update plot method to use bsa that can prove emptiness fixing a plotting bug where semialgebraic complex would try to plot an empty set update NCC polyhedron to fix type error that happened when plot method was trying to plot an empty set. --- .../igp/parametric.sage | 6 +-- .../spam/basic_semialgebraic.py | 49 ++++++++++--------- 2 files changed, 29 insertions(+), 26 deletions(-) diff --git a/cutgeneratingfunctionology/igp/parametric.sage b/cutgeneratingfunctionology/igp/parametric.sage index d193a62c5..06b7b71df 100644 --- a/cutgeneratingfunctionology/igp/parametric.sage +++ b/cutgeneratingfunctionology/igp/parametric.sage @@ -2342,7 +2342,7 @@ class SemialgebraicComplex(SageObject): 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()): + 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 @@ -2387,9 +2387,9 @@ class SemialgebraicComplex(SageObject): if g: self.graph += g if animated: - artists.append(self.graph) + artists.append(self.graph) if animated: - return animate(artists) + 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)