Integration of polynomials over 2D polygons #12673

ArifAhmed1995 · 2017-05-25T19:49:18Z

Module for integrating uni/bivariate polynomials over 2-polytopes.

certik · 2017-05-31T13:42:53Z

Hi @ArifAhmed1995. Thanks for starting to work on this. The Travis failure seems unrelated (#12697).

Would you mind creating some jupyter notebook with examples (and plots) how to use this?

When you create new PRs, or just want my feedback, please tag me using @certik, then I'll get an email and can comment on your work.

Let me know if you have any questions.

ArifAhmed1995 · 2017-05-31T13:50:31Z

I'm close to getting the 2D prototype working. Once that's done I'll create the notebook and add it to the examples folder.

certik · 2017-05-31T13:52:56Z

Perfect!

certik

Also a test fails -- it could be perhaps caused by the 1/2 thing, per my other comment.

certik · 2017-06-03T00:44:53Z

sympy/integrals/intpoly.py

+        s = 0
+        for i in point.values():
+            s += i ** 2
+        return s**(1/2)


This is a common mistake, see http://docs.sympy.org/dev/gotchas.html#python-numbers-vs-sympy-numbers, you have to use S(1)/2.

certik · 2017-06-03T00:45:48Z

sympy/integrals/intpoly.py

+                                    intersect[1] - x0[gens[1]]))
+            if is_vertex(intersect):
+                if isinstance(expr, Expr):
+                    value += distance_origin *\


space before \

ArifAhmed1995 · 2017-06-07T13:30:58Z

@certik The tests should pass now. I've made the notebooks and have discussed the existing limitations of the methods therein.
We should discuss the final 2D API now so that I can implement it in this week (before 11th).

certik · 2017-06-08T17:44:28Z

@ndotsu, here is @ArifAhmed1995 latest code.

gxyd · 2017-06-09T08:41:03Z

sympy/integrals/intpoly.py

+    xl, yl = list(), list()
+
+    for vertex in poly.vertices:
+        xl.append(vertex.x)


It would be more efficient if you use xl = list(map((lambda vertex: vertex.x), poly.vertices)), similarly for yl.

gxyd · 2017-06-09T08:52:27Z

sympy/integrals/intpoly.py

+    """
+    if dims is None:
+        if isinstance(expr, Expr):
+            dims = list(expr.free_symbols)


Can we have tuple(expr.free_symbols)? I don't think you are using mutability of dims.

gxyd · 2017-06-09T08:53:03Z

sympy/integrals/intpoly.py

+        if isinstance(expr, Expr):
+            dims = list(expr.free_symbols)
+        else:
+            dims = [x, y]


dims = tuple([x, y])

gxyd · 2017-06-09T09:03:37Z

I would say, you should use var = S.Zero instead of python int assignment var = 0. Unless you intentionally don't want to use python int.

gxyd · 2017-06-09T09:07:46Z

sympy/integrals/intpoly.py

+    return ()
+
+
+def is_vertex(ent):


This sounds like it should be some class-level variable is_Vertex or is_vertex.

But how to do it for say (2, 3) or (1, 2, 3) ? Plus , this function may be expanded in the future.

Ah, leave it then. May be you can have it for any sequence (why just tuple?).

gxyd · 2017-06-09T09:11:04Z

sympy/integrals/intpoly.py

+    expr: Denotes a polynomial(SymPy expression)
+    """
+    from sympy.plotting.plot import plot3d, plot
+    gens = list(expr.free_symbols)


no need for a list, can work with gens = expr.free_symbols.

gxyd · 2017-06-09T09:17:19Z

sympy/integrals/intpoly.py

+        else:
+            poly_dict[1] = expr
+    else:
+        poly_dict = {0: expr}


better do poly_dict[0] = expr, instead of a complete variable re-assignment.

ArifAhmed1995 · 2017-06-09T11:18:34Z

@gxyd Thanks for the review ! I've made all but one of the changes.

gxyd · 2017-06-12T11:57:43Z

sympy/integrals/intpoly.py

        if isinstance(expr, Expr):
            dims = tuple(expr.free_symbols)
        else:
            dims = (x, y)

+    if kwargs['clockwise'] is True:


You can use

clockwise = kwargs.get('clockwise', True) if clockwise: poly = ...

Ah yes. I've made the change locally. There's a deficiency in the clockwise sort so I'll push the total changes after correcting that.

gxyd · 2017-06-12T12:01:58Z

sympy/integrals/intpoly.py

+        second = (b.x - center.x) * (b.x - center.x) + (b.y - center.y) * (b.y - center.y)
+        return -1 if first > second else 1
+
+    return Polygon(*sorted(vertices, key = cmp_to_key(compareTo)))


Looks better if you have key=cmp_to_key(...) i.e no space. https://www.python.org/dev/peps/pep-0008/#other-recommendations

certik

In general, please do not use floating point numbers, but exact numbers, I commented on the places where I noticed it.

certik · 2017-06-14T05:17:16Z

sympy/integrals/tests/test_intpoly.py

+    assert polytope_integrate(Polygon(Point(0, 0), Point(0, 1),
+                                      Point(1, 1), Point(1, 0)), x * y) == 1/4
+    assert polytope_integrate(Polygon(Point(0, 3), Point(5, 3), Point(1, 1)),
+                              6*x**2 - 40*y) == float(-935)/3


Why is this a floating point number? We should try to compare to an exact number, i.e., -S(935)/3.

certik · 2017-06-14T05:17:32Z

sympy/integrals/tests/test_intpoly.py

+    assert polytope_integrate(Polygon(Point(0, 0), Point(0, 2),
+                                      Point(4, 0)), 1, dims=(x, y)) == 4
+    assert polytope_integrate(Polygon(Point(0, 0), Point(0, 1),
+                                      Point(1, 1), Point(1, 0)), x * y) == 1/4


1/4 -> S(1)/4

certik · 2017-06-14T05:18:06Z

sympy/integrals/tests/test_intpoly.py

+    assert best_origin((2, 1), 3, l1, expr1) == (0, 3.0)
+    assert best_origin((2, 0), 3, l2, x ** 7) == (1.5, 0)
+    assert best_origin((0, 2), 3, l3, x ** 7) == (0, 1.5)
+    assert best_origin((1, 1), 2, l4, x ** 7 * y ** 3) == (0, 2.0)


Should return an exact number, not a floating point.

certik · 2017-06-14T05:18:14Z

sympy/integrals/tests/test_intpoly.py

+
+    l1 = Segment2D(Point(0, 3), Point(1, 1))
+    l2 = Segment2D(Point(1.5, 0), Point(1.5, 3))
+    l3 = Segment2D(Point(0, 1.5), Point(3, 1.5))


1.5 -> S(3)/2

certik · 2017-06-14T05:24:20Z

Also there seems to be a test failure at Travis:

File "/home/travis/virtualenv/python2.7.9/lib/python2.7/site-packages/sympy/integrals/intpoly.py", line 30, in sympy.integrals.intpoly.polytope_integrate
Failed example:
    polytope_integrate(poly, expr)
Exception raised:
    Traceback (most recent call last):
      File "/opt/python/2.7.9/lib/python2.7/doctest.py", line 1315, in __run
        compileflags, 1) in test.globs
      File "<doctest sympy.integrals.intpoly.polytope_integrate[6]>", line 1, in <module>
        polytope_integrate(poly, expr)
      File "/home/travis/virtualenv/python2.7.9/lib/python2.7/site-packages/sympy/integrals/intpoly.py", line 33, in polytope_integrate
        if kwargs['dims'] is None:
    KeyError: 'dims'

ArifAhmed1995 · 2017-06-18T10:33:48Z

@certik I updated the code. The failing test runs fine locally :

============================= test process starts ==============================
executable:         /home/arif7/anaconda3/bin/python  (3.6.0-final-0) [CPython]
architecture:       64-bit
cache:              yes
ground types:       python 
random seed:        16576427
hash randomization: on (PYTHONHASHSEED=1314471289)

sympy/assumptions/tests/test_refine.py[8] ........                          [OK]

________________________________ slowest tests _________________________________
test_pow - Took 180.055 seconds
================= tests finished: 8 passed, in 184.17 seconds ==================

ArifAhmed1995 · 2017-06-18T12:21:23Z

@certik As mentioned in my blog post there are two major issues :

The current implementation of Polygons in the geometry module does not account for polygons with intersecting sides. I got the following error when running tests for polytopes with intersecting sides :

sympy.geometry.exceptions.GeometryError: Polygon has intersecting sides.

Should I go ahead and modify the existing polygon class to accommodate this class of polygons ?

The computation of best origin seems to be expensive. I'll try to optimize it by including some inexpensive preliminary checks and simply return first vertex if difference between exponents of variables is not too much. What that difference should be has to be decided by plotting some graphs and then set value to be that intersection point.

certik · 2017-06-27T13:45:19Z

@ArifAhmed1995 I thought @ndotsu already commented on this --- but I don't see his answer now. Let me find it.

Regarding the remaining Travis failure (https://travis-ci.org/sympy/sympy/jobs/244146796), that is not related to this PR.

ArifAhmed1995 · 2017-06-28T22:55:48Z

@certik I've added the dynamic programming algorithm. This is a scalable technique and much faster than the standard one(200 - 300 times faster for higher power polynomials). I'll work on the intersecting polygon case now.
At the end of this week, I'll explain all the changes in a blog post. Maybe I can still fine-tune the new technique.

certik · 2017-07-08T14:19:49Z

Thanks @gxyd for the review and @ArifAhmed1995 for fixing it. If tests pass, this is now ok to merge with me.

I have one last request. The commit d76be087beba0e974a72db92d4478a2ce9582d27 adds a binary notebook (with binary images) into the git history. I would prefer if you can remove it from that commit.

In general, try to learn to create logical commits, so in the past, you should have created just one commit that adds the (binary) notebook, and not to lump it with some other changes in the code. Then it would be trivial to simply remove that commit now. The same thing when you removed the binary parts of the notebook --- you lumped that change with some other changes.

You can edit the commits and remove it the binary notebook. Do you know how to do that? If not, you can either learn that (I think it's a very useful skill, you use git rebase -i and split the commits which touch the notebook, then you squash them), or you can also just squash all of the commits into just few, making sure the binary notebook is removed.

The dictionary representations for best origin and direction parameter(w.r.t hyperplane) are changed to tuple/list form. This is because of discrepancies arising due to {x: a, y: b} and {y: b, x: a} referring to the same object. This change may be temporary or permanent depending on how scalable and robust it proves to be for the final 2D API and the 3D use case.

1. Use map function. 2. Use S.Zero instead of Python zero. 3. Use Tuple instead of List.

intpoly(): 1.decompose(): Returned dictionary also has constant with value 0. 2.integration_reduction(): Compute dot product directly --> only one recursive call. 3.polytope_integrate(): Add support for hyperplane representation. test_intpoly(): 1. Add tests for hyperplane representation. 2. Add tests(two XFAILs) from the paper by Chen et al.(2015)

The user can enter a list of polynomials to integrate over a polytope along with a maximum degree as input. Since all the monomial integrals are calculated once and re-used, the process is faster as compared to the case where each polynomial is computed individually.

On further testing the standard algorithm performed faster than the dynamic one. However, if the terms of the polynomial are such that every monomial's partial derivatives are present previously then dynamic is faster. Therefore, it is used for the multiple case, that is when a maximum degree is supplied and the values for a list of polynomials is to be calculated.

ArifAhmed1995 · 2017-07-08T15:46:14Z

@certik Is it fine now ?

certik · 2017-07-08T17:49:27Z

It's better, but the commit 6ddbad1ca551cd651414b0a0ef14fe40d8403069 still adds binary data. To make it simple, just remove the notebook from 6ddbad1ca551cd651414b0a0ef14fe40d8403069 and just add the notebook in 1bbf991bb4121060f1bce4d431226a89e70a21f2, without the binary images.

1. Add module docstring. 2. Change main function docstring. 3. Add comment for XFAIL tests. 4. Modify decompose function, use parameter as boolean 5. Modify hyperplane_parameters Add notebook

ArifAhmed1995 · 2017-07-08T21:38:28Z

@certik This is done now, I hope. I will work on the intersecting polygon issue for one more day and then start on the 3D use-case. If the intersecting issue isn't solved today, I'll come back to it some other time.

certik

This looks good, +1 from me. Thanks for fixing everything @ArifAhmed1995.

ndotsu · 2017-07-09T07:52:03Z

Arif, The intersecting polygon case should work as is on the examples in the paper - no special treatment is needed so long as the code works for any nonconvex polygon (algorithm only `sees` the edges of the polygon).

…

-suku.

On 7/8/17 11:38 PM, Arif Ahmed wrote: @certik <https://github.com/certik> This is done now, I hope. I will work on the intersecting polygon issue for one more day and then start on the 3D use-case. If the intersecting issue isn't solved today, I'll come back to it some other time. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#12673 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AbwYKgmwNjhLpg35hTMu3__KLn89EGXRks5sL_bYgaJpZM4Nmzx0>.

-- N. Sukumar, Professor Work Phone: (530)754-6415 Department of Civil & Environmental Engg. FAX : (530)752-7872 3159 Ghausi Hall, One Shields Avenue E-mail: nsukumar@ucdavis.edu University of California *GO BLAZERS* Davis, CA 95616 WEB: http://dilbert.engr.ucdavis.edu/~suku

ArifAhmed1995 · 2017-07-09T09:02:53Z

@ndotsu Yes, it will work once the Polygon class in SymPy can handle intersecting polygons without explicitly supplying those intersection points as input. I'm working on it right now.

It's not an issue with the current implementation but with the Polygon class.

ndotsu · 2017-07-09T09:28:39Z

Okay, very good. In case you have not already, do take a look at Section 5.4 in the paper for a discussion on this case.

On 7/9/17 11:02 AM, Arif Ahmed wrote: @ndotsu <https://github.com/ndotsu> Yes, it will work once the |Polygon| class in SymPy can handle intersecting polygons without explicitly supplying those intersection points as input. I'm working on it right now. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#12673 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AbwYKmPd92DdCXaNho_INj5wMWvlPGgFks5sMJdBgaJpZM4Nmzx0>.

-- N. Sukumar, Professor Work Phone: (530)754-6415 Department of Civil & Environmental Engg. FAX : (530)752-7872 3159 Ghausi Hall, One Shields Avenue E-mail: nsukumar@ucdavis.edu University of California *GO BLAZERS* Davis, CA 95616 WEB: http://dilbert.engr.ucdavis.edu/~suku

ArifAhmed1995 · 2017-07-10T23:04:15Z

@smichr I've opened #12931 to resolve the above issue of implicit intersections.

certik mentioned this pull request May 31, 2017

test_heurisch.py fails on Travis #12697

Open

ArifAhmed1995 force-pushed the intpoly branch from f39d353 to f43b56f Compare June 2, 2017 17:29

certik requested changes Jun 3, 2017

View reviewed changes

ArifAhmed1995 force-pushed the intpoly branch 3 times, most recently from 70137ba to 25a4fe1 Compare June 7, 2017 08:55

gxyd reviewed Jun 9, 2017

View reviewed changes

ArifAhmed1995 force-pushed the intpoly branch from c31e06a to 857463f Compare June 10, 2017 21:45

gxyd reviewed Jun 12, 2017

View reviewed changes

certik reviewed Jun 14, 2017

View reviewed changes

ArifAhmed1995 force-pushed the intpoly branch 2 times, most recently from 0671e50 to 1093843 Compare June 18, 2017 05:22

ArifAhmed1995 added 15 commits July 8, 2017 20:47

Complete prototype for 2D use case

fb7f881

Update 2D use case, add tests

a86e5f1

Rectify code, PEP8 fixes and add more tests

238863b

Revert change since list(set(x)) is arbitrary

3858223

Remove implicit import statement

1335ea3

Improve code :

2a4f403

1. Use map function. 2. Use S.Zero instead of Python zero. 3. Use Tuple instead of List.

Add clockwise sort and minor change to reduction function

b74644a

Complete Dynamic implementation

76e085b

Add dynamic programming algorithm

66bf2aa

Remove standard algorithm

ae62870

Major correction in DP algorithm

da832c2

ArifAhmed1995 force-pushed the intpoly branch from 31cf30a to 88e8b64 Compare July 8, 2017 15:39

ArifAhmed1995 force-pushed the intpoly branch from 88e8b64 to 1bbf991 Compare July 8, 2017 17:35

Improve code :

2e0c904

1. Add module docstring. 2. Change main function docstring. 3. Add comment for XFAIL tests. 4. Modify decompose function, use parameter as boolean 5. Modify hyperplane_parameters Add notebook

ArifAhmed1995 force-pushed the intpoly branch from 1bbf991 to 2e0c904 Compare July 8, 2017 18:22

certik approved these changes Jul 9, 2017

View reviewed changes

certik merged commit b2b4bfb into sympy:master Jul 9, 2017

ArifAhmed1995 deleted the intpoly branch July 27, 2017 06:51

Integration of polynomials over 2D polygons #12673

Integration of polynomials over 2D polygons #12673

Conversation

ArifAhmed1995 commented May 25, 2017 • edited

certik commented May 31, 2017

ArifAhmed1995 commented May 31, 2017

certik commented May 31, 2017

certik left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ArifAhmed1995 commented Jun 7, 2017

certik commented Jun 8, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gxyd commented Jun 9, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ArifAhmed1995 commented Jun 9, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gxyd Jun 12, 2017 • edited

Choose a reason for hiding this comment

certik left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

certik commented Jun 14, 2017

ArifAhmed1995 commented Jun 18, 2017

ArifAhmed1995 commented Jun 18, 2017

certik commented Jun 27, 2017

ArifAhmed1995 commented Jun 28, 2017

certik commented Jul 8, 2017

ArifAhmed1995 commented Jul 8, 2017

certik commented Jul 8, 2017

ArifAhmed1995 commented Jul 8, 2017

certik left a comment

Choose a reason for hiding this comment

ndotsu commented Jul 9, 2017 via email

ArifAhmed1995 commented Jul 9, 2017 • edited

ndotsu commented Jul 9, 2017 via email

ArifAhmed1995 commented Jul 10, 2017

ArifAhmed1995 commented May 25, 2017 •

edited

gxyd commented Jun 9, 2017 •

edited

gxyd Jun 12, 2017 •

edited

ArifAhmed1995 commented Jul 9, 2017 •

edited