Skip to content

Semialgebraic complex plot and other plotting methods #83

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ComboProblem wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Semialgebraic complex plot and other plotting methods #83

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
14d044f
Update parametric.sage
ComboProblem Sep 6, 2023
1bdb425
Update basic_semialgebraic.py
ComboProblem Sep 6, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
123 changes: 81 additions & 42 deletions cutgeneratingfunctionology/igp/parametric.sage
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -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()
Expand All @@ -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):
Expand Down
49 changes: 26 additions & 23 deletions cutgeneratingfunctionology/spam/basic_semialgebraic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -1015,13 +1017,14 @@ def eq_poly(self):

Together, ``eq_poly``, ``lt_poly``, and ``le_poly`` describe ``self``.
"""
## add tests
Copy link
Owner

@mkoeppe mkoeppe Sep 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why was this change (test for is_empty) necessary again? Is something wrong with the minimized_constraints?

Copy link
Contributor Author

@ComboProblem ComboProblem Sep 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't remember exactly. I think it was there as a fix from before we knew that BasicSemialgebraicSet_section‎ wasn't stable and when we had fixed that. I've done a few quick tests and it seems like the test for emptiness is not longer necessary.

Copy link
Owner

@mkoeppe mkoeppe Sep 15, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's back out this change then

Copy link
Contributor Author

@ComboProblem ComboProblem Sep 19, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've backed out those changes.

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"""
Expand All @@ -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"""
Expand All @@ -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):
Expand Down Expand Up @@ -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

Copy link
Owner

@mkoeppe mkoeppe Sep 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

keep this empty line

Copy link
Contributor Author

@ComboProblem ComboProblem Sep 19, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changed.

def __copy__(self):
r"""
Make a copy of ``self``.
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down