# -*- coding: utf-8 -*-
#
import math
import numpy
import sympy


def check_degree_1d(
        quadrature, exact, max_degree, tol=1.0e-14
        ):
    val = quadrature(
        lambda x: [x**degree for degree in range(max_degree+1)]
        ).flatten()
    exact_val = numpy.array([exact(degree) for degree in range(max_degree+1)])
    eps = numpy.finfo(float).eps
    # check relative error
    # Allow 1e1 over machine precision.
    alpha = abs(exact_val) * tol + (1e1+tol+exact_val)*eps
    # check where the error is larger than alpha
    is_larger = (exact_val - val) > alpha
    return numpy.where(is_larger)[0] - 1 if any(is_larger) else max_degree


def check_degree(
        quadrature, exact, exponents_creator, max_degree, tol=1.0e-14
        ):
    for degree in range(max_degree+1):
        val = quadrature(
            lambda x, deg=degree: [
                sympy.prod([x[i]**k[i] for i in range(len(k))])
                for k in exponents_creator(deg)
                ])
        exact_val = numpy.array([exact(k) for k in exponents_creator(degree)])
        eps = numpy.finfo(float).eps
        # check relative error
        # The allowance is quite large here, 1e5 over machine precision.
        # Some test fail if lowered, though.
        # TODO increase precision
        alpha = abs(exact_val) * tol + (1.0e5+tol+exact_val)*eps
        is_larger = abs(exact_val - val) > alpha
        if any(is_larger):
            return degree - 1

    return max_degree


def partition(n, d):
    '''Create all nonnegative tuples of length d which sum up to n.
    '''
    # <https://stackoverflow.com/a/45348441/353337>
    def rec(n, d, depth=0):
        if d == depth:
            return [[]]
        return [
            item + [i]
            for i in range(n+1)
            for item in rec(n-i, d, depth=depth+1)
            ]
    return [[n-sum(p)] + p for p in rec(n, d-1)]


def integrate_monomial_over_unit_circle(k):
    '''The integral

    I = \\int_0^2pi \\cos(phi)**k[0] \\sin(phi)**k[1]

    equals 0 if any of k[0], k[1] is odd. If both are even, we make use of

    I = 4 \\int_0^pi/2 \\cos(phi)**k[0] \\sin(phi)**k[1]
      = 2 B(0.5*(k[0]+1), 0.5*(k[1]+1))

    with B(x, y) being the Beta function. It has the representation

    B(x, y) = Gamma(x) Gamma(y) / Gamma(x + y).
    '''
    if any(numpy.array(k) % 2 == 1):
        return 0.0

    return 2 * math.exp(
        + math.lgamma(0.5 * (k[0] + 1))
        + math.lgamma(0.5 * (k[1] + 1))
        - math.lgamma(0.5 * (k[0] + k[1]) + 1)
        )

@@ -749,24 +749,98 @@
 %0A        ):%0A
+    exponents = numpy.concatenate(%5B%0A        exponents_creator(degree)%0A    
     for degr
@@ -864,17 +864,16 @@
 egree+1)
-:
 %0A       
@@ -877,66 +877,106 @@
     
+%5D)%0A%0A    exact_
 val
+s
  = 
-quadrature(%0A            lambda x, deg=degree: %5B%0A    
+numpy.array(%5Bexact(k) for k in exponents%5D)%0A%0A    def fun(x):%0A        return %5B%0A
     
@@ -1039,28 +1039,24 @@
             
- 
-   
 for k in exp
@@ -1065,161 +1065,52 @@
 ents
-_creator(deg)%0A                %5D)%0A        exact_
+%0A            %5D%0A%0A    
 val
+s
  = 
-numpy.array(%5Bexact(k) for k in exponents_creator(degree)%5D)%0A        eps = numpy.finfo(float).eps%0A    
+quadrature(fun)%0A%0A
     
@@ -1132,20 +1132,16 @@
 e error%0A
-    
     # Th
@@ -1205,20 +1205,16 @@
 on.%0A    
-  
-  
 # Some t
@@ -1216,16 +1216,17 @@
 ome test
+s
  fail if
@@ -1243,20 +1243,16 @@
 though.%0A
-    
     # TO
@@ -1273,24 +1273,53 @@
 ecision%0A    
+eps = numpy.finfo(float).eps%0A
     alpha = 
@@ -1323,32 +1323,33 @@
  = abs(exact_val
+s
 ) * tol + (1.0e5
@@ -1354,32 +1354,33 @@
 e5+tol+exact_val
+s
 )*eps%0A        is
@@ -1373,20 +1373,16 @@
 eps%0A    
-   
- 
 is_large
@@ -1390,38 +1390,40 @@
  = abs(exact_val
+s
  - val
+s
 ) %3E alpha%0A      
@@ -1416,27 +1416,34 @@
 %3E alpha%0A
-    
+%0A
     if 
+not numpy.
 any(is_l
@@ -1462,31 +1462,105 @@
     
-    return degree
+return max_degree%0A%0A    k = numpy.where(is_larger)%5B0%5D%0A    degree = numpy.sum(exponents%5Bk%5B0%5D%5D)
  - 1%0A
-%0A
     
@@ -1558,36 +1558,32 @@
  - 1%0A    return 
-max_
 degree%0A%0A%0Adef par

test all at once