{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some background can be found on\n",
    "http://oeis.org/wiki/User:Peter_Luschny/BellTransform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%load_ext sage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "oeis_search_on = False # global option for the notebook\n",
    "\n",
    "def OEIS_ID(M, num=4):\n",
    "    global oeis_search_on\n",
    "    if not oeis_search_on: return ''\n",
    "    L = [M[i,j] for i in range(M.nrows()) for j in range(M.ncols())]  \n",
    "    V = [l for l in L if (l <> 0)]\n",
    "    if V == []: return []\n",
    "    s = ''\n",
    "    for v in V: s += v.str() + ' '  \n",
    "    R = oeis(s, max_results=num) \n",
    "    W = [str(r)[0:7] for r in R]\n",
    "    t = ''\n",
    "    for w in W: t += 'id:'+w+'|'\n",
    "    return t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def bell_transform(n, a):\n",
    "    row = []\n",
    "    fn = factorial(n)\n",
    "    for k in (0..n):\n",
    "        result = 0\n",
    "        for p in Partitions(n, length=k):\n",
    "            factorial_product = 1\n",
    "            power_factorial_product = 1\n",
    "            for part, count in p.to_exp_dict().iteritems():\n",
    "                factorial_product *= factorial(count)\n",
    "                power_factorial_product *= factorial(part)**count\n",
    "            coefficient = fn//(factorial_product*power_factorial_product)\n",
    "            result += coefficient*prod([a[i-1] for i in p])\n",
    "        row.append(result)\n",
    "    return row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def bell_matrix(generator, dim):\n",
    "    G = [generator(k) for k in srange(dim)]\n",
    "    row = lambda n: bell_transform(n, G) \n",
    "    return matrix(ZZ, [row(n)+[0]*(dim-n-1) for n in srange(dim)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "B = bell_matrix(factorial, 8)\n",
    "print B; print OEIS_ID(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "B = bell_matrix(factorial, 8)\n",
    "print B; print OEIS_ID(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "B = bell_matrix(factorial, 8).inverse()\n",
    "print B; print OEIS_ID(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The unsigned inverse matrix.\n",
    "def inverse_bell_matrix(generator, dim): \n",
    "    M = bell_matrix(generator, dim).inverse()\n",
    "    return matrix(ZZ, dim, lambda n,k: (-1)^(n-k)*M[n,k])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "B = inverse_bell_matrix(factorial, 8)\n",
    "print B; print OEIS_ID(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "B = bell_matrix(lambda n: n+1, 8).inverse()\n",
    "print B; print OEIS_ID(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "B = bell_matrix(lambda n: n+1, 8)\n",
    "print B; print OEIS_ID(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# [A023531, A048993, A075497, A075498, A075499, A075500, ...]\n",
    "for n in range(7):\n",
    "    B = bell_matrix(lambda k: n^k, 7)\n",
    "    print \"Bell matrix generated by power function\", n\n",
    "    print OEIS_ID(B)    \n",
    "    print B "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# [A023531, A132393, A039683, A051141, A051142, A051150, ...]\n",
    "for n in range(7):\n",
    "    B = bell_matrix(lambda k: n^k, 7).inverse()\n",
    "    print \"Inverse Bell matrix generated by power function\", n\n",
    "    print OEIS_ID(B)    \n",
    "    print B "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "multifactorial = lambda a,b: lambda n: prod(a*k + b for k in (0..n-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[multifactorial(2,2)(n) for n in range(8)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A132393, \n",
    "for a in (1..4):\n",
    "    for b in (1..a):\n",
    "        B = bell_matrix(multifactorial(a,b), 7)\n",
    "        print \"Bell matrix generated by multifactorial\", (a,b)\n",
    "        print OEIS_ID(B)    \n",
    "        print B "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for a in (1..4):\n",
    "    for b in (1..a):\n",
    "        B = inverse_bell_matrix(multifactorial(a,b), 7)\n",
    "        print \"Inverse Bell matrix generated by multifactorial\", (a,b)\n",
    "        print OEIS_ID(B)\n",
    "        print B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "shifted_factorial = lambda sh: lambda n: factorial(n + sh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# [A132393, A111596 (Lah), A136656, ...]\n",
    "for sh in range(4):\n",
    "    B = bell_matrix(shifted_factorial(sh), 7)\n",
    "    print \"Bell matrix generated by shifted factorial\", sh\n",
    "    print OEIS_ID(B)\n",
    "    print B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A048993, A111596 (Lah),\n",
    "for sh in range(2):\n",
    "    B = inverse_bell_matrix(shifted_factorial(sh), 7)\n",
    "    print \"Inverse Bell matrix generated by shifted factorial\", sh\n",
    "    print OEIS_ID(B)\n",
    "    print B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "risingfactorial = lambda n: lambda k: rising_factorial(n, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A023531, A132393, A111596, A046089, A049352, A049353, A049374, A134141, ... \n",
    "for n in (0..7):\n",
    "    B = bell_matrix(risingfactorial(n), 7)\n",
    "    print \"Bell matrix generated by rising factorial(n, k)\", n\n",
    "    print OEIS_ID(B)\n",
    "    print B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for n in (0..7):\n",
    "    B = inverse_bell_matrix(risingfactorial(n), 7)\n",
    "    print \"Inverse Bell matrix generated by rising factorial(n, k)\", n\n",
    "    print OEIS_ID(B)\n",
    "    print B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "monotonicfactorial = lambda r: lambda n: rising_factorial(r, n)/factorial(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for n in (0..4):\n",
    "    B = bell_matrix(monotonicfactorial(n), 7)\n",
    "    print \"Bell matrix generated by monotonic factorial(r, n)\", n\n",
    "    print OEIS_ID(B)\n",
    "    print B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A061356\n",
    "for n in (0..4):\n",
    "    B = inverse_bell_matrix(monotonicfactorial(n), 7)\n",
    "    print \"Inverse Bell matrix generated by monotonic factorial(r, n)\", n\n",
    "    print OEIS_ID(B)\n",
    "    print B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fallingfactorial = lambda n: lambda k: falling_factorial(n, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A023531, (A049403, A122848), A049404, A049410, A049424, A049411, ....  \n",
    "for n in (0..5):\n",
    "    B = bell_matrix(fallingfactorial(n), 7)\n",
    "    print \"Bell matrix generated by falling factorial(n, k)\", n\n",
    "    print OEIS_ID(B)\n",
    "    print B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for n in (0..7):\n",
    "    B = inverse_bell_matrix(fallingfactorial(n), 7)\n",
    "    print \"Inverse Bell matrix generated by falling factorial(n, k)\", n\n",
    "    print OEIS_ID(B)\n",
    "    print B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bell_second_order(generator, n):\n",
    "    G = [generator(k) for k in range(n)]\n",
    "    row = lambda n: bell_transform(n, G)\n",
    "    S = [sum(row(k)) for k in range(n)]\n",
    "    return bell_transform(n, S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A264430\n",
    "[bell_second_order(bell_number, n)  for n in range(8)] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A187761\n",
    "[sum(s) for s in [bell_second_order(lambda k: 1, n)  for n in range(10)]] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A265023\n",
    "[sum(s) for s in [bell_second_order(lambda k: (-1)^k, n)  for n in range(10)]] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Generalizes A000248\n",
    "[bell_second_order(lambda k: k+1, n)  for n in range(8)] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bell_transform_coeff(n, k, a):\n",
    "    @cached_function\n",
    "    def B(n, k):\n",
    "        if k==0: return a[0]^n\n",
    "        return sum(binomial(n-1,j-1)*a[j]*B(n-j,k-1) for j in (0..n-k+1))\n",
    "    return B(n, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def partial_bell_polynomial(n, k):\n",
    "    v = var(['x_'+str(i) for i in (0..n+1)])\n",
    "    return bell_transform_coeff(n, k, v).expand() \n",
    "\n",
    "for n in (0..4):\n",
    "    print [partial_bell_polynomial(n,k) for k in (0..n)] \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def bell_polynomial(n):\n",
    "    return sum(partial_bell_polynomial(n,k).subs(x_0=0) for k in (0..n))\n",
    "\n",
    "for n in (0..4): print bell_polynomial(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Alternatively:\n",
    "def bell_polynomial1(n):\n",
    "    X = var(['x_'+str(i) for i in (0..n+1)])\n",
    "    @cached_function\n",
    "    def B(n, k):\n",
    "        if k==0: return k^n\n",
    "        return sum(binomial(n-1,j-1)*B(n-j,k-1)*X[j] for j in (0..n-k+1)).expand()\n",
    "    return sum(B(n,k) for k in (0..n)) \n",
    "\n",
    "for n in (0..4): print bell_polynomial1(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def univariate_bell_polynomial(n):\n",
    "    p = bell_polynomial(n).subs(x_0=0)\n",
    "    q = p({p.variables()[i]:x for i in range(len(p.variables()))})\n",
    "    R = PolynomialRing(QQ, 'x')\n",
    "    return R(q)\n",
    "\n",
    "for n in (0..6): print univariate_bell_polynomial(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Alternatively:\n",
    "def univariate_bell_polynomial(n):\n",
    "    polynom = x^n\n",
    "    fn = factorial(n)\n",
    "    for k in (0..n-1):\n",
    "        result = 0 \n",
    "        for p in Partitions(n, length=k):\n",
    "            factorial_product = 1\n",
    "            power_factorial_product = 1\n",
    "            for part, count in p.to_exp_dict().iteritems():\n",
    "                factorial_product *= factorial(count)\n",
    "                power_factorial_product *= factorial(part)**count\n",
    "            coefficient = fn//(factorial_product*power_factorial_product)\n",
    "            result += coefficient  \n",
    "        polynom += result*x^k\n",
    "    return polynom\n",
    "\n",
    "for n in (0..6): print univariate_bell_polynomial(n)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  },
  "name": "BellTransform.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 0
}