From d1fa42b90283351fc2bdaeef52f42302cfc536d8 Mon Sep 17 00:00:00 2001 From: Arif Ahmed Date: Thu, 3 Aug 2017 01:10:40 +0530 Subject: [PATCH 1/7] Basic 3D prototype with constant expr --- sympy/integrals/intpoly.py | 134 ++++++++++++++++++++++---- sympy/integrals/tests/test_intpoly.py | 20 +++- 2 files changed, 136 insertions(+), 18 deletions(-) diff --git a/sympy/integrals/intpoly.py b/sympy/integrals/intpoly.py index fc72bd206feb..031b78fbdec5 100644 --- a/sympy/integrals/intpoly.py +++ b/sympy/integrals/intpoly.py @@ -14,13 +14,90 @@ from functools import cmp_to_key from sympy.core import S, diff, Expr, Symbol - -from sympy.geometry import Segment2D, Polygon, Point -from sympy.abc import x, y +from sympy.geometry import Segment2D, Polygon, Point, Point2D, Point3D +from sympy.abc import x, y, z from sympy.polys.polytools import LC, gcd_list, degree_list +def polytope_integrate3d(poly, expr): + expr = S(expr) + hp_params = hyperplane_parameters(poly) + facets = poly + return main_integrate3d(expr, facets, hp_params) + + +def main_integrate3d(expr, facets, hp_params): + dims = (x, y, z) + dim_length = len(dims) + integral_value = S.Zero + polynomials = decompose(expr) + for deg in polynomials: + poly_contribute = S.Zero + facet_count = 0 + for i, facet in enumerate(facets): + hp = hp_params[i] + if hp[1] == S.Zero: + continue + poly_contribute += polygon_integrate(facet, i, facets, expr, deg) *\ + (hp[1] / norm(tuple(hp[0]))) + facet_count += 1 + poly_contribute /= (dim_length + deg) + integral_value += poly_contribute + return integral_value + + +def polygon_integrate(facet, index, facets, expr, degree): + if expr == S.Zero: + return S.Zero + result = S.Zero + m = len(facets) + x0 = facet[0] + for i in range(len(facet)): + side = (facet[i], facet[(i + 1) % len(facet)]) + for j in range(m): + if j != index and len(set(side).intersection(set(facets[j]))) == 2: + result += distance_to_side(x0, side) *\ + lineseg_integrate(facet, i, side, expr, degree) + + expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] + diff(expr, z) * x0[2] + result += polygon_integrate(facet, index, facets, expr, degree - 1) + result /= (degree + 2) + return result + + +def distance_to_side(point, line_seg): + x1, x2 = line_seg + num = norm(((point[1] - x1[1]) * (point[2] - x2[2]) - + (point[1] - x2[1]) * (point[2] - x1[2]), + (point[0] - x2[0]) * (point[2] - x1[2]) - + (point[0] - x1[0]) * (point[2] - x2[2]), + (point[0] - x1[0]) * (point[1] - x2[1]) - + (point[0] - x2[0]) * (point[1] - x1[1]))) + denom = norm((x2[0] - x1[0], x2[1] - x1[1], x2[2] - x1[2])) + return num/denom + + +def lineseg_integrate(polygon, i, line_seg, expr, degree): + if expr == S.Zero: + return S.Zero + result = S.Zero + x0 = line_seg[0] + distance = norm(tuple([line_seg[1][i] - line_seg[0][i] for i in range(3)])) + #if isinstance(expr, Expr): + # expr_dict = {x: line_seg[0][0], + # y: line_seg[0][1], + # z: line_seg[0][2]} + # result += distance * expr.subs(expr_dict) + #else: + result += distance * expr + + expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] + diff(expr, z) * x0[2] + result += lineseg_integrate(polygon, i, line_seg, expr, degree) + result /= (degree + 1) + return result + + def polytope_integrate(poly, expr, **kwargs): """Integrates homogeneous functions over polytopes. @@ -265,6 +342,8 @@ def left_integral(m, index, facets, x0, expr, gens): value += distance_origin * expr.subs(expr_dict) else: value += distance_origin * expr + else: + integration_reduction_dynamic() return value @@ -370,22 +449,40 @@ def hyperplane_parameters(poly): >>> hyperplane_parameters(Polygon(Point(0, 3), Point(5, 3), Point(1, 1))) [((0, 1), 3), ((1, -2), -1), ((-2, -1), -3)] """ - vertices = list(poly.vertices) + [poly.vertices[0]] # Close the polygon. - params = [None] * (len(vertices) - 1) - for i in range(len(vertices) - 1): - v1 = vertices[i] - v2 = vertices[i + 1] + if isinstance(poly, Polygon): + vertices = list(poly.vertices) + [poly.vertices[0]] # Close the polygon. + params = [None] * (len(vertices) - 1) - a1 = v1[1] - v2[1] - a2 = v2[0] - v1[0] - b = v2[0] * v1[1] - v2[1] * v1[0] + for i in range(len(vertices) - 1): + v1 = vertices[i] + v2 = vertices[i + 1] - factor = gcd_list([a1, a2, b]) + a1 = v1[1] - v2[1] + a2 = v2[0] - v1[0] + b = v2[0] * v1[1] - v2[1] * v1[0] - b = S(b)/factor - a = (S(a1)/factor, S(a2)/factor) - params[i] = (a, b) + factor = gcd_list([a1, a2, b]) + b = S(b) / factor + a = (S(a1) / factor, S(a2) / factor) + params[i] = (a, b) + else: + params = [None] * len(poly) + for i, polygon in enumerate(poly): + v1, v2, v3 = polygon[:3] + v2 = [v2[j] - v1[j] for j in range(0, 3)] + v3 = [v3[j] - v1[j] for j in range(0, 3)] + a = [v3[1] * v2[2] - v3[2] * v2[1], + v3[2] * v2[0] - v3[0] * v2[2], + v3[0] * v2[1] - v3[1] * v2[0]] + + b = sum([a[j] * v1[j] for j in range(0, 3)]) + fac = gcd_list(a) + if fac is S.Zero: + fac = 1 + a = [j / fac for j in a] + b = b / fac + params[i] = (a, b) return params @@ -665,9 +762,12 @@ def norm(point): """ half = S(1)/2 if isinstance(point, tuple): - return (point[0] ** 2 + point[1] ** 2) ** half + return sum([coord ** 2 for coord in point]) ** half elif isinstance(point, Point): - return (point.x ** 2 + point.y ** 2) ** half + if isinstance(point, Point2D): + return (point.x ** 2 + point.y ** 2) ** half + else: + return (point.x ** 2 + point.y ** 2 + point.z) ** half elif isinstance(point, dict): return sum(i**2 for i in point.values()) ** half diff --git a/sympy/integrals/tests/test_intpoly.py b/sympy/integrals/tests/test_intpoly.py index 1e0cc2973143..3bd6d9f77728 100644 --- a/sympy/integrals/tests/test_intpoly.py +++ b/sympy/integrals/tests/test_intpoly.py @@ -5,7 +5,7 @@ from sympy.core import S from sympy.integrals.intpoly import (decompose, best_origin, - polytope_integrate) + polytope_integrate, polytope_integrate3d) from sympy.geometry.line import Segment2D from sympy.geometry.polygon import Polygon @@ -148,6 +148,24 @@ def test_polytope_integrate(): assert result_dict[expr2] == 13062161/27 assert result_dict[expr3] == 1946257153/924 + # Tests for 3D polytopes + unit_cube = [[(0, 1, 0), (1, 1, 0), (1, 1, 1), (0, 1, 1)], + [(0, 1, 1), (1, 1, 1), (1, 0, 1), (0, 0, 1)], + [(1, 1, 1), (1, 1, 0), (1, 0, 0), (1, 0, 1)], + [(0, 0, 1), (1, 0, 1), (1, 0, 0), (0, 0, 0)], + [(0, 1, 1), (0, 0, 1), (0, 0, 0), (0, 1, 0)], + [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)]] + + cube2 = [[(0, 6, 6), (6, 6, 6), (3, 6, 0), (0, 6, 0)], + [(3, 6, 0), (6, 6, 6), (6, 0, 6), (3, 0, 0)], + [(0, 6, 6), (0, 0, 6), (6, 0, 6), (6, 6, 6)], + [(0, 0, 0), (3, 0, 0), (6, 0, 6), (0, 0, 6)], + [(0, 6, 6), (0, 6, 0), (0, 0, 0), (0, 0, 6)], + [(0, 0, 0), (0, 6, 0), (3, 6, 0), (3, 0, 0)]] + + assert polytope_integrate3d(cube2, 1) == S(-162) + assert polytope_integrate3d(unit_cube, 1) == 1 + @XFAIL def test_polytopes_intersecting_sides(): From 535ff1bf194b3e84976743f283a68f76efaa8810 Mon Sep 17 00:00:00 2001 From: Arif Ahmed Date: Thu, 3 Aug 2017 15:38:56 +0530 Subject: [PATCH 2/7] Polynomials work as well --- sympy/integrals/intpoly.py | 19 ++++++++----------- sympy/integrals/tests/test_intpoly.py | 27 ++++++++++++++------------- 2 files changed, 22 insertions(+), 24 deletions(-) diff --git a/sympy/integrals/intpoly.py b/sympy/integrals/intpoly.py index 031b78fbdec5..cf8c7d7254f2 100644 --- a/sympy/integrals/intpoly.py +++ b/sympy/integrals/intpoly.py @@ -84,16 +84,15 @@ def lineseg_integrate(polygon, i, line_seg, expr, degree): result = S.Zero x0 = line_seg[0] distance = norm(tuple([line_seg[1][i] - line_seg[0][i] for i in range(3)])) - #if isinstance(expr, Expr): - # expr_dict = {x: line_seg[0][0], - # y: line_seg[0][1], - # z: line_seg[0][2]} - # result += distance * expr.subs(expr_dict) - #else: - result += distance * expr - + if isinstance(expr, Expr): + expr_dict = {x: line_seg[1][0], + y: line_seg[1][1], + z: line_seg[1][2]} + result += distance * expr.subs(expr_dict) + else: + result += distance * expr expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] + diff(expr, z) * x0[2] - result += lineseg_integrate(polygon, i, line_seg, expr, degree) + result += lineseg_integrate(polygon, i, line_seg, expr, degree - 1) result /= (degree + 1) return result @@ -342,8 +341,6 @@ def left_integral(m, index, facets, x0, expr, gens): value += distance_origin * expr.subs(expr_dict) else: value += distance_origin * expr - else: - integration_reduction_dynamic() return value diff --git a/sympy/integrals/tests/test_intpoly.py b/sympy/integrals/tests/test_intpoly.py index 3bd6d9f77728..ed3c6cf0d84f 100644 --- a/sympy/integrals/tests/test_intpoly.py +++ b/sympy/integrals/tests/test_intpoly.py @@ -10,7 +10,7 @@ from sympy.geometry.line import Segment2D from sympy.geometry.polygon import Polygon from sympy.geometry.point import Point -from sympy.abc import x, y +from sympy.abc import x, y, z from sympy.utilities.pytest import raises, XFAIL @@ -149,22 +149,23 @@ def test_polytope_integrate(): assert result_dict[expr3] == 1946257153/924 # Tests for 3D polytopes - unit_cube = [[(0, 1, 0), (1, 1, 0), (1, 1, 1), (0, 1, 1)], - [(0, 1, 1), (1, 1, 1), (1, 0, 1), (0, 0, 1)], - [(1, 1, 1), (1, 1, 0), (1, 0, 0), (1, 0, 1)], - [(0, 0, 1), (1, 0, 1), (1, 0, 0), (0, 0, 0)], - [(0, 1, 1), (0, 0, 1), (0, 0, 0), (0, 1, 0)], - [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)]] + cube = [[(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)], + [(0, 5, 5), (5, 5, 5), (5, 0, 5), (0, 0, 5)], + [(5, 5, 5), (5, 5, 0), (5, 0, 0), (5, 0, 5)], + [(0, 0, 5), (5, 0, 5), (5, 0, 0), (0, 0, 0)], + [(0, 5, 5), (0, 0, 5), (0, 0, 0), (0, 5, 0)], + [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 5, 0)]] cube2 = [[(0, 6, 6), (6, 6, 6), (3, 6, 0), (0, 6, 0)], - [(3, 6, 0), (6, 6, 6), (6, 0, 6), (3, 0, 0)], - [(0, 6, 6), (0, 0, 6), (6, 0, 6), (6, 6, 6)], - [(0, 0, 0), (3, 0, 0), (6, 0, 6), (0, 0, 6)], - [(0, 6, 6), (0, 6, 0), (0, 0, 0), (0, 0, 6)], - [(0, 0, 0), (0, 6, 0), (3, 6, 0), (3, 0, 0)]] + [(3, 6, 0), (6, 6, 6), (6, 0, 6), (3, 0, 0)], + [(0, 6, 6), (0, 0, 6), (6, 0, 6), (6, 6, 6)], + [(0, 0, 0), (3, 0, 0), (6, 0, 6), (0, 0, 6)], + [(0, 6, 6), (0, 6, 0), (0, 0, 0), (0, 0, 6)], + [(0, 0, 0), (0, 6, 0), (3, 6, 0), (3, 0, 0)]] + assert polytope_integrate3d(cube, x**2 + y**2 + x*y + z**2) == S(15625)/4 assert polytope_integrate3d(cube2, 1) == S(-162) - assert polytope_integrate3d(unit_cube, 1) == 1 + assert polytope_integrate3d(cube, 1) == 125 @XFAIL From 6dcd1f6ff5adcdaec217d922ff0b9975720d788f Mon Sep 17 00:00:00 2001 From: Arif Ahmed Date: Fri, 11 Aug 2017 01:31:00 +0530 Subject: [PATCH 3/7] Begin hyperplane implementation --- sympy/integrals/intpoly.py | 324 +++++++++++++++----------- sympy/integrals/tests/test_intpoly.py | 8 +- 2 files changed, 189 insertions(+), 143 deletions(-) diff --git a/sympy/integrals/intpoly.py b/sympy/integrals/intpoly.py index cf8c7d7254f2..bbdddb9e0a15 100644 --- a/sympy/integrals/intpoly.py +++ b/sympy/integrals/intpoly.py @@ -14,17 +14,99 @@ from functools import cmp_to_key from sympy.core import S, diff, Expr, Symbol -from sympy.geometry import Segment2D, Polygon, Point, Point2D, Point3D +from sympy.geometry import Segment2D, Segment3D, Polygon, Point, Point2D, \ + Point3D from sympy.abc import x, y, z from sympy.polys.polytools import LC, gcd_list, degree_list -def polytope_integrate3d(poly, expr): +def polytope_integrate(poly, expr, **kwargs): + """Integrates homogeneous functions over polytopes. + + This function accepts the polytope in `poly` (currently only polygons are + implemented) and the function in `expr` (currently only + univariate/bivariate polynomials are implemented) and returns the exact + integral of `expr` over `poly`. + Parameters + ========== + poly : The input Polygon. + expr : The input polynomial. + + Optional Parameters: + -------------------- + clockwise : Binary value to sort input points of the polygon clockwise. + max_degree : The maximum degree of any monomial of the input polynomial. + Examples + ======== + >>> from sympy.abc import x, y + >>> from sympy.geometry.polygon import Polygon + >>> from sympy.geometry.point import Point + >>> from sympy.integrals.intpoly import polytope_integrate + >>> polygon = Polygon(Point(0,0), Point(0,1), Point(1,1), Point(1,0)) + >>> polys = [1, x, y, x*y, x**2*y, x*y**2] + >>> expr = x*y + >>> polytope_integrate(polygon, expr) + 1/4 + >>> polytope_integrate(polygon, polys, max_degree=3) + {1: 1, x: 1/2, y: 1/2, x*y: 1/4, x*y**2: 1/6, x**2*y: 1/6} + """ + clockwise = kwargs.get('clockwise', False) + max_degree = kwargs.get('max_degree', None) + + if clockwise is True: + poly = point_sort(poly) + expr = S(expr) - hp_params = hyperplane_parameters(poly) - facets = poly - return main_integrate3d(expr, facets, hp_params) + + if isinstance(poly, Polygon): + # For Vertex Representation + hp_params = hyperplane_parameters(poly) + facets = poly.sides + elif len(poly[0]) == 2: + # For Hyperplane Representation + plen = len(poly) + if len(poly[0][0]) == 2:# 2D case + intersections = [intersection(poly[(i - 1) % plen], poly[i], + "plane2D") + for i in range(0, plen)] + hp_params = poly + lints = len(intersections) + facets = [Segment2D(intersections[i], intersections[(i + 1) % lints]) + for i in range(0, lints)] + else:# 3D case + intersections = [intersection(poly[i], poly, "plane3D") + for i in range(0, plen)] + facets = [point_sort(polygon, poly[i]) for i, polygon in + enumerate(intersections)] + hp_params = poly + else: + # For 3D case. + hp_params = hyperplane_parameters(poly) + facets = poly + + if max_degree is not None: + result = {} + if not isinstance(expr, list): + raise TypeError('Input polynomials must be list of expressions') + result_dict = main_integrate(0, facets, hp_params, max_degree) + for poly in expr: + if poly not in result: + if poly is S.Zero: + result[S.Zero] = S.Zero + continue + integral_value = S.Zero + monoms = decompose(poly, separate=True) + for monom in monoms: + if monom.is_number: + integral_value += result_dict[1] * monom + else: + coeff = LC(monom) + integral_value += result_dict[monom / coeff] * coeff + result[poly] = integral_value + return result + + return main_integrate(expr, facets, hp_params) def main_integrate3d(expr, facets, hp_params): @@ -97,81 +179,6 @@ def lineseg_integrate(polygon, i, line_seg, expr, degree): return result -def polytope_integrate(poly, expr, **kwargs): - """Integrates homogeneous functions over polytopes. - - This function accepts the polytope in `poly` (currently only polygons are - implemented) and the function in `expr` (currently only - univariate/bivariate polynomials are implemented) and returns the exact - integral of `expr` over `poly`. - Parameters - ========== - poly : The input Polygon. - expr : The input polynomial. - - Optional Parameters: - clockwise : Binary value to sort input points of the polygon clockwise. - max_degree : The maximum degree of any monomial of the input polynomial. - Examples - ======== - >>> from sympy.abc import x, y - >>> from sympy.geometry.polygon import Polygon - >>> from sympy.geometry.point import Point - >>> from sympy.integrals.intpoly import polytope_integrate - >>> polygon = Polygon(Point(0,0), Point(0,1), Point(1,1), Point(1,0)) - >>> polys = [1, x, y, x*y, x**2*y, x*y**2] - >>> expr = x*y - >>> polytope_integrate(polygon, expr) - 1/4 - >>> polytope_integrate(polygon, polys, max_degree=3) - {1: 1, x: 1/2, y: 1/2, x*y: 1/4, x*y**2: 1/6, x**2*y: 1/6} - """ - clockwise = kwargs.get('clockwise', False) - max_degree = kwargs.get('max_degree', None) - - if clockwise is True and isinstance(poly, Polygon): - poly = clockwise_sort(poly) - - expr = S(expr) - - if isinstance(poly, Polygon): - # For Vertex Representation - hp_params = hyperplane_parameters(poly) - facets = poly.sides - else: - # For Hyperplane Representation - plen = len(poly) - intersections = [intersection(poly[(i - 1) % plen], poly[i]) - for i in range(0, plen)] - hp_params = poly - lints = len(intersections) - facets = [Segment2D(intersections[i], intersections[(i + 1) % lints]) - for i in range(0, lints)] - - if max_degree is not None: - result = {} - if not isinstance(expr, list): - raise TypeError('Input polynomials must be list of expressions') - result_dict = main_integrate(0, facets, hp_params, max_degree) - for polys in expr: - if polys not in result: - if polys is S.Zero: - result[S.Zero] = S.Zero - continue - integral_value = S.Zero - monoms = decompose(polys, separate=True) - for monom in monoms: - if monom.is_number: - integral_value += result_dict[1] * monom - else: - coeff = LC(monom) - integral_value += result_dict[monom / coeff] * coeff - result[polys] = integral_value - return result - - return main_integrate(expr, facets, hp_params) - - def main_integrate(expr, facets, hp_params, max_degree=None): """Function to translate the problem of integrating univariate/bivariate polynomials over a 2-Polytope to integrating over it's boundary facets. @@ -184,6 +191,7 @@ def main_integrate(expr, facets, hp_params, max_degree=None): hp_params : Hyperplane Parameters of the facets Optional Parameters: + -------------------- max_degree : The maximum degree of any monomial of the input polynomial. >>> from sympy.abc import x, y @@ -325,7 +333,7 @@ def left_integral(m, index, facets, x0, expr, gens): for j in range(0, m): intersect = () if j == (index - 1) % m or j == (index + 1) % m: - intersect = intersection(facets[index], facets[j]) + intersect = intersection(facets[index], facets[j], "segment2D") if intersect: distance_origin = norm(tuple(map(lambda x, y: x - y, intersect, x0))) @@ -467,22 +475,27 @@ def hyperplane_parameters(poly): params = [None] * len(poly) for i, polygon in enumerate(poly): v1, v2, v3 = polygon[:3] - v2 = [v2[j] - v1[j] for j in range(0, 3)] - v3 = [v3[j] - v1[j] for j in range(0, 3)] - a = [v3[1] * v2[2] - v3[2] * v2[1], - v3[2] * v2[0] - v3[0] * v2[2], - v3[0] * v2[1] - v3[1] * v2[0]] - - b = sum([a[j] * v1[j] for j in range(0, 3)]) - fac = gcd_list(a) + normal = cross_product(v1, v2, v3) + b = sum([normal[j] * v1[j] for j in range(0, 3)]) + fac = gcd_list(normal) if fac is S.Zero: fac = 1 - a = [j / fac for j in a] + normal = [j / fac for j in normal] b = b / fac - params[i] = (a, b) + params[i] = (normal, b) return params +def cross_product(v1, v2, v3): + """Returns the cross-product of vectors (v2 - v1) and (v3 - v1) + """ + v2 = [v2[j] - v1[j] for j in range(0, 3)] + v3 = [v3[j] - v1[j] for j in range(0, 3)] + return [v3[1] * v2[2] - v3[2] * v2[1], + v3[2] * v2[0] - v3[0] * v2[2], + v3[0] * v2[1] - v3[1] * v2[0]] + + def best_origin(a, b, lineseg, expr): """Helper method for polytope_integrate. Returns a point on the lineseg whose vector inner product with the @@ -639,7 +652,7 @@ def decompose(expr, separate=False): expr : Polynomial(SymPy expression) Optional Parameters : - + --------------------- separate : If True then simply return a list of the constituent monomials If not then break up the polynomial into constituent homogeneous polynomials. @@ -691,7 +704,7 @@ def decompose(expr, separate=False): return poly_dict -def clockwise_sort(poly): +def point_sort(poly, normal=None, clockwise=True): """Returns the same polygon with points sorted in clockwise order. Note that it's necessary for input points to be sorted in some order @@ -700,46 +713,70 @@ def clockwise_sort(poly): Parameters ========== - poly: 2-Polytope + poly: 2D or 3D Polygon + + Optional Parameters: + --------------------- + normal : The normal of the plane which the 3-Polytope is a part of. + clockwise : Returns points sorted in clockwise order if True and + anti-clockwise if False. Examples ======== - >>> from sympy.integrals.intpoly import clockwise_sort + >>> from sympy.integrals.intpoly import point_sort >>> from sympy.geometry.point import Point - >>> from sympy.geometry.polygon import Polygon - >>> clockwise_sort(Polygon(Point(0, 0), Point(1, 0), Point(1, 1))) - Triangle(Point2D(1, 1), Point2D(1, 0), Point2D(0, 0)) + >>> point_sort([Point(0, 0), Point(1, 0), Point(1, 1)]) + [Point2D(1, 1), Point2D(1, 0), Point2D(0, 0)] """ - n = len(poly.vertices) - vertices = list(poly.vertices) - center = Point(sum(map(lambda vertex: vertex.x, poly.vertices)) / n, - sum(map(lambda vertex: vertex.y, poly.vertices)) / n) + n = len(poly) + if n < 2: + return poly + + order = S(1) if clockwise else S(-1) + dim = len(poly[0]) + if dim == 2: + center = Point(sum(map(lambda vertex: vertex.x, poly)) / n, + sum(map(lambda vertex: vertex.y, poly)) / n) + else: + center = Point(sum(map(lambda vertex: vertex.x, poly)) / n, + sum(map(lambda vertex: vertex.y, poly)) / n, + sum(map(lambda vertex: vertex.z, poly)) / n) - def compareTo(a, b): + def compare(a, b): if a.x - center.x >= S.Zero and b.x - center.x < S.Zero: - return S(-1) + return -order elif a.x - center.x < S.Zero and b.x - center.x >= S.Zero: - return S(1) + return order elif a.x - center.x == S.Zero and b.x - center.x == S.Zero: if a.y - center.y >= S.Zero or b.y - center.y >= S.Zero: - return S(-1) if a.y > b.y else S(1) - return S(-1) if b.y > a.y else S(1) + return -order if a.y > b.y else order + return -order if b.y > a.y else order det = (a.x - center.x) * (b.y - center.y) -\ (b.x - center.x) * (a.y - center.y) if det < S.Zero: - return S(-1) + return -order elif det > S.Zero: - return S(1) + return order first = (a.x - center.x) * (a.x - center.x) +\ (a.y - center.y) * (a.y - center.y) second = (b.x - center.x) * (b.x - center.x) +\ (b.y - center.y) * (b.y - center.y) - return S(-1) if first > second else S(1) + return -order if first > second else order - return Polygon(*sorted(vertices, key=cmp_to_key(compareTo))) + def compare3d(a, b): + det = cross_product(center, a, b) + dot_product = sum([det[i] * normal[i] for i in range(0, 3)]) + if dot_product < S.Zero: + return -order + elif dot_product > S.Zero: + return order + + if dim == 2: + return sorted(poly, key=cmp_to_key(compare)) + return sorted(poly, key=cmp_to_key(compare3d)) def norm(point): @@ -747,7 +784,7 @@ def norm(point): Parameters ========== - point: This denotes a point in the dimensional space. + point: This denotes a point in the dimension_al spac_e. Examples ======== @@ -769,9 +806,8 @@ def norm(point): return sum(i**2 for i in point.values()) ** half -def intersection(lineseg_1, lineseg_2): - """Returns intersection between lines of which - the input line segments are a part of. +def intersection(geom_1, geom_2, intersection_type): + """Returns intersection between geometric objects. Note that this function is meant for use in integration_reduction and at that point in the calling function the lines denoted by the segments surely intersect within segment boundaries. Coincident lines @@ -779,7 +815,7 @@ def intersection(lineseg_1, lineseg_2): Parameters ========== - lineseg_1, lineseg_2: The input line segments + geom_1, geom_2: The input line segments Examples ======== @@ -788,34 +824,44 @@ def intersection(lineseg_1, lineseg_2): >>> from sympy.geometry.line import Segment2D >>> l1 = Segment2D(Point(1, 1), Point(3, 5)) >>> l2 = Segment2D(Point(2, 0), Point(2, 5)) - >>> intersection(l1, l2) + >>> intersection(l1, l2, "segment2D") (2, 3) """ - if isinstance(lineseg_1, Segment2D): - x1, y1 = lineseg_1.points[0] - x2, y2 = lineseg_1.points[1] - x3, y3 = lineseg_2.points[0] - x4, y4 = lineseg_2.points[1] - else: - a1x, a1y = S(lineseg_1[0][0]), S(lineseg_1[0][1]) - a2x, a2y = S(lineseg_2[0][0]), S(lineseg_2[0][1]) - b1, b2 = S(lineseg_1[1]), S(lineseg_2[1]) - - denom = a1x * a2y - a2x * a1y + if intersection_type[:-2] == "segment": + if intersection_type == "segment2D": + x1, y1 = geom_1.points[0] + x2, y2 = geom_1.points[1] + x3, y3 = geom_2.points[0] + x4, y4 = geom_2.points[1] + elif intersection_type == "segment3D": + x1, y1, z1 = geom_1.points[0] + x2, y2, z2 = geom_1.points[1] + x3, y3, z3 = geom_2.points[0] + x4, y4, z4 = geom_2.points[1] + + denom = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4) if denom: - return (S(b1 * a2y - b2 * a1y) / denom, - S(b2 * a1x - b1 * a2x) / denom) - return () + t1 = x1 * y2 - y1 * x2 + t2 = x3 * y4 - x4 * y3 + return (S(t1 * (x3 - x4) - t2 * (x1 - x2)) / denom, + S(t1 * (y3 - y4) - t2 * (y1 - y2)) / denom) + if intersection_type[:-2] == "plane": + if intersection_type == "plane2D":# Intersection of hyperplanes + a1x, a1y = geom_1[0] + a2x, a2y = geom_2[0] + b1, b2 = S(geom_1[1]), S(geom_2[1]) + + denom = a1x * a2y - a2x * a1y + if denom: + return (S(b1 * a2y - b2 * a1y) / denom, + S(b2 * a1x - b1 * a2x) / denom) + elif intersection_type == "plane3D": + a1x, a1y, a1z = geom_1[0] + facets = [] + #for polygon in geom_2: - denom = (x1 - x2)*(y3 - y4) - (y1 - y2)*(x3 - x4) - - if denom: - t1 = x1*y2 - y1*x2 - t2 = x3*y4 - x4*y3 - return(S(t1*(x3 - x4) - t2*(x1 - x2))/denom, - S(t1*(y3 - y4) - t2*(y1 - y2))/denom) - return () + return () def is_vertex(ent): diff --git a/sympy/integrals/tests/test_intpoly.py b/sympy/integrals/tests/test_intpoly.py index ed3c6cf0d84f..4e792b86cdd4 100644 --- a/sympy/integrals/tests/test_intpoly.py +++ b/sympy/integrals/tests/test_intpoly.py @@ -5,7 +5,7 @@ from sympy.core import S from sympy.integrals.intpoly import (decompose, best_origin, - polytope_integrate, polytope_integrate3d) + polytope_integrate) from sympy.geometry.line import Segment2D from sympy.geometry.polygon import Polygon @@ -163,9 +163,9 @@ def test_polytope_integrate(): [(0, 6, 6), (0, 6, 0), (0, 0, 0), (0, 0, 6)], [(0, 0, 0), (0, 6, 0), (3, 6, 0), (3, 0, 0)]] - assert polytope_integrate3d(cube, x**2 + y**2 + x*y + z**2) == S(15625)/4 - assert polytope_integrate3d(cube2, 1) == S(-162) - assert polytope_integrate3d(cube, 1) == 125 + assert polytope_integrate(cube, x**2 + y**2 + x*y + z**2) == S(15625)/4 + assert polytope_integrate(cube2, 1) == S(-162) + assert polytope_integrate(cube, 1) == 125 @XFAIL From 838ad2e72b2180d1ca9bb91cea88bfb76ab8852f Mon Sep 17 00:00:00 2001 From: Arif Ahmed Date: Mon, 14 Aug 2017 02:45:51 +0530 Subject: [PATCH 4/7] Write doctests, tweak code, remove H-representation --- sympy/integrals/intpoly.py | 115 +++++++++++++++++++++++++++++-------- 1 file changed, 91 insertions(+), 24 deletions(-) diff --git a/sympy/integrals/intpoly.py b/sympy/integrals/intpoly.py index bbdddb9e0a15..b672b4853022 100644 --- a/sympy/integrals/intpoly.py +++ b/sympy/integrals/intpoly.py @@ -14,8 +14,7 @@ from functools import cmp_to_key from sympy.core import S, diff, Expr, Symbol -from sympy.geometry import Segment2D, Segment3D, Polygon, Point, Point2D, \ - Point3D +from sympy.geometry import Segment2D, Polygon, Point, Point2D from sympy.abc import x, y, z from sympy.polys.polytools import LC, gcd_list, degree_list @@ -60,30 +59,29 @@ def polytope_integrate(poly, expr, **kwargs): expr = S(expr) if isinstance(poly, Polygon): - # For Vertex Representation + # For Vertex Representation(2D case) hp_params = hyperplane_parameters(poly) facets = poly.sides elif len(poly[0]) == 2: # For Hyperplane Representation plen = len(poly) - if len(poly[0][0]) == 2:# 2D case + if len(poly[0][0]) == 2: # 2D case intersections = [intersection(poly[(i - 1) % plen], poly[i], "plane2D") for i in range(0, plen)] hp_params = poly lints = len(intersections) - facets = [Segment2D(intersections[i], intersections[(i + 1) % lints]) + facets = [Segment2D(intersections[i], + intersections[(i + 1) % lints]) for i in range(0, lints)] - else:# 3D case - intersections = [intersection(poly[i], poly, "plane3D") - for i in range(0, plen)] - facets = [point_sort(polygon, poly[i]) for i, polygon in - enumerate(intersections)] - hp_params = poly + else: + raise NotImplementedError("Integration for H-representation 3D" + "case not implemented yet.") else: # For 3D case. hp_params = hyperplane_parameters(poly) facets = poly + return main_integrate3d(expr, facets, hp_params) if max_degree is not None: result = {} @@ -110,6 +108,30 @@ def polytope_integrate(poly, expr, **kwargs): def main_integrate3d(expr, facets, hp_params): + """Function to translate the problem of integrating uni/bi/trivariate + polynomials over a 3-Polytope to integrating over it's faces. + This is done using Generalized Stokes Theorem and Euler Theorem. + + Parameters + =========== + expr : The input polynomial + facets : Faces of the 3-Polytope + hp_params : Hyperplane Parameters of the facets + + >>> from sympy.abc import x, y + >>> from sympy.integrals.intpoly import main_integrate3d,\ + hyperplane_parameters + >>> triangle = [(0, 0, 0), (1, 1, 1), (0, 0, 2)] + >>> cube = [[(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)],\ + [(0, 5, 5), (5, 5, 5), (5, 0, 5), (0, 0, 5)],\ + [(5, 5, 5), (5, 5, 0), (5, 0, 0), (5, 0, 5)],\ + [(0, 0, 5), (5, 0, 5), (5, 0, 0), (0, 0, 0)],\ + [(0, 5, 5), (0, 0, 5), (0, 0, 0), (0, 5, 0)],\ + [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 5, 0)]] + >>> hp_params = hyperplane_parameters(cube) + >>> main_integrate3d(1, cube, hp_params) + 125 + """ dims = (x, y, z) dim_length = len(dims) integral_value = S.Zero @@ -122,7 +144,7 @@ def main_integrate3d(expr, facets, hp_params): if hp[1] == S.Zero: continue poly_contribute += polygon_integrate(facet, i, facets, expr, deg) *\ - (hp[1] / norm(tuple(hp[0]))) + (hp[1] / norm(tuple(hp[0]))) facet_count += 1 poly_contribute /= (dim_length + deg) integral_value += poly_contribute @@ -130,6 +152,28 @@ def main_integrate3d(expr, facets, hp_params): def polygon_integrate(facet, index, facets, expr, degree): + """Helper function to integrate the input uni/bi/trivariate polynomial + over a certain face of the 3-Polytope. + + Parameters + =========== + facet : A particular face of the 3-Polytope over which `expr` is integrated + index : The index of `facet` in `facets` + facets : Faces of the 3-Polytope + expr : The input polynomial + degree : Degree of `expr` + + >>> from sympy.abc import x, y + >>> from sympy.integrals.intpoly import polygon_integrate + >>> cube = [[(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)],\ + [(0, 5, 5), (5, 5, 5), (5, 0, 5), (0, 0, 5)],\ + [(5, 5, 5), (5, 5, 0), (5, 0, 0), (5, 0, 5)],\ + [(0, 0, 5), (5, 0, 5), (5, 0, 0), (0, 0, 0)],\ + [(0, 5, 5), (0, 0, 5), (0, 0, 0), (0, 5, 0)],\ + [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 5, 0)]] + >>> polygon_integrate(cube[0], 0, cube, 1, 0) + 25 + """ if expr == S.Zero: return S.Zero result = S.Zero @@ -149,6 +193,19 @@ def polygon_integrate(facet, index, facets, expr, degree): def distance_to_side(point, line_seg): + """Helper function to compute the distance between given 3D point and + a line segment. + + Parameters + =========== + point : 3D Point + line_seg : Line Segment + + >>> from sympy.integrals.intpoly import distance_to_side + >>> point = (0, 0, 0) + >>> distance_to_side(point, [(0, 0, 1), (0, 1, 0)]) + sqrt(2)/2 + """ x1, x2 = line_seg num = norm(((point[1] - x1[1]) * (point[2] - x2[2]) - (point[1] - x2[1]) * (point[2] - x1[2]), @@ -157,15 +214,30 @@ def distance_to_side(point, line_seg): (point[0] - x1[0]) * (point[1] - x2[1]) - (point[0] - x2[0]) * (point[1] - x1[1]))) denom = norm((x2[0] - x1[0], x2[1] - x1[1], x2[2] - x1[2])) - return num/denom + return S(num)/denom + +def lineseg_integrate(polygon, index, line_seg, expr, degree): + """Helper function to compute the line integral of `expr` over `line_seg` -def lineseg_integrate(polygon, i, line_seg, expr, degree): + Parameters + =========== + polygon : Face of a 3-Polytope + index : index of line_seg in polygon + line_seg : Line Segment + + >>> from sympy.integrals.intpoly import lineseg_integrate + >>> polygon = [(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)] + >>> line_seg = [(0, 5, 0), (5, 5, 0)] + >>> lineseg_integrate(polygon, 0, line_seg, 1, 0) + 5 + """ if expr == S.Zero: return S.Zero result = S.Zero x0 = line_seg[0] - distance = norm(tuple([line_seg[1][i] - line_seg[0][i] for i in range(3)])) + distance = norm(tuple([line_seg[1][i] - line_seg[0][i] for i in + range(3)])) if isinstance(expr, Expr): expr_dict = {x: line_seg[1][0], y: line_seg[1][1], @@ -174,7 +246,7 @@ def lineseg_integrate(polygon, i, line_seg, expr, degree): else: result += distance * expr expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] + diff(expr, z) * x0[2] - result += lineseg_integrate(polygon, i, line_seg, expr, degree - 1) + result += lineseg_integrate(polygon, index, line_seg, expr, degree - 1) result /= (degree + 1) return result @@ -455,7 +527,7 @@ def hyperplane_parameters(poly): [((0, 1), 3), ((1, -2), -1), ((-2, -1), -3)] """ if isinstance(poly, Polygon): - vertices = list(poly.vertices) + [poly.vertices[0]] # Close the polygon. + vertices = list(poly.vertices) + [poly.vertices[0]] # Close the polygon params = [None] * (len(vertices) - 1) for i in range(len(vertices) - 1): @@ -538,6 +610,7 @@ def best_origin(a, b, lineseg, expr): (0, 3.0) """ a1, b1 = lineseg.points[0] + def x_axis_cut(ls): """Returns the point where the input line segment intersects the x-axis. @@ -847,7 +920,7 @@ def intersection(geom_1, geom_2, intersection_type): return (S(t1 * (x3 - x4) - t2 * (x1 - x2)) / denom, S(t1 * (y3 - y4) - t2 * (y1 - y2)) / denom) if intersection_type[:-2] == "plane": - if intersection_type == "plane2D":# Intersection of hyperplanes + if intersection_type == "plane2D": # Intersection of hyperplanes a1x, a1y = geom_1[0] a2x, a2y = geom_2[0] b1, b2 = S(geom_1[1]), S(geom_2[1]) @@ -856,12 +929,6 @@ def intersection(geom_1, geom_2, intersection_type): if denom: return (S(b1 * a2y - b2 * a1y) / denom, S(b2 * a1x - b1 * a2x) / denom) - elif intersection_type == "plane3D": - a1x, a1y, a1z = geom_1[0] - facets = [] - #for polygon in geom_2: - - return () def is_vertex(ent): From 504cf3de63d5cea77124456d05b55bfd14eaebf9 Mon Sep 17 00:00:00 2001 From: Arif Ahmed Date: Tue, 15 Aug 2017 00:34:30 +0530 Subject: [PATCH 5/7] Documentation fixes --- sympy/integrals/intpoly.py | 87 +++++++++++++++++++++----------------- 1 file changed, 48 insertions(+), 39 deletions(-) diff --git a/sympy/integrals/intpoly.py b/sympy/integrals/intpoly.py index b672b4853022..7b44b7f48408 100644 --- a/sympy/integrals/intpoly.py +++ b/sympy/integrals/intpoly.py @@ -1,6 +1,6 @@ """ Module to implement integration of uni/bivariate polynomials over -2D Polytopes(Polygons). +2D Polytopes and uni/bi/trivariate polynomials over 3D Polytopes. Uses evaluation techniques as described in Chin et al(2015)[1] @@ -21,12 +21,11 @@ def polytope_integrate(poly, expr, **kwargs): - """Integrates homogeneous functions over polytopes. + """Integrates polynomials over 2/3-Polytopes. - This function accepts the polytope in `poly` (currently only polygons are - implemented) and the function in `expr` (currently only - univariate/bivariate polynomials are implemented) and returns the exact - integral of `expr` over `poly`. + This function accepts the polytope in `poly` and the function in `expr` + (uni/bi/trivariate polynomials are implemented) and returns + the exact integral of `expr` over `poly`. Parameters ========== poly : The input Polygon. @@ -34,7 +33,7 @@ def polytope_integrate(poly, expr, **kwargs): Optional Parameters: -------------------- - clockwise : Binary value to sort input points of the polygon clockwise. + clockwise : Binary value to sort input points of the 2-Polytope clockwise. max_degree : The maximum degree of any monomial of the input polynomial. Examples ======== @@ -54,7 +53,11 @@ def polytope_integrate(poly, expr, **kwargs): max_degree = kwargs.get('max_degree', None) if clockwise is True: - poly = point_sort(poly) + if isinstance(poly, Polygon): + poly = point_sort(poly) + else: + raise TypeError("clockwise=True works for only 2-Polytope" + "V-representation input") expr = S(expr) @@ -63,9 +66,9 @@ def polytope_integrate(poly, expr, **kwargs): hp_params = hyperplane_parameters(poly) facets = poly.sides elif len(poly[0]) == 2: - # For Hyperplane Representation + # For Hyperplane Representation(2D case) plen = len(poly) - if len(poly[0][0]) == 2: # 2D case + if len(poly[0][0]) == 2: intersections = [intersection(poly[(i - 1) % plen], poly[i], "plane2D") for i in range(0, plen)] @@ -78,7 +81,7 @@ def polytope_integrate(poly, expr, **kwargs): raise NotImplementedError("Integration for H-representation 3D" "case not implemented yet.") else: - # For 3D case. + # For Vertex Representation(3D case) hp_params = hyperplane_parameters(poly) facets = poly return main_integrate3d(expr, facets, hp_params) @@ -109,8 +112,8 @@ def polytope_integrate(poly, expr, **kwargs): def main_integrate3d(expr, facets, hp_params): """Function to translate the problem of integrating uni/bi/trivariate - polynomials over a 3-Polytope to integrating over it's faces. - This is done using Generalized Stokes Theorem and Euler Theorem. + polynomials over a 3-Polytope to integrating over its faces. + This is done using Generalized Stokes's Theorem and Euler's Theorem. Parameters =========== @@ -123,11 +126,11 @@ def main_integrate3d(expr, facets, hp_params): hyperplane_parameters >>> triangle = [(0, 0, 0), (1, 1, 1), (0, 0, 2)] >>> cube = [[(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)],\ - [(0, 5, 5), (5, 5, 5), (5, 0, 5), (0, 0, 5)],\ - [(5, 5, 5), (5, 5, 0), (5, 0, 0), (5, 0, 5)],\ - [(0, 0, 5), (5, 0, 5), (5, 0, 0), (0, 0, 0)],\ - [(0, 5, 5), (0, 0, 5), (0, 0, 0), (0, 5, 0)],\ - [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 5, 0)]] + [(0, 5, 5), (5, 5, 5), (5, 0, 5), (0, 0, 5)],\ + [(5, 5, 5), (5, 5, 0), (5, 0, 0), (5, 0, 5)],\ + [(0, 0, 5), (5, 0, 5), (5, 0, 0), (0, 0, 0)],\ + [(0, 5, 5), (0, 0, 5), (0, 0, 0), (0, 5, 0)],\ + [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 5, 0)]] >>> hp_params = hyperplane_parameters(cube) >>> main_integrate3d(1, cube, hp_params) 125 @@ -166,11 +169,11 @@ def polygon_integrate(facet, index, facets, expr, degree): >>> from sympy.abc import x, y >>> from sympy.integrals.intpoly import polygon_integrate >>> cube = [[(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)],\ - [(0, 5, 5), (5, 5, 5), (5, 0, 5), (0, 0, 5)],\ - [(5, 5, 5), (5, 5, 0), (5, 0, 0), (5, 0, 5)],\ - [(0, 0, 5), (5, 0, 5), (5, 0, 0), (0, 0, 0)],\ - [(0, 5, 5), (0, 0, 5), (0, 0, 0), (0, 5, 0)],\ - [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 5, 0)]] + [(0, 5, 5), (5, 5, 5), (5, 0, 5), (0, 0, 5)],\ + [(5, 5, 5), (5, 5, 0), (5, 0, 0), (5, 0, 5)],\ + [(0, 0, 5), (5, 0, 5), (5, 0, 0), (0, 0, 0)],\ + [(0, 5, 5), (0, 0, 5), (0, 0, 0), (0, 5, 0)],\ + [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 5, 0)]] >>> polygon_integrate(cube[0], 0, cube, 1, 0) 25 """ @@ -253,8 +256,8 @@ def lineseg_integrate(polygon, index, line_seg, expr, degree): def main_integrate(expr, facets, hp_params, max_degree=None): """Function to translate the problem of integrating univariate/bivariate - polynomials over a 2-Polytope to integrating over it's boundary facets. - This is done using Generalized Stokes Theorem and Euler Theorem. + polynomials over a 2-Polytope to integrating over its boundary facets. + This is done using Generalized Stokes's Theorem and Euler's Theorem. Parameters =========== @@ -429,7 +432,7 @@ def integration_reduction_dynamic(facets, index, a, b, expr, x0, monomial_values, monom_index): """The same integration_reduction function which uses a dynamic programming approach to compute terms by using the values of the integral - the gradient of previous terms. + of previously computed terms. Parameters =========== @@ -486,8 +489,8 @@ def integration_reduction_dynamic(facets, index, a, b, expr, def gradient_terms(binomial_power=0): - """Returns a list of all the possible - monomials between 0 and y**binomial_power + """Returns a list of all the possible monomials between + 0 and y**binomial_power Parameters =========== @@ -569,7 +572,8 @@ def cross_product(v1, v2, v3): def best_origin(a, b, lineseg, expr): - """Helper method for polytope_integrate. + """Helper method for polytope_integrate. Currently not used in the main + algorithm. Returns a point on the lineseg whose vector inner product with the divergence of `expr` yields an expression with the least maximum total power. @@ -778,11 +782,13 @@ def decompose(expr, separate=False): def point_sort(poly, normal=None, clockwise=True): - """Returns the same polygon with points sorted in clockwise order. + """Returns the same polygon with points sorted in clockwise or + anti-clockwise order. Note that it's necessary for input points to be sorted in some order - (clockwise or anti-clockwise) for the algorithm to work. As a convention - algorithm has been implemented keeping clockwise orientation in mind. + (clockwise or anti-clockwise) for the integration algorithm to work. + As a convention algorithm has been implemented keeping clockwise + orientation in mind. Parameters ========== @@ -800,7 +806,6 @@ def point_sort(poly, normal=None, clockwise=True): >>> from sympy.geometry.point import Point >>> point_sort([Point(0, 0), Point(1, 0), Point(1, 1)]) [Point2D(1, 1), Point2D(1, 0), Point2D(0, 0)] - """ n = len(poly) if n < 2: @@ -865,7 +870,6 @@ def norm(point): >>> from sympy.geometry.point import Point >>> norm(Point(2, 7)) sqrt(53) - """ half = S(1)/2 if isinstance(point, tuple): @@ -881,10 +885,12 @@ def norm(point): def intersection(geom_1, geom_2, intersection_type): """Returns intersection between geometric objects. - Note that this function is meant for use in integration_reduction - and at that point in the calling function the lines denoted by the - segments surely intersect within segment boundaries. Coincident lines - are taken to be non-intersecting. + + Note that this function is meant for use in integration_reduction and + at that point in the calling function the lines denoted by the segments + surely intersect within segment boundaries. Coincident lines are taken + to be non-intersecting. Also, the hyperplane intersection for 2D case is + also implemented. Parameters ========== @@ -899,7 +905,10 @@ def intersection(geom_1, geom_2, intersection_type): >>> l2 = Segment2D(Point(2, 0), Point(2, 5)) >>> intersection(l1, l2, "segment2D") (2, 3) - + >>> p1 = ((-1, 0), 0) + >>> p2 = ((0, 1), 1) + >>> intersection(p1, p2, "plane2D") + (0, 1) """ if intersection_type[:-2] == "segment": if intersection_type == "segment2D": From 4793a02a883f1de4a599d543763f806afb305b6d Mon Sep 17 00:00:00 2001 From: Arif Ahmed Date: Fri, 18 Aug 2017 15:41:25 +0530 Subject: [PATCH 6/7] Use economical representation and update notebook --- .../notebooks/IntegrationOverPolytopes.ipynb | 234 ++++++++++++++---- sympy/integrals/intpoly.py | 77 +++--- sympy/integrals/tests/test_intpoly.py | 22 +- 3 files changed, 243 insertions(+), 90 deletions(-) diff --git a/examples/notebooks/IntegrationOverPolytopes.ipynb b/examples/notebooks/IntegrationOverPolytopes.ipynb index eb077de53c2e..d3cb721a183e 100644 --- a/examples/notebooks/IntegrationOverPolytopes.ipynb +++ b/examples/notebooks/IntegrationOverPolytopes.ipynb @@ -20,7 +20,7 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false + "collapsed": true }, "outputs": [], "source": [ @@ -43,15 +43,22 @@ "collapsed": true }, "source": [ - "### polytope_integrate(poly, expr, **kwargs) :\n", - " Pre-processes the input data for integrating univariate/bivariate polynomials over 2-Polytopes.\n", + "### polytope_integrate(poly, expr, **kwargs)\n", + " Integrates polynomials over 2/3-Polytopes.\n", + "\n", + " This function accepts the polytope in `poly` and the function in `expr` (uni/bi/trivariate polynomials are\n", + " implemented) and returns the exact integral of `expr` over `poly`.\n", " \n", - " poly(Polygon) : 2-Polytope\n", - " expr(SymPy expression) : uni/bi-variate polynomial\n", + " Parameters\n", + " ---------------------------------------\n", + " 1. poly(Polygon) : 2/3-Polytope\n", + " 2. expr(SymPy expression) : uni/bi-variate polynomial for 2-Polytope and uni/bi/tri-variate for 3-Polytope\n", " \n", " Optional Parameters\n", - " clockwise(Boolean) : If user is not sure about orientation of vertices and wants to clockwise sort the points.\n", - " max_degree(Integer) : Maximum degree of any monomial of the input polynomial.\n", + " ---------------------------------------\n", + " 1. clockwise(Boolean) : If user is not sure about orientation of vertices of the 2-Polytope and wants\n", + " to clockwise sort the points.\n", + " 2. max_degree(Integer) : Maximum degree of any monomial of the input polynomial. This would require \n", " \n", " #### Examples :" ] @@ -59,9 +66,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "triangle = Polygon(Point(0,0), Point(1,1), Point(1,0))\n", @@ -84,14 +89,164 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### integration_reduction(facets, index, a, b, expr, dims, degree) :\n", - " facets : List of facets that decide the region enclose by 2-Polytope.\n", - " index : The index of the facet with respect to which the integral is supposed to be found.\n", - " a, b : Hyperplane parameters corresponding to facets.\n", - " expr : Uni/Bi-variate Polynomial\n", - " dims : List of symbols denoting axes\n", - " degree : Degree of the homogeneous polynoimal(expr)\n", + "### main_integrate3d(expr, facets, vertices, hp_params)\n", + " Function to translate the problem of integrating uni/bi/tri-variate\n", + " polynomials over a 3-Polytope to integrating over its faces.\n", + " This is done using Generalized Stokes's Theorem and Euler's Theorem.\n", + " \n", + " Parameters\n", + " ------------------\n", + " 1. expr : The input polynomial\n", + " 2. facets : Faces of the 3-Polytope(expressed as indices of `vertices`)\n", + " 3. vertices : Vertices that constitute the Polytope\n", + " 4. hp_params : Hyperplane Parameters of the facets\n", + " \n", + " #### Examples:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cube = [[(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)],\n", + " [2, 6, 7, 3], [3, 7, 5, 1], [7, 6, 4, 5], [1, 5, 4, 0], [3, 1, 0, 2], [0, 4, 6, 2]]\n", + "vertices = cube[0]\n", + "faces = cube[1:]\n", + "hp_params = hyperplane_parameters(faces, vertices)\n", + "main_integrate3d(1, faces, vertices, hp_params)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### polygon_integrate(facet, index, facets, vertices, expr, degree)\n", + " Helper function to integrate the input uni/bi/trivariate polynomial\n", + " over a certain face of the 3-Polytope.\n", + " \n", + " Parameters\n", + " ------------------\n", + " facet : Particular face of the 3-Polytope over which `expr` is integrated\n", + " index : The index of `facet` in `facets`\n", + " facets : Faces of the 3-Polytope(expressed as indices of `vertices`)\n", + " vertices : Vertices that constitute the facet\n", + " expr : The input polynomial\n", + " degree : Degree of `expr`\n", + " \n", + " #### Examples:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0),\n", + " (5, 0, 5), (5, 5, 0), (5, 5, 5)],\n", + " [2, 6, 7, 3], [3, 7, 5, 1], [7, 6, 4, 5], [1, 5, 4, 0],\n", + " [3, 1, 0, 2], [0, 4, 6, 2]]\n", + "facet = cube[1]\n", + "facets = cube[1:]\n", + "vertices = cube[0]\n", + "print(\"Area of polygon < [(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)] > : \", polygon_integrate(facet, 0, facets, vertices, 1, 0))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### distance_to_side(point, line_seg)\n", + "\n", + " Helper function to compute the distance between given 3D point and\n", + " a line segment.\n", + "\n", + " Parameters\n", + " -----------------\n", + " point : 3D Point\n", + " line_seg : Line Segment\n", + " \n", + "#### Examples:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "point = (0, 0, 0)\n", + "distance_to_side(point, [(0, 0, 1), (0, 1, 0)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### lineseg_integrate(polygon, index, line_seg, expr, degree)\n", + " Helper function to compute the line integral of `expr` over `line_seg`\n", + "\n", + " Parameters\n", + " -------------\n", + " polygon : Face of a 3-Polytope\n", + " index : index of line_seg in polygon\n", + " line_seg : Line Segment\n", + "#### Examples :" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "polygon = [(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)]\n", + "line_seg = [(0, 5, 0), (5, 5, 0)]\n", + "print(lineseg_integrate(polygon, 0, line_seg, 1, 0))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### main_integrate(expr, facets, hp_params, max_degree=None)\n", + "\n", + " Function to translate the problem of integrating univariate/bivariate\n", + " polynomials over a 2-Polytope to integrating over its boundary facets.\n", + " This is done using Generalized Stokes's Theorem and Euler's Theorem.\n", + "\n", + " Parameters\n", + " --------------------\n", + " expr : The input polynomial\n", + " facets : Facets(Line Segments) of the 2-Polytope\n", + " hp_params : Hyperplane Parameters of the facets\n", + "\n", + " Optional Parameters:\n", + " --------------------\n", + " max_degree : The maximum degree of any monomial of the input polynomial.\n", " \n", + "#### Examples:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "triangle = Polygon(Point(0, 3), Point(5, 3), Point(1, 1))\n", + "facets = triangle.sides\n", + "hp_params = hyperplane_parameters(triangle)\n", + "print(main_integrate(x**2 + y**2, facets, hp_params))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### integration_reduction(facets, index, a, b, expr, dims, degree)\n", " This is a helper function for polytope_integrate. It relates the result of the integral of a polynomial over a\n", " d-dimensional entity to the result of the same integral of that polynomial over the (d - 1)-dimensional \n", " facet[index].\n", @@ -103,15 +258,22 @@ " this function is a helper one and works for a facet which bounds the polytope(i.e. the intersection point with the\n", " other facets is required), not for an independent line.\n", " \n", + " Parameters\n", + " ------------------\n", + " facets : List of facets that decide the region enclose by 2-Polytope.\n", + " index : The index of the facet with respect to which the integral is supposed to be found.\n", + " a, b : Hyperplane parameters corresponding to facets.\n", + " expr : Uni/Bi-variate Polynomial\n", + " dims : List of symbols denoting axes\n", + " degree : Degree of the homogeneous polynoimal(expr)\n", + " \n", " #### Examples:" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "facets = [Segment2D(Point(0, 0), Point(1, 1)), Segment2D(Point(1, 1), Point(1, 0)), Segment2D(Point(0, 0), Point(1, 0))]\n", @@ -136,9 +298,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "triangle = Polygon(Point(0,0), Point(1,1), Point(1,0))\n", @@ -164,9 +324,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "print(\"Best origin for x**3*y on x + y = 3 : \", best_origin((1,1), 3, Segment2D(Point(0, 3), Point(3, 0)), x**3*y))\n", @@ -190,9 +348,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "print(decompose(1 + x + x**2 + x*y))\n", @@ -216,9 +372,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "print(norm((1, 2)))\n", @@ -243,9 +397,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "print(intersection(Segment2D(Point(0, 0), Point(2, 2)), Segment2D(Point(1, 0), Point(0, 1))))\n", @@ -268,7 +420,7 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false + "collapsed": true }, "outputs": [], "source": [ @@ -295,9 +447,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "hexagon = Polygon(Point(0, 0), Point(-sqrt(3) / 2, 0.5),\n", @@ -326,9 +476,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "expr = x**2\n", @@ -364,7 +512,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.0" + "version": "3.5.2" } }, "nbformat": 4, diff --git a/sympy/integrals/intpoly.py b/sympy/integrals/intpoly.py index 7b44b7f48408..5e5118791921 100644 --- a/sympy/integrals/intpoly.py +++ b/sympy/integrals/intpoly.py @@ -56,8 +56,8 @@ def polytope_integrate(poly, expr, **kwargs): if isinstance(poly, Polygon): poly = point_sort(poly) else: - raise TypeError("clockwise=True works for only 2-Polytope" - "V-representation input") + raise TypeError("clockwise=True works for only 2-Polytope" + "V-representation input") expr = S(expr) @@ -82,9 +82,10 @@ def polytope_integrate(poly, expr, **kwargs): "case not implemented yet.") else: # For Vertex Representation(3D case) - hp_params = hyperplane_parameters(poly) - facets = poly - return main_integrate3d(expr, facets, hp_params) + vertices = poly[0] + facets = poly[1:] + hp_params = hyperplane_parameters(facets, vertices) + return main_integrate3d(expr, facets, vertices, hp_params) if max_degree is not None: result = {} @@ -110,29 +111,29 @@ def polytope_integrate(poly, expr, **kwargs): return main_integrate(expr, facets, hp_params) -def main_integrate3d(expr, facets, hp_params): - """Function to translate the problem of integrating uni/bi/trivariate +def main_integrate3d(expr, facets, vertices, hp_params): + """Function to translate the problem of integrating uni/bi/tri-variate polynomials over a 3-Polytope to integrating over its faces. This is done using Generalized Stokes's Theorem and Euler's Theorem. Parameters =========== expr : The input polynomial - facets : Faces of the 3-Polytope + facets : Faces of the 3-Polytope(expressed as indices of `vertices`) + vertices : Vertices that constitute the Polytope hp_params : Hyperplane Parameters of the facets >>> from sympy.abc import x, y - >>> from sympy.integrals.intpoly import main_integrate3d,\ + >>> from sympy.integrals.intpoly import main_integrate3d, \ hyperplane_parameters - >>> triangle = [(0, 0, 0), (1, 1, 1), (0, 0, 2)] - >>> cube = [[(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)],\ - [(0, 5, 5), (5, 5, 5), (5, 0, 5), (0, 0, 5)],\ - [(5, 5, 5), (5, 5, 0), (5, 0, 0), (5, 0, 5)],\ - [(0, 0, 5), (5, 0, 5), (5, 0, 0), (0, 0, 0)],\ - [(0, 5, 5), (0, 0, 5), (0, 0, 0), (0, 5, 0)],\ - [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 5, 0)]] - >>> hp_params = hyperplane_parameters(cube) - >>> main_integrate3d(1, cube, hp_params) + >>> cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0),\ + (5, 0, 5), (5, 5, 0), (5, 5, 5)],\ + [2, 6, 7, 3], [3, 7, 5, 1], [7, 6, 4, 5], [1, 5, 4, 0],\ + [3, 1, 0, 2], [0, 4, 6, 2]] + >>> vertices = cube[0] + >>> faces = cube[1:] + >>> hp_params = hyperplane_parameters(faces, vertices) + >>> main_integrate3d(1, faces, vertices, hp_params) 125 """ dims = (x, y, z) @@ -146,7 +147,8 @@ def main_integrate3d(expr, facets, hp_params): hp = hp_params[i] if hp[1] == S.Zero: continue - poly_contribute += polygon_integrate(facet, i, facets, expr, deg) *\ + poly_contribute += polygon_integrate(facet, i, facets, + vertices, expr, deg) *\ (hp[1] / norm(tuple(hp[0]))) facet_count += 1 poly_contribute /= (dim_length + deg) @@ -154,43 +156,47 @@ def main_integrate3d(expr, facets, hp_params): return integral_value -def polygon_integrate(facet, index, facets, expr, degree): +def polygon_integrate(facet, index, facets, vertices, expr, degree): """Helper function to integrate the input uni/bi/trivariate polynomial over a certain face of the 3-Polytope. Parameters =========== - facet : A particular face of the 3-Polytope over which `expr` is integrated + facet : Particular face of the 3-Polytope over which `expr` is integrated index : The index of `facet` in `facets` - facets : Faces of the 3-Polytope + facets : Faces of the 3-Polytope(expressed as indices of `vertices`) + vertices : Vertices that constitute the facet expr : The input polynomial degree : Degree of `expr` >>> from sympy.abc import x, y >>> from sympy.integrals.intpoly import polygon_integrate - >>> cube = [[(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)],\ - [(0, 5, 5), (5, 5, 5), (5, 0, 5), (0, 0, 5)],\ - [(5, 5, 5), (5, 5, 0), (5, 0, 0), (5, 0, 5)],\ - [(0, 0, 5), (5, 0, 5), (5, 0, 0), (0, 0, 0)],\ - [(0, 5, 5), (0, 0, 5), (0, 0, 0), (0, 5, 0)],\ - [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 5, 0)]] - >>> polygon_integrate(cube[0], 0, cube, 1, 0) + >>> cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0),\ + (5, 0, 5), (5, 5, 0), (5, 5, 5)],\ + [2, 6, 7, 3], [3, 7, 5, 1], [7, 6, 4, 5], [1, 5, 4, 0],\ + [3, 1, 0, 2], [0, 4, 6, 2]] + >>> facet = cube[1] + >>> facets = cube[1:] + >>> vertices = cube[0] + >>> polygon_integrate(facet, 0, facets, vertices, 1, 0) 25 """ if expr == S.Zero: return S.Zero result = S.Zero m = len(facets) - x0 = facet[0] + x0 = vertices[facet[0]] for i in range(len(facet)): side = (facet[i], facet[(i + 1) % len(facet)]) for j in range(m): if j != index and len(set(side).intersection(set(facets[j]))) == 2: + side = (vertices[side[0]], vertices[side[1]]) result += distance_to_side(x0, side) *\ lineseg_integrate(facet, i, side, expr, degree) expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] + diff(expr, z) * x0[2] - result += polygon_integrate(facet, index, facets, expr, degree - 1) + result += polygon_integrate(facet, index, facets, vertices, expr, + degree - 1) result /= (degree + 2) return result @@ -513,7 +519,7 @@ def gradient_terms(binomial_power=0): return terms -def hyperplane_parameters(poly): +def hyperplane_parameters(poly, vertices=None): """A helper function to return the hyperplane parameters of which the facets of the polygon are a part of. Currently works for only 2-Polytopes. @@ -521,6 +527,9 @@ def hyperplane_parameters(poly): ========== poly : The input Polygon + Optional Parameters + ---------------------- + vertices : Vertex indices of faces of 3-Polytope Examples ======== >>> from sympy.geometry.point import Point @@ -549,7 +558,7 @@ def hyperplane_parameters(poly): else: params = [None] * len(poly) for i, polygon in enumerate(poly): - v1, v2, v3 = polygon[:3] + v1, v2, v3 = [vertices[vertex] for vertex in polygon[:3]] normal = cross_product(v1, v2, v3) b = sum([normal[j] * v1[j] for j in range(0, 3)]) fac = gcd_list(normal) @@ -997,4 +1006,4 @@ def plot_polynomial(expr): if len(gens) == 2: plot3d(expr) else: - plot(expr) + plot(expr) \ No newline at end of file diff --git a/sympy/integrals/tests/test_intpoly.py b/sympy/integrals/tests/test_intpoly.py index 4e792b86cdd4..abfe0358ae6a 100644 --- a/sympy/integrals/tests/test_intpoly.py +++ b/sympy/integrals/tests/test_intpoly.py @@ -149,19 +149,15 @@ def test_polytope_integrate(): assert result_dict[expr3] == 1946257153/924 # Tests for 3D polytopes - cube = [[(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)], - [(0, 5, 5), (5, 5, 5), (5, 0, 5), (0, 0, 5)], - [(5, 5, 5), (5, 5, 0), (5, 0, 0), (5, 0, 5)], - [(0, 0, 5), (5, 0, 5), (5, 0, 0), (0, 0, 0)], - [(0, 5, 5), (0, 0, 5), (0, 0, 0), (0, 5, 0)], - [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 5, 0)]] - - cube2 = [[(0, 6, 6), (6, 6, 6), (3, 6, 0), (0, 6, 0)], - [(3, 6, 0), (6, 6, 6), (6, 0, 6), (3, 0, 0)], - [(0, 6, 6), (0, 0, 6), (6, 0, 6), (6, 6, 6)], - [(0, 0, 0), (3, 0, 0), (6, 0, 6), (0, 0, 6)], - [(0, 6, 6), (0, 6, 0), (0, 0, 0), (0, 0, 6)], - [(0, 0, 0), (0, 6, 0), (3, 6, 0), (3, 0, 0)]] + cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0), + (5, 0, 5), (5, 5, 0), (5, 5, 5)], + [2, 6, 7, 3], [3, 7, 5, 1], [7, 6, 4, 5], [1, 5, 4, 0], + [3, 1, 0, 2], [0, 4, 6, 2]] + + cube2 = [[(0, 0, 0), (0, 6, 6), (6, 6, 6), (3, 6, 0), + (0, 6, 0), (6, 0, 6), (3, 0, 0), (0, 0, 6)], + [1, 2, 3, 4], [3, 2, 5, 6], [1, 7, 5, 2], [0, 6, 5, 7], + [1, 4, 0, 7], [0, 4, 3, 6]] assert polytope_integrate(cube, x**2 + y**2 + x*y + z**2) == S(15625)/4 assert polytope_integrate(cube2, 1) == S(-162) From ad1cacbba9baf0b88006bed651a182489025e194 Mon Sep 17 00:00:00 2001 From: Arif Ahmed Date: Tue, 22 Aug 2017 19:32:54 +0530 Subject: [PATCH 7/7] Add 3D tests from Chin et al(2015) --- sympy/integrals/intpoly.py | 67 ++++++++++++++------------- sympy/integrals/tests/test_intpoly.py | 43 +++++++++++------ 2 files changed, 64 insertions(+), 46 deletions(-) diff --git a/sympy/integrals/intpoly.py b/sympy/integrals/intpoly.py index 5e5118791921..fbdad5e4d03c 100644 --- a/sympy/integrals/intpoly.py +++ b/sympy/integrals/intpoly.py @@ -134,7 +134,7 @@ def main_integrate3d(expr, facets, vertices, hp_params): >>> faces = cube[1:] >>> hp_params = hyperplane_parameters(faces, vertices) >>> main_integrate3d(1, faces, vertices, hp_params) - 125 + -125 """ dims = (x, y, z) dim_length = len(dims) @@ -147,8 +147,8 @@ def main_integrate3d(expr, facets, vertices, hp_params): hp = hp_params[i] if hp[1] == S.Zero: continue - poly_contribute += polygon_integrate(facet, i, facets, - vertices, expr, deg) *\ + pi = polygon_integrate(facet, hp, i, facets, vertices, expr, deg) + poly_contribute += pi *\ (hp[1] / norm(tuple(hp[0]))) facet_count += 1 poly_contribute /= (dim_length + deg) @@ -156,7 +156,7 @@ def main_integrate3d(expr, facets, vertices, hp_params): return integral_value -def polygon_integrate(facet, index, facets, vertices, expr, degree): +def polygon_integrate(facet, hp_param, index, facets, vertices, expr, degree): """Helper function to integrate the input uni/bi/trivariate polynomial over a certain face of the 3-Polytope. @@ -178,9 +178,10 @@ def polygon_integrate(facet, index, facets, vertices, expr, degree): >>> facet = cube[1] >>> facets = cube[1:] >>> vertices = cube[0] - >>> polygon_integrate(facet, 0, facets, vertices, 1, 0) - 25 + >>> polygon_integrate(facet, [(0, 1, 0), 5], 0, facets, vertices, 1, 0) + -25 """ + if expr == S.Zero: return S.Zero result = S.Zero @@ -191,19 +192,20 @@ def polygon_integrate(facet, index, facets, vertices, expr, degree): for j in range(m): if j != index and len(set(side).intersection(set(facets[j]))) == 2: side = (vertices[side[0]], vertices[side[1]]) - result += distance_to_side(x0, side) *\ - lineseg_integrate(facet, i, side, expr, degree) + result += distance_to_side(x0, side, hp_param[0]) *\ + lineseg_integrate(facet, i, side, expr, degree) - expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] + diff(expr, z) * x0[2] - result += polygon_integrate(facet, index, facets, vertices, expr, - degree - 1) + expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] +\ + diff(expr, z) * x0[2] + result += polygon_integrate(facet, hp_param, index, facets, vertices, + expr, degree - 1) result /= (degree + 2) return result -def distance_to_side(point, line_seg): - """Helper function to compute the distance between given 3D point and - a line segment. +def distance_to_side(point, line_seg, A): + """Helper function to compute the signed distance between given 3D point + and a line segment. Parameters =========== @@ -212,18 +214,19 @@ def distance_to_side(point, line_seg): >>> from sympy.integrals.intpoly import distance_to_side >>> point = (0, 0, 0) - >>> distance_to_side(point, [(0, 0, 1), (0, 1, 0)]) - sqrt(2)/2 + >>> distance_to_side(point, [(0, 0, 1), (0, 1, 0)], (1, 0, 0)) + -sqrt(2)/2 """ x1, x2 = line_seg - num = norm(((point[1] - x1[1]) * (point[2] - x2[2]) - - (point[1] - x2[1]) * (point[2] - x1[2]), - (point[0] - x2[0]) * (point[2] - x1[2]) - - (point[0] - x1[0]) * (point[2] - x2[2]), - (point[0] - x1[0]) * (point[1] - x2[1]) - - (point[0] - x2[0]) * (point[1] - x1[1]))) - denom = norm((x2[0] - x1[0], x2[1] - x1[1], x2[2] - x1[2])) - return S(num)/denom + rev_normal = [-1 * S(i)/norm(A) for i in A] + vector = [x2[i] - x1[i] for i in range(0, 3)] + vector = [vector[i]/norm(vector) for i in range(0, 3)] + + n_side = cross_product((0, 0, 0), rev_normal, vector) + vectorx0 = [line_seg[0][i] - point[i] for i in range(0, 3)] + dot_product = sum([vectorx0[i] * n_side[i] for i in range(0, 3)]) + + return dot_product def lineseg_integrate(polygon, index, line_seg, expr, degree): @@ -254,7 +257,8 @@ def lineseg_integrate(polygon, index, line_seg, expr, degree): result += distance * expr.subs(expr_dict) else: result += distance * expr - expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] + diff(expr, z) * x0[2] + expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] +\ + diff(expr, z) * x0[2] result += lineseg_integrate(polygon, index, line_seg, expr, degree - 1) result /= (degree + 1) return result @@ -478,8 +482,6 @@ def integration_reduction_dynamic(facets, index, a, b, expr, return expr if not expr.is_number: - a, b = (S(a[0]), S(a[1])), S(b) - x_index = monom_index - max_y_degree + \ x_degree - 2 if x_degree >= 1 else 0 y_index = monom_index - 1 if y_degree >= 1 else 0 @@ -572,12 +574,13 @@ def hyperplane_parameters(poly, vertices=None): def cross_product(v1, v2, v3): """Returns the cross-product of vectors (v2 - v1) and (v3 - v1) + That is : (v2 - v1) X (v3 - v1) """ v2 = [v2[j] - v1[j] for j in range(0, 3)] v3 = [v3[j] - v1[j] for j in range(0, 3)] - return [v3[1] * v2[2] - v3[2] * v2[1], - v3[2] * v2[0] - v3[0] * v2[2], - v3[0] * v2[1] - v3[1] * v2[0]] + return [v3[2] * v2[1] - v3[1] * v2[2], + v3[0] * v2[2] - v3[2] * v2[0], + v3[1] * v2[0] - v3[0] * v2[1]] def best_origin(a, b, lineseg, expr): @@ -881,7 +884,7 @@ def norm(point): sqrt(53) """ half = S(1)/2 - if isinstance(point, tuple): + if isinstance(point, (list, tuple)): return sum([coord ** 2 for coord in point]) ** half elif isinstance(point, Point): if isinstance(point, Point2D): @@ -1006,4 +1009,4 @@ def plot_polynomial(expr): if len(gens) == 2: plot3d(expr) else: - plot(expr) \ No newline at end of file + plot(expr) diff --git a/sympy/integrals/tests/test_intpoly.py b/sympy/integrals/tests/test_intpoly.py index abfe0358ae6a..7f8df683bc2c 100644 --- a/sympy/integrals/tests/test_intpoly.py +++ b/sympy/integrals/tests/test_intpoly.py @@ -12,7 +12,7 @@ from sympy.geometry.point import Point from sympy.abc import x, y, z -from sympy.utilities.pytest import raises, XFAIL +from sympy.utilities.pytest import XFAIL def test_decompose(): @@ -144,25 +144,40 @@ def test_polytope_integrate(): expr3 = x**10 + x**9*y + x**8*y**2 + x**5*y**5 polys.extend((expr1, expr2, expr3)) result_dict = polytope_integrate(tri, polys, max_degree=10) - assert result_dict[expr1] == 615780107/594 - assert result_dict[expr2] == 13062161/27 - assert result_dict[expr3] == 1946257153/924 + assert result_dict[expr1] == S(615780107)/594 + assert result_dict[expr2] == S(13062161)/27 + assert result_dict[expr3] == S(1946257153)/924 # Tests for 3D polytopes - cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0), - (5, 0, 5), (5, 5, 0), (5, 5, 5)], - [2, 6, 7, 3], [3, 7, 5, 1], [7, 6, 4, 5], [1, 5, 4, 0], - [3, 1, 0, 2], [0, 4, 6, 2]] - - cube2 = [[(0, 0, 0), (0, 6, 6), (6, 6, 6), (3, 6, 0), + cube1 = [[(0, 0, 0), (0, 6, 6), (6, 6, 6), (3, 6, 0), (0, 6, 0), (6, 0, 6), (3, 0, 0), (0, 0, 6)], [1, 2, 3, 4], [3, 2, 5, 6], [1, 7, 5, 2], [0, 6, 5, 7], [1, 4, 0, 7], [0, 4, 3, 6]] + assert polytope_integrate(cube1, 1) == S(162) - assert polytope_integrate(cube, x**2 + y**2 + x*y + z**2) == S(15625)/4 - assert polytope_integrate(cube2, 1) == S(-162) - assert polytope_integrate(cube, 1) == 125 - + # 3D Test cases in Chin et al(2015) + cube2 = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0), + (5, 0, 5), (5, 5, 0), (5, 5, 5)], + [3, 7, 6, 2], [1, 5, 7, 3], [5, 4, 6, 7], [0, 4, 5, 1], + [2, 0, 1, 3], [2, 6, 4, 0]] + + cube3 = [[(0, 0, 0), (5, 0, 0), (5, 4, 0), (3, 2, 0), (3, 5, 0), + (0, 5, 0), (0, 0, 5), (5, 0, 5), (5, 4, 5), (3, 2, 5), + (3, 5, 5), (0, 5, 5)], + [6, 11, 5, 0], [1, 7, 6, 0], [5, 4, 3, 2, 1, 0], [11, 10, 4, 5], + [10, 9, 3, 4], [9, 8, 2, 3], [8, 7, 1, 2], [7, 8, 9, 10, 11, 6]] + + cube4 = [[(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1), + (S(1) / 4, S(1) / 4, S(1) / 4)], + [0, 2, 1], [1, 3, 0], [4, 2, 3], [4, 3, 1], + [0, 1, 2], [2, 4, 1], [0, 3, 2]] + + assert polytope_integrate(cube2, x ** 2 + y ** 2 + x * y + z ** 2) ==\ + S(15625)/4 + assert polytope_integrate(cube3, x ** 2 + y ** 2 + x * y + z ** 2) ==\ + S(33835) / 12 + assert polytope_integrate(cube4, x ** 2 + y ** 2 + x * y + z ** 2) ==\ + S(37) / 960 @XFAIL def test_polytopes_intersecting_sides():