-
Notifications
You must be signed in to change notification settings - Fork 0
GSoC 2015 Application Luv Agarwal: Cylindrical Algebraic Decomposition
Hello, I am Luv Agarwal a sophomore at IIIT - Hyderabad, India pursuing a degree in B-Tech in Computer Science with MS by research in Computational and Natural Sciences. I work on ubuntu 14.04 and use vim and sublime interchangeably as my text editor. I have done coding in a low level language ARM and various high level languages (C, C++, PHP, Python, javascript). I have been using git (with bitbucket) for around 1.5 years, github and mercurial for around a year and quite familiar with their workflow.
Name : Luv Agarwal
University : International institute of Information Technology - Hyderabad
email : agarwal.iiit@gmail.com
github username : luviiit
merged https://github.com/sympy/sympy/pull/8931
unmerged https://github.com/sympy/sympy/pull/8858
Some terms used later:
- CAD - Cylindrical Algebraic Decomposition
- Projection | Projection operator - A component of CAD algorithm.
Solving real multivariable polynomial systems is a crucial feature for a Computer Algebra System in a long term. Currently, we don't have support for handling such systems.
There are several algorithms that are used to deal with such systems.
CAD algorithm is the best algorithm known till date to deal with real polynomial systems in general. In time, there are different modifications done in CAD in order to improve it's performance to solve different problems or workout some special cases. Some of these enhancements are as follows:
- Hong's improved Projection operator (This operator is complete and is used in place of the original projection operator introduced by Collin)
- George E. Collins and Hoon Hong's Partial Cylindrical Algebraic Decomposition for Quantifier Elimination.
- Adam W. Strzebonski Cylindrical Algebraic Decomposition using validated numerics - used to reduce the algebraic number computations
As Wolfram states "As The primary method used by the Wolfram Language for solving real polynomial systems and real quantifier elimination is the CAD algorithm. There are, however, simpler methods applicable in special cases."
In Wolfram language CAD implementation uses the Collins-Hong algorithm (Partial CAD as mentioned above) with Brown-McCallum projection operator for a special condition - if set of input polynomials is a well-oriented set and Hong projection operator in general. While contruncting CAD, we need to deal with plenty of algebraic number computations (computationally intensive) which can sometimes be omitted by using Strzebonski's genealogy-based method using validated numerics backed up by exact algebraic number.
The project comprise of implementation of following things:
- Cylindrical Algebraic Decomposition with Hong's improved projection operator.
- Hong's partial CAD derived from complete CAD^ to perform quantifier elimination.
- A platform to solve real polynomial inequalities using CAD.
Given n polynomials in r variables, the objective of CAD algorithm is to divide complete r dimensional space in cells such that each of the n polynomials is invariant on each cell. In other words each input polynomial has the same sign over complete cell and therefore each cell can be denoted by a sample point which can be used to calculate the sign of each polynomial over every point in that cell.
Now I will now explain the algorithms in some more detail in following way:
1. CAD in two parts:
1.1 a general CAD algorithm assuming all algebraic numbers are rational numbers so that I can represent all the numbers exactly in numerical form and don't have to deal with algebraic number manipulations which will help focus more on understanding the principles of the algorithm.
1.2 workout for removing the assumption in 1.1
2. partial CAD using complete CAD^ to perform quantifier elimination.
3. Inequalities solving
Input: a set of polynomials in r variables.
Output: r-dimensional cells each represented by sample point and definition
1. compute the coarsest squarefree basis of input polynomials using basis (to make sure that each polynomial doesn't have a repeated root)
2. Apply projection operator iteratively r-1 times so that we end up with projection polynomials (proj(A), proj^2(A), ... proj^(r-1)(A)) where proj^(k+1)(A) = proj(proj^k(A)).
3. Start from proj^(r-1)(A).
Let the roots of these polynomials be a[1], a[2], ..., a[n]
Set cells = [(b[1], ), (a[1], ), (b[2], ), (a[2], ), ..., (a[n], ), (b[n+1], )]
where b[i] is any rational number between a[i-1] and a[i] for i in {2..n}
b[1] is any rational number less than a[1]
b[n+1] is any rational number greater than a[n]
Now each cell is a one-dimensional point.
4. Now define each cell using define_base.
5. Extension phase:
level = 1
for p in reverse(projection_polynomials):
level += 1
A = p
cells_new = []
# cells = [(b[1][1], .., b[1][level]), (b[2][1], .., b[2][level]), ..., (b[m][1], .., b[m][level])]
for cell in cells:
A = A(cell) # substitute the cell into A
# Now A is polynomial in variable x[level]
A = basis(A)
r = roots(A)
cells_new.extend([cell.append(root) for root in r])
cells = cells_new
now define all these new cells using define
Now we have the cells and their definitions.
1. compute the coarsest squarefree basis of input polynomials using basis(to make sure that each polynomial doesn't have a repeated root)
2. Apply projection operator iteratively r-1 times so that we end up with projection polynomials (proj(A), proj^2(A), ... proj^(r-1)(A)) where proj^(k+1)(A) = proj(proj^k(A)).
3. Start from proj^(r-1)(A).
Let the roots of these polynomials be a[1], a[2], ..., a[n]
+ if not root.is_rational: represent it with polynomial over `Q` and isolating interval.
Set cells = [(b[1], ), (a[1], ), (b[2], ), (a[2], ), ..., (a[n], ), (b[n+1], )]
where b[i] is any rational number between a[i-1] and a[i] for i in {2..n}
b[1] is any rational number less than a[1]
b[n+1] is any rational number greater than a[n]
Now each cell is a one-dimensional point.
4. Now define each cell using define_base.
5. Extension phase:
level = 1
for p in reverse(projection_polynomials):
level += 1
A = p
cells_new = []
# cells = [(b[1][1], .., b[1][level]), (b[2][1], .., b[2][level]), ..., (b[m][1], .., b[m][level])]
for cell in cells:
+ # cell = (b<sub>i,1</sub>, .., b<sub>i,level</sub>)
+ alpha, B[1], B[2], ..., B[level] = simple(cell)
+ # B's are the polynomials over Q
+ # alpha is an algebraic number represented by an polynomial over Q and an isolating interval such that B[j](alpha) = b[i][j]
A = A(cell) # substitute the cell into A
# Now A is polynomial in variable x[level]
+ # Now A is a polynomial over Q(alpha)
A = basis(A)
r = roots(A)
cells_new.extend([cell.append(root) for root in r])
cells = cells_new
now define all these new cells using define
Now we have the cells and their definitions.
ii) Partial CAD to perform quantifier elimination:
a) After the construction of children of the cell try to determine the truth value of children
Function: `evaltv`
b) Try to determine the truth value of cell if a) is successful. If it is determined then try to determine the truth value of as many ancestors as possible dropping the subtree of the cell whose truth value is determined.
Function: `prptv`
b) Intead of simply running `bfs` (as we did in extension phase of above algorithm) we will choose cells such that it will result in `dropping of more subtrees`. For example we can choose cells based on their level (extension level), type of cell (sector or section), et al.
Function: `choose`
TODO
- cylindrical algebraic decomposition
>>> cylindrically_decompose(x**2 + y**2 < 1, {x, y})
And(-1 < x, x < 1, -sqrt(-x**2 + 1) < y, y < sqrt(x**2 - 1))
>>> cylindrically_decompose(x**2 + y**2 > 1, {x, y})
Or(And(Or(y < -sqrt(-x**2 + 1), y > sqrt(-x**2 + 1)), x <= 1, x >= -1), x < -1, x > 1)
TODO