Skip to content

Commit

Permalink
...
Browse files Browse the repository at this point in the history
  • Loading branch information
dddejan committed Nov 16, 2015
1 parent 6ac07fb commit aad7439
Show file tree
Hide file tree
Showing 2 changed files with 51 additions and 70 deletions.
95 changes: 40 additions & 55 deletions examples/cad/cad.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,53 +41,38 @@ def cylinder_notify(self, cylinder, assignment):
print "Cylinder:\n", cylinder
pass


# Get all the reductums (including the polynomial itself), but not the constant
def get_reductums(f, x):
R = []
while f.var() == x:
R.append(f)
f = f.reductum()
return R

# Polynomials with sign conditions are considered a conjunction
CONJUNCTIVE = 0
# Polynomials with sign conditions are considered a disjunction
DISJUNCTIVE = 1

# Do cylindrical algebraic decomposition (CAD)
class CAD:

# Map from variables to polynomials
projection_map = {}

# Signs of relevant polynomials, mapping top variables to polynomials (for lifting)
sign_condition_map = {}

# Set of variables we're working with
variables = []

# Notification
cylinder_notify = CylinderNotify()

# Polynomials with sign conditions are considered a conjunction
CONJUCTIVE = 0
# Polynomials with sign conditions are considered a disjunction
DISJUNCTIVE = 1

# Type of CAD
type = CONJUCTIVE


# Initialize
def __init__(self, variable_list):
# Set the variable order
polypy.variable_order.set(variable_list)
# Map from variables to sets of polynomials
def __init__(self, variable_list, type = CONJUNCTIVE):
# Init variables
self.projection_map = {}
self.sign_condition_map = {}
self.variables = variable_list
self.cylinder_notify = CylinderNotify()
self.type = type
# Initialize the maps
for x in variable_list:
self.projection_map[x] = set()
self.sign_condition_map[x] = set()

# Variables
self.variables = variable_list

# Get all the reductums (including the polynomial itself), but not the constant
@staticmethod
def get_reductums(f, x):
R = []
while f.var() == x:
R.append(f)
f = f.reductum()
return R

# Set the variable order
polypy.variable_order.set(variable_list)

# Add a polynomial to CAD
def add_polynomial(self, f, sign_condition = None):
# Factor the polynomial
Expand All @@ -104,7 +89,7 @@ def add_polynomial(self, f, sign_condition = None):

# Check if the assignment respects the polynomials with top variable x
def check_assignment(self, x, assignment):
if self.type == CAD.CONJUCTIVE:
if self.type == CONJUNCTIVE:
for p, sgn_condition in self.sign_condition_map[x]:
if not p.sgn_check(assignment, sgn_condition):
return False
Expand All @@ -121,29 +106,26 @@ def add_polynomials(self, polynomials):
self.add_polynomial(f)

# Project. Go down the variable stack and project:
# [1] coeff(f) for f in poly[x]
# [2] psc(g, g') for f in poly[x], g in R(f, x)
# [3] psc(g1, g2) for f1, f2 in poly[x], g1 in R(f1, x), g2 in R(f2, x)
def project(self):
for x in reversed(self.variables):
# Project variable x
x_poly_set = self.projection_map[x]
# [1] coeff(f) for f in poly[x]
for f in self.projection_map[x]:
for f in x_poly_set:
self.add_polynomials(f.coefficients())
# [2] psc(g, g') for f in poly[x], g in R(f, x)
for f in self.projection_map[x]:
for g in CAD.get_reductums(f, x):
for f in x_poly_set:
for g in get_reductums(f, x):
g_d = f.derivative()
if (g_d.var() == x):
self.add_polynomials(g.psc(g_d))
# [3] psc(g1, g2) for f1, f2 in poly[x], g1 in R(f1, x), g2 in R(f2, x)
for (f1, f2) in itertools.combinations(self.projection_map[x], 2):
f1_R = CAD.get_reductums(f1, x)
f2_R = CAD.get_reductums(f2, x)
for (f1, f2) in itertools.combinations(x_poly_set, 2):
f1_R = get_reductums(f1, x)
f2_R = get_reductums(f2, x)
for (g1, g2) in itertools.product(f1_R, f2_R):
self.add_polynomials(g1.psc(g2))

self.add_polynomial(g1.psc(g2))
# Lift the first variable, update the assignment and lift recursively
def lift_first_var(self, variables, assignment, cylinder):
# We've tried all variables, this assignment checks out
Expand Down Expand Up @@ -190,7 +172,12 @@ def lift(self):
assignment = polypy.Assignment()
cylinder = Cylinder()
self.lift_first_var(self.variables, assignment, cylinder)


# Run the CAD construction
def run(self):
self.project()
self.lift()

# Print internal state
def print_state(self):
print "Variables:", self.variables
Expand All @@ -208,7 +195,5 @@ def print_state(self):
# Setup CAD
cad = CAD([x, y])
cad.add_polynomial(x**2 + y**2 - 1, polypy.SGN_GE_0)
# Project
cad.project()
# Lift
cad.lift()
# Run CAD
cad.run()
26 changes: 11 additions & 15 deletions examples/cad/plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,37 +6,33 @@
import matplotlib.pyplot as plt

# 2D plotting of polynomials
class PolyPlot2D(CylinderNotify):

# The variables
x = None
y = None

# List of polynomials with sign conditions
polynomials = None

class PolyPlot2D(cad.CylinderNotify):

# Initialize
def __init__(self, x, y):
self.x = x
self.y = y
self.cad = cad.CAD([x, y])
self.polynomials = []
self.cylinders = []

# Add a polynomial
def add_polynomial(self, f, sign_condition):
self.polynomials.append((f, sign_condition))
self.cad.add_polynomial(f, sign_condition)

# Show the plot
def show(self):
# Do CAD
# Run the CAD and collect cyllinders
mycad = cad.CAD([x, y])
mycad.project()
mycad.lift()
mycad.cylinder_notify = self
for (f, sign_condition) in self.polynomials:
mycad.add_polynomial(f, sign_condition)
mycad.run()

# Notifications on sylinders
def cylinder_notify(self, cylinder):
pass
def cylinder_notify(self, cylinder, assignment):
self.cylinders.append(cylinder)
print "Cylinder:\n", cylinder

if __name__ == "__main__":
# Some variables
Expand Down

0 comments on commit aad7439

Please sign in to comment.