import itertools
import math
import functools
import collections
import fractions
[docs]def nCk(n,k):
r'''
Compute the binomial coefficients :math:`\binom nk = \frac{n!}{k!(n-k)!}`
>>> nCk(0,0)
1
>>> nCk(1,0)
1
>>> nCk(4, 2)
6
>>> [[nCk(n,k) for k in range(n+1)] for n in range(7)] == [
... [1],
... [1, 1],
... [1, 2, 1],
... [1, 3, 3, 1],
... [1, 4, 6, 4, 1],
... [1, 5, 10, 10, 5, 1],
... [1, 6, 15, 20, 15, 6, 1] ]
True
'''
f = math.factorial
return f(n) // f(k) // f(n-k)
def tostr(expr):
'''
Convert quadrature result to a human readable string.
>>> q = ([((('x', 0),), 1), ((('x', 1),), 1)], fractions.Fraction(1, 2))
>>> tostr(q)
'1/2 (x0 + x1)'
'''
expression, multiple = expr
terms = []
for t,c in expression:
s = ''
if c > 1:
s += str(c)
s += ' '
s += ''.join(map(str, sum(t,())))
terms.append(s)
s = ' + '.join(terms)
if not multiple == 1:
s = '{}/{} ({})'.format(multiple.numerator, multiple.denominator, s)
return s
[docs]def quadrature(term, verts):
'''
Compute the quadrature of a homogenous polynomial given by term='xxy'
over the vertices verts=(0,1,2) of a simplex.
>>> quadrature('', (0,1))
([((), 1)], Fraction(1, 1))
>>> quadrature('x', (0,1))
([((('x', 0),), 1), ((('x', 1),), 1)], Fraction(1, 2))
>>> quadrature('x', (0,1,2))
([((('x', 0),), 1), ((('x', 1),), 1), ((('x', 2),), 1)], Fraction(1, 3))
>>> quadrature('xx', (0,1)) == (
... [((('x', 0), ('x', 0)), 1),
... ((('x', 0), ('x', 1)), 1),
... ((('x', 1), ('x', 1)), 1)], fractions.Fraction(1, 3))
True
>>> quadrature('xy', (0,1,2)) == (
... [((('x', 0), ('y', 0)), 2),
... ((('x', 0), ('y', 1)), 1),
... ((('x', 0), ('y', 2)), 1),
... ((('x', 1), ('y', 0)), 1),
... ((('x', 1), ('y', 1)), 2),
... ((('x', 1), ('y', 2)), 1),
... ((('x', 2), ('y', 0)), 1),
... ((('x', 2), ('y', 1)), 1),
... ((('x', 2), ('y', 2)), 2)], fractions.Fraction(1, 12))
True
'''
k = len(verts)-1 # dimension of simplex
q = len(term) # order of term
P = itertools.permutations(term)
V = itertools.combinations_with_replacement(verts, q)
Xp = itertools.product(P, V)
X = (tuple(sorted(zip(*x))) for x in Xp)
counter = collections.Counter(X)
gcd = functools.reduce(fractions.gcd, counter.values())
return (sorted([(c, counter[c] // gcd) for c in counter ]),
fractions.Fraction(gcd, math.factorial(q)*nCk(k+q, q)))
def q(term, verts):
return tostr(quadrature(term, verts))
[docs]def test_line_segment():
'''
>>> def q(t): return tostr(quadrature(t, (0,1)))
>>> q('x')
'1/2 (x0 + x1)'
>>> q('xx')
'1/3 (x0x0 + x0x1 + x1x1)'
>>> q('xxx')
'1/4 (x0x0x0 + x0x0x1 + x0x1x1 + x1x1x1)'
>>> q('xy')
'1/6 (2 x0y0 + x0y1 + x1y0 + 2 x1y1)'
>>> q('xxy')
'1/12 (3 x0x0y0 + x0x0y1 + 2 x0x1y0 + 2 x0x1y1 + x1x1y0 + 3 x1x1y1)'
'''
pass
[docs]def test_triangle():
'''
>>> def q(t): return tostr(quadrature(t, (0,1,2)))
>>> q('x')
'1/3 (x0 + x1 + x2)'
>>> q('xx')
'1/6 (x0x0 + x0x1 + x0x2 + x1x1 + x1x2 + x2x2)'
>>> q('xxx')
'1/10 (x0x0x0 + x0x0x1 + x0x0x2 + x0x1x1 + x0x1x2 + x0x2x2 + x1x1x1 + x1x1x2 + x1x2x2 + x2x2x2)'
>>> q('xy')
'1/12 (2 x0y0 + x0y1 + x0y2 + x1y0 + 2 x1y1 + x1y2 + x2y0 + x2y1 + 2 x2y2)'
'''
pass
[docs]def test_tetrahedron():
'''
>>> def q(t): return tostr(quadrature(t, (0,1,2,3)))
>>> q('x')
'1/4 (x0 + x1 + x2 + x3)'
>>> q('xx')
'1/10 (x0x0 + x0x1 + x0x2 + x0x3 + x1x1 + x1x2 + x1x3 + x2x2 + x2x3 + x3x3)'
>>> q('xxx')
'1/20 (x0x0x0 + x0x0x1 + x0x0x2 + x0x0x3 + x0x1x1 + x0x1x2 + x0x1x3 + x0x2x2 + x0x2x3 + x0x3x3 + x1x1x1 + x1x1x2 + x1x1x3 + x1x2x2 + x1x2x3 + x1x3x3 + x2x2x2 + x2x2x3 + x2x3x3 + x3x3x3)'
>>> q('xy')
'1/20 (2 x0y0 + x0y1 + x0y2 + x0y3 + x1y0 + 2 x1y1 + x1y2 + x1y3 + x2y0 + x2y1 + 2 x2y2 + x2y3 + x3y0 + x3y1 + x3y2 + 2 x3y3)'
'''
pass