# Source code for sfepy.discrete.simplex_cubature

"""
by Andreas Kloeckner.
"""
from __future__ import division

from __future__ import absolute_import
import numpy as nm
import six
from six.moves import range
from functools import reduce

[docs]def generate_decreasing_nonnegative_tuples_summing_to(n, length, min=0,
max=None):
if length == 0:
yield ()
elif length == 1:
if n <= max:
#print "MX", n, max
yield (n,)
else:
return
else:
if max is None or n < max:
max = n

for i in range(min, max+1):
#print "SIG", sig, i
for remainder in generate_decreasing_nonnegative_tuples_summing_to(
n-i, length-1, min, i):
yield (i,) + remainder

[docs]def generate_permutations(original):
"""
Generate all permutations of the list original'.

Nicked from http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/252178
"""
if len(original) <=1:
yield original
else:
for perm in generate_permutations(original[1:]):
for i in range(len(perm)+1):
#nb str[0:1] works in both string and list contexts
yield perm[:i] + original[0:1] + perm[i:]

[docs]def generate_unique_permutations(original):
"""
Generate all unique permutations of the list original'.
"""

for perm in generate_permutations(original):
yield perm

[docs]def wandering_element(length, wanderer=1, landscape=0):
for i in range(length):
yield i*(landscape,) + (wanderer,) + (length-1-i)*(landscape,)

[docs]def factorial(n):
from operator import mul
assert n == int(n)
return reduce(mul, (i for i in range(1,n+1)), 1)

def _extended_euclidean(q, r):
"""
Return a tuple (p, a, b) such that p = aq + br,
where p is the greatest common divisor.
"""

# see [Davenport], Appendix, p. 214

if abs(q) < abs(r):
p, a, b = _extended_euclidean(r, q)
return p, b, a

Q = 1, 0
R = 0, 1

while r:
quot, t = divmod(q, r)
T = Q - quot*R, Q - quot*R
q, r = r, t
Q, R = R, T

return q, Q, Q

def _gcd(q, r):
return _extended_euclidean(q, r)

def _simplify_fraction(a_b):
(a, b) = a_b
gcd = _gcd(a,b)
return (a//gcd, b//gcd)

[docs]def get_simplex_cubature(order, dimension):
r"""
Cubature on an M{n}-simplex.

cf.
A. Grundmann and H.M. Moeller,
Invariant integration formulas for the n-simplex by combinatorial methods,
SIAM J. Numer. Anal.  15 (1978), 282--290.

This cubature rule has both negative and positive weights.  It is exact for
polynomials up to order :math:2s+1, where :math:s is given as *order*.
The integration domain is the unit simplex

.. math::

T_n := \{(x_1, \dots, x_n): x_i \ge -1, \sum_i x_i \le -1\}
"""
s = order
n = dimension
d = exact_to = 2*s+1

points_to_weights = {}

for i in range(s+1):
weight = (-1)**i * 2**(-2*s) \
* (d + n-2*i)**d \
/ factorial(i) \
/ factorial(d+n-i)

for t in generate_decreasing_nonnegative_tuples_summing_to(s-i, n+1):
for beta in generate_unique_permutations(t):
denominator = d+n-2*i
point = tuple(
_simplify_fraction((2*beta_i+1, denominator))
for beta_i in beta)

points_to_weights[point] = points_to_weights.get(point, 0) \
+ weight

vertices = [-1 * nm.ones((n,))] \
+ [nm.array(x)
for x in wandering_element(n, landscape=-1, wanderer=1)]

pos_points = []
pos_weights = []
neg_points = []
neg_weights = []

dim_factor = 2**n
for p, w in six.iteritems(points_to_weights):
real_p = reduce(add, (a/b*v for (a,b),v in zip(p, vertices)))
if w > 0:
pos_points.append(real_p)
pos_weights.append(dim_factor*w)
else:
neg_points.append(real_p)
neg_weights.append(dim_factor*w)

points = nm.array(pos_points + neg_points)
weights = nm.array(pos_weights + neg_weights)

return points, weights, exact_to