# sympy/sympy

### Subversion checkout URL

You can clone with HTTPS or Subversion.

# Add the implementation of the incomplite Bell polynomials.#1181

Merged
merged 9 commits into from
+84 −30

### 3 participants

They are also called the second kind of Bell polynomials or "partial"
Bell polynomials. They can be found in Riordan's formulation of di Bruno's
formula for computing an arbitrary order derivative of the composition of
two functions.

P.S.
It's used for series composition.
I gradually begin to transfer the parts of the code from series-prototype necessary in any case.

sympy/functions/combinatorial/numbers.py
 @@ -23,6 +23,7 @@ def _product(a, b): # Dummy symbol used for computing polynomial sequences _sym = Symbol('x') +_symbols = Function('_x')
 goodok added a note Mar 27, 2012 It is worth to use IndexedBase, but the module 'sympy.tensor.indexed' uses 'exp' from sympy.functions, which in turns uses this combinatoric module. to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/functions/combinatorial/numbers.py
 @@ -307,8 +308,26 @@ class bell(Function): n /___ \ k - 1 / k-1 k = 1 + The second kind of Bell polynomials are sometimes called "partial" Bell + polynomials or incomplite Bell polynomials) are defined as:
 ness01 added a note Mar 27, 2012 missing opening bracked, and "incomplete" instead of incomplite goodok added a note Mar 27, 2012 Thanks, will insert these corrections in a few hours. to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/functions/combinatorial/numbers.py
 @@ -307,8 +308,26 @@ class bell(Function): n /___ \ k - 1 / k-1 k = 1 + The second kind of Bell polynomials are sometimes called "partial" Bell + polynomials or incomplite Bell polynomials) are defined as: + + ___ / x \ j / x \ j + \ n! | 1 | 1 | 2 | 2 + B /x , ..., x \ = ) ------------- | --- | | --- | ... + \ 1 n - k + 1/ /___ j ! j ! ... \ 1! / \ 2! / + n,k 1 2 + + j + j + ... = k + 1 2 + + j + 2 * j + 3 * j... = k + 1 2 3 +
 ness01 added a note Mar 27, 2012 Maybe Latex that? goodok added a note Mar 27, 2012 This is the main part of this pull request :) I oriented on the rest of this docstring, but think we can use LaTeX as sphinx can parse it. to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/functions/combinatorial/numbers.py
 ((20 lines not shown)) * bell(n) gives the nth Bell number, B_n * bell(n, x) gives the nth Bell polynomial, B_n(x) + * bell(n, k, (x_1, x_2, ..., x_{n-k+1})) gives the second kind of Bell
 ness01 added a note Mar 27, 2012 "Bell polynomials of the second kind" to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/functions/combinatorial/numbers.py
 ((14 lines not shown)) + where + B_{0,0} = 1; + B_{n,0} = 0; for n=>1 + B_{0,k} = 0; for k=>1 + """ + if (n==0) and (k==0): + return S.One + elif (n==0) or (k==0): + return S.Zero + s = S.Zero + + a = 1 + for m in xrange(1, n-k+2): + s += a*_symbols(m)*bell._bell_incomplete_poly(n-m, k-1) + a = a*(n-m)//m + return s.expand()
 ness01 added a note Mar 27, 2012 I wonder if it might be worth turning this into `polys` code. goodok added a note Mar 27, 2012 Don't mind. As and the above ordinary Bell polynomial existed in master. As you see, in my implementation the unnecessary summation of sympy expressions (polynomial) is used, instead of to keep coefficients of polynomials in dictionary or list and work with them. It is possible, I think, to separate the pure combinatoric calculation of the coefficients, and the obtaining the resulted polynomial. The only problem I see how to form the recursion in more convenient form for usage. Will think. goodok added a note Mar 27, 2012 btw, ordinary Bell polynomial exists and in the mpmath library. goodok added a note Mar 27, 2012 Removed `:math:` About the usage of Polys, I tried, but I tried to use IndexedBase object as abstract generators to poly (in spite of the mentioned modules intersection `exp`), but IndexedBase is non-commutative. I have added the special tests, which I want to be passed. In particular with the implementation light abstract commutative indexed object (`IndexedLight`). ```>>> A = IndexedLight('A') >>> A[1] A[1] >>> Poly(A[1], A) >>> Poly(A[1], A, domain='ZZ[A[1]]') >>> Poly(A[1], A) * A[1] >>> Poly(A[1]**2, A, domain='ZZ[A[1]]')``` it works as desired! So I think that it is possible to include IndexedLight in the core somewhere, so that there is no conflict between the modules as sympy.tensor I have tried to implement it but the bell functions failures. It turns that the problems with commutativity: ```>>> Poly(A[1], A) * A[1] >>> Poly(A[1]**2, A, domain='ZZ[A[1]]') >>> A[1]*Poly(A[1], A) >>> A[1]*Poly(A[1], A, domain='ZZ[A[1]]')``` That is ```>>> Poly(x, x) * x # works Poly(x**2, x, domain='ZZ') >>> x*Poly(x, x) # doesn't works x*Poly(x, x, domain='ZZ')``` Even when I reverse `symbols[m]*self._bell()` to `self._bell()*symbols[m]` it doesn't works, I guess that Mul sorts arguments by hashes. Additionally I have added the tests with tuples instead of variables. This is use-case to avoid substitution, when the bells polynomial is used for calculation of the value of them themselves, in sequences compositions (infinity), and in this case the calculations running directly. When Poly will be used, then it is impossible to use numbers as generators. But now I guess, that it may be to worth to use substitution still. As conclusion, due to the (3) I can't implement Poly now. Only one line can be uncommented for it now. to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/functions/combinatorial/numbers.py
 @@ -361,13 +382,46 @@ def _bell_poly(n, prev): s += a * prev[k-1] return (_sym * s).expand() + @staticmethod + #@assoc_recurrence_memo([[S.One]])
 goodok added a note Mar 27, 2012 It worth to use caching but this "assoc_recurrence_memo" wrapper is not used anywhere in SymPy, and not tested. The aim of this wrapper not only the caching itself, but also to evolve the recursion to circle, to avoid stack overflow. It can be easyToFix issue for someone else after this pushed in. to join this conversation on GitHub. Already have an account? Sign in to comment

SymPy Bot Summary: There were test failures.

@goodok: Please fix the test failures.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYv4sRDA

Interpreter: /usr/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 4a4e88e
branch hash: 7a50637a1ea68bc8264fc45b95fa68498702b699

Automatic review by SymPy Bot.

 goodok `Add the implementation of the incomplete Bell polynomials.` ```They are also called the second kind of Bell polynomials or "partial" Bell polynomials. They can be found in Riordan's formulation of di Bruno's formula for computing an arbitrary order derivative of the composition of two functions.``` `6dac5e5` goodok `Optimize binomial calculation for the incomplete Bell polynomials.` `de86ee7` goodok `Bell polynomials: fix grammar.` `741aade` goodok `Bell polynomials: convert docstring's formulas to latex.` `fcbe5ce`

SymPy Bot Summary: All tests have passed.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYo-wQDA

Interpreter: /usr/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 4a4e88e
branch hash: c401d390057c4a6fcdc79be068641450bbed14a9

Automatic review by SymPy Bot.

 goodok `Bell polynomials: Add tests, remove substitution.` `2d049dc` goodok `Bell polynomials: remove :math:` `5cf3401`

SymPy Bot Summary: All tests have passed.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY0YMRDA

Interpreter: /usr/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 4a4e88e
branch hash: c401d390057c4a6fcdc79be068641450bbed14a9

Automatic review by SymPy Bot.

 goodok `Bell polynomials: test commutativity.` `238bf72`

Regarding polys: I had thought something along the lines of how orthopoly.py is structured should be adopted for all polynomial generation code. But don't worry about it, I think there is no harm adding this later...

SymPy Bot Summary: There were test failures.

@goodok: Please fix the test failures.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY0oMRDA

Interpreter: /usr/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 4a4e88e
branch hash: 238bf72

Automatic review by SymPy Bot.

Thank you for review, Tom.

So postpone with the Poly, (I will play with Poly once more although)

The only one (mentioned) problem is about `@assoc_recurrence_memo`, but as I say this can be postponed too,
I detected that the caching itself is work due to the `@cacheit` `__new__` of the one of the base classes.
But as I mentioned it worth do deal with the `@assoc_recurrence_memo` to envelope recursion to circle.

sympy/functions/combinatorial/numbers.py
 @@ -361,13 +363,46 @@ def _bell_poly(n, prev): s += a * prev[k-1] return (_sym * s).expand() + @staticmethod + #@assoc_recurrence_memo([[S.One]]) + def _bell_incomplete_poly(n, k, symbols): + r""" + The second kind of Bell polynomials (incomplete Bell polynomials). + + Calculated by recurrence formula: + + .. math:: B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1}) = + \sum_{m=1}^{n-k+1} + \x_m \binom{n-1}{m-1} B_{n-m,k-1}(x_1, x_2, \dotsc, x_{n-m-k}) + + where + B_{0,0} = 1; + B_{n,0} = 0; for n=>1
 Collaborator smichr added a note Mar 27, 2012 do you mean `n >= 1`? and the same for k? goodok added a note Mar 27, 2012 Oh, yes, to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/functions/combinatorial/numbers.py
 ((17 lines not shown)) + B_{0,0} = 1; + B_{n,0} = 0; for n=>1 + B_{0,k} = 0; for k=>1 + + """ + if (n==0) and (k==0): + return S.One + elif (n==0) or (k==0): + return S.Zero + #s = Poly(S.Zero, symbols) + s = S.Zero + a = S.One + for m in xrange(1, n-k+2): + s += a*bell._bell_incomplete_poly(n-m, k-1, symbols)*symbols[m] + a = a*(n-m)/m + return s.expand()
 Collaborator smichr added a note Mar 27, 2012 will `expand_mul` suffice here? or simplify's `._mexpand` (for multinomials and muls)? goodok added a note Mar 27, 2012 `expand_mul` will exactly help. I was hoping that Poly embedding will remove `expand` at all. Collaborator smichr added a note Mar 27, 2012 `expand_mul` will exactly help. I was hoping that Poly embedding will remove `expand` at all. If you use polys, it will expand by default. to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/functions/combinatorial/tests/test_comb_numbers.py
 ((7 lines not shown)) from sympy.utilities.pytest import XFAIL x = Symbol('x') +class IndexedLight(Expr):
 Collaborator smichr added a note Mar 27, 2012 I don't understand why these are being added in the test file. What are they doing that can't be done by something already present in sympy? goodok added a note Mar 27, 2012 IndexedBase is non-commutative. So I can't use it here. At the same time I wanted that the case of infinity indexed object will pass always in the future (even if Poly will be embeded). But now I see that Poly fails and for ordinary symbols. Will remove this test. to join this conversation on GitHub. Already have an account? Sign in to comment
 goodok `Bell polynomials: remove IndexedLight from the tests, and other fixes.` `6b1cdec`

This looks very good to me now. I still think we should eventually replace `_bell_incomplete_poly` by a non-recursive version based on polys, but this is not a priority. Running the tests again.

SymPy Bot Summary: There were test failures.

@goodok: Please fix the test failures.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYh_QQDA

Interpreter: /usr/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 023f517
branch hash: 6b1cdec

Automatic review by SymPy Bot.

Only mathml failure, which can be ignored. I'm +1 on this.

we should eventually replace _bell_incomplete_poly by a non-recursive version based on polys

Thanks, I only add to the comments of code about it, and also that the caching mechanism must specialized to cache only coefficients.

Now

``````>>> X = symbols("x:101"); Y = symbols("y:101");
>>> bell(100, 3, X)
2 sec
>>> bell(100, 3, X)
0 sec
>>> bell(100, 3, Y)
again 2 sec
``````

Updated example

sympy/functions/combinatorial/numbers.py
 ((20 lines not shown)) + B_{0,0} = 1; + B_{n,0} = 0; for n>=1 + B_{0,k} = 0; for k>=1 + + """ + if (n==0) and (k==0): + return S.One + elif (n==0) or (k==0): + return S.Zero + #s = Poly(S.Zero, symbols) + s = S.Zero + a = S.One + for m in xrange(1, n-k+2): + s += a*bell._bell_incomplete_poly(n-m, k-1, symbols)*symbols[m] + a = a*(n-m)/m + return expand_mul(s)
 ness01 added a note Mar 28, 2012 Actually this current implementation never uses `symbols[0]`, does it? Is this intended? Collaborator smichr added a note Mar 28, 2012 It should. It would be un-sympy like to drop a symbol that is passed in symbols. Makint this loop on `xrange(n - k + 1)` would then necessitate (at least) the change indicated above at line 330. goodok added a note Mar 28, 2012 No, it is not. Must be fixed to keep the correlation of functions' interface as it is described in the docstring. to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/functions/combinatorial/numbers.py
 @@ -327,6 +327,8 @@ class bell(Function): 846749014511809332450147 >>> bell(4, Symbol('t')) t**4 + 6*t**3 + 7*t**2 + t + >>> bell(6, 2, symbols('x:6'))
 Collaborator smichr added a note Mar 28, 2012 ``````bell(6, 2, symbols('x1:7')) or >>> bell(6, 2, symbols('x:6')) 6*x0*x4 + 15*x1*x3 + 10*x2**2 `````` (See note below regarding `symbols[0]`.) to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/functions/combinatorial/numbers.py
 ((18 lines not shown)) + + where + B_{0,0} = 1; + B_{n,0} = 0; for n>=1 + B_{0,k} = 0; for k>=1 + + """ + if (n==0) and (k==0): + return S.One + elif (n==0) or (k==0): + return S.Zero + #s = Poly(S.Zero, symbols) + s = S.Zero + a = S.One + for m in xrange(1, n-k+2): + s += a*bell._bell_incomplete_poly(n-m, k-1, symbols)*symbols[m]
 ness01 added a note Mar 28, 2012 if you write `s += bell._bell_incomplete_poly(n-m, k-1, symbols)*Poly(a*symbols[m], symbols)` then the code works. However this is very slow. This is probably because the Polys code is not especially good for multivariate polynomials. Or perhaps one has to work at a lower level to get this to work fast. goodok added a note Mar 28, 2012 Thank you. Yes, the lower level seems to be a perspective direction (for this recursive method), with some internal structure. 2-d kind Bell polynomials, as yo can see, have property that the "order" of monomials (if we weight-ed the x_1, x_2 variables) are the same `\$j_1+2j_2+3j_2+\dotsb=n\$`. So the structure is sparse. But may be exists another combinatorial methods (even though recursion by one index, that is 1D), But I did'n find. to join this conversation on GitHub. Already have an account? Sign in to comment

SymPy Bot Summary: All tests have passed.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY9ckRDA

Interpreter: /usr/local/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: b08aa6b
branch hash: 6b1cdec

Automatic review by SymPy Bot.

 goodok `Bell polynomials: fix zero-indexed.` `7ed7161`

SymPy Bot Summary: All tests have passed.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYuuQQDA

Interpreter: /usr/local/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: b08aa6b
branch hash: 7ed7161

Automatic review by SymPy Bot.

I'm +1 on this.

I will push it tomorrow. Then it will be possible to add the issue about the caching.

SymPy Bot Summary: There were test failures.

@goodok: Please fix the test failures.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY2IMRDA

Interpreter: /usr/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 7517d07
branch hash: 7ed7161

Automatic review by SymPy Bot.

(only failure is mathml, so looks good)

SymPy Bot Summary: All tests have passed.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYv-QQDA

Interpreter: /usr/bin/python3 (3.2.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 7517d07
branch hash: 7ed7161

Automatic review by SymPy Bot.

(only failure is mathml, so looks good)

I thought that we must use the special options for this test "#doctest: +NORMALIZE_WHITESPACE -NESS01"

But I think I've seen that you are not alone whom uses the new building of Python.

merged commit `a97fd41` into from
Commits on Mar 27, 2012
1. goodok authored
```They are also called the second kind of Bell polynomials or "partial"
Bell polynomials. They can be found in Riordan's formulation of di Bruno's
formula for computing an arbitrary order derivative of the composition of
two functions.```
2. goodok authored
3. goodok authored
4. goodok authored
5. goodok authored
6. goodok authored
7. goodok authored
8. goodok authored
Commits on Mar 28, 2012
1. goodok authored