alg/symbolic_metamodeling/1-_Introduction_to_Meijer_G-functions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Meijer $G$-functions\n",
    "\n",
    "This notebook provides introductory material on the Meijer-$G$ functions based on our paper *\"Demystifying Black-box Models with Symbolic Metamodels\"* submitted to **NeurIPS 2019** by *Ahmed M. Alaa and Mihaela van der Schaar*. The notebook provides the formal definition of Meijer-$G$ functions, along with examples for how to use our Python implementation of the $G$-function in the context of symbolic computation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What are Meijer $G$-functions?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **Meijer G-functions** are a family of univariate functions, each of them determined by finitely many indexes. The $G$-function is a linear combination of certain special functions of standard type, and hence it is often use to solve indefinite integrals [1]. Despite their useful properties, $G$-functions are not very well known in the mathematical community. In this notebook, we provide a brief background on Meijer $G$-functions, presenting its formal definition along with its basic properties. In addition, we show how to use our python implementation of $G$-functions in symbolic computations. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formally, the Meijer $G$-function is a univariate function given by a line integral in the complex plane [2, 3]:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$G^{m,n}_{p, q}\\left(^{a_1,\\ldots,a_p}_{b_1,\\ldots,b_q}\\,\\Big|\\, x\\right) = \\frac{1}{2\\pi i} \\int_{\\mathcal{L}} \\frac{\\prod^m_{j=1} \\Gamma(b_j - s)\\prod^n_{j=1} \\Gamma(1 - a_j + s)}{\\prod^q_{j=m+1} \\Gamma(1 - b_j + s)\\prod^p_{j=n+1} \\Gamma(a_j + s)}\\, x^s\\, ds$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $\\Gamma(.)$ is the Gamma function, $0 \\leq m \\leq q$, $0 \\leq n \\leq p$, where $m$, $n$, $p$ and $q$ are integer numbers, $a_k - b_j \\neq 1, 2, 3,.\\,.\\,.$ for $k = 1, 2,\\,.\\,.\\,., n$, $j = 1, 2\\,.\\,.\\,., m$ and $x \\neq 0$. This means that no pole of any $\\Gamma(b_j - s)$ coincides with any pole of any $\\Gamma(1 - a_k + s)$. The integration path $\\mathcal{L}$ separates the poles of the factors $\\Gamma(b_j - s)$ from those of the factors $\\Gamma(1 - a_k + s)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/FigA1.png\", width=640>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown in the Figure above, there are three cases in which integration path $\\mathcal{L}$ can be chosen to separate the poles and the zeros as follows [4]:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- $L$ goes from $-i\\infty$ to $i\\infty$.\n",
    "- $L$ is a loop that starts at $i\\infty$ on a line parallel to the $+$ve real axis, encircles the poles of the factors $\\Gamma(b_j - s)$ and returns to $i\\infty$ on another line parallel to the $+$ve real axis.\n",
    "- $L$ is a loop that starts at $i\\infty$ on a line parallel to the $-$ve real axis, encircles the poles of the factors $\\Gamma(1 - a_k + s)$ and returns to $\\infty$ on another line parallel to the $-$ve real axis.\n",
    "\n",
    "When more than one of Cases 1, 2, and 3 is applicable the same value is obtained for the $G$-function. The conditions for convergence of the line integral can be established by applying Stirling's asymptotic approximation to the gamma functions $\\Gamma(b_j - s)$ and $\\Gamma(1 - a_k + s)$  in the integrand. As a consequence of the definition above, the Meijer $G$-function is an analytic function of $x$ (with possible exception of the origin $x$ = 0 and the unit circle $|x| = 1$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to use Meijer $G$-functions in PySymbolic?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**PySymbolic** is a symbolic computation library that we developed in order to implement symbolic metamodeling. This is basically a wrapper over **SymPy** that provides all the relevant functionalities needed for metamodeling, in addition to a high-level API that makes symbolic comptations easier. \n",
    "\n",
    "To see how can we evaluate Meijer $G$-functions in PySymbolic, we first import the **MeijerG** class from the **models.special_functions** module as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pysymbolic.models.special_functions import MeijerG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to create an instance of a Meijer $G$-function, we only need to provide it with the poles and the zeros through two input arguments as follows:\n",
    "    \n",
    "**MeijerG(theta=[a_1,..., a_n,..., a_p, b_1,...,b_m,...,b_q, scale], order=[m, n, p, q])**\n",
    "\n",
    "Here the arguments **theta** and **order** correspond to the actual values of the poles and zeros, and the indexes $m$, $n$, $p$ and $q$, respectively. The **scale** variable corresponds to a real number by which the input $x$ is scaled. To demonstrate the above, let us create an object holding the following Meijer $G$-function: \n",
    "\n",
    "$G^{0,1}_{3, 1}\\left(^{2, 2, 2}_{1}\\,\\Big|\\, 5x\\right)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create the corresponding MeijerG object, we define the inputs as theta = [2,2,2,1,5] since $a_1=a_2=a_3 = 2$, and $b_1=1$, whereas the scaling factor is 5. Thus, we have the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MeijerG(theta=[2, 2, 2, 1, 5], order=[0, 1, 3, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now to recover the expression above in a symbolic format, we use the \"expression\" method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "My_meijerg_object = MeijerG(theta=[2, 2, 2, 1, 5], order=[0, 1, 3, 1])\n",
    "\n",
    "My_meijerg_object.expression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see above, it turns out that this particular selection of poles and zeros reduces the Meijer $G$-function to an identity function! This is the power of Meijer $G$-functions: for different selections of poles and zeros, they reduce to familiar functional forms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the **theta** and **order** inputs, the **MeijerG** object takes two other optional arguments:\n",
    "\n",
    "- **evaluation_mode**: When evaluating the value numerical value of a Meijer $G$-function at a specific value of $x$, this argument determines the way numerical evaluation is carried out. The options include ['eval','numpy', 'theano', 'cython']. The 'eval' option applies direct symbolic computation, which is the slowest approach. 'numpy' and 'theano' apply numeric evaluation using numpy arrays or theano graphs, whereas 'cython' creates an intermediate C code to enable fast computation of the Meijer $G$-function. The default setting is 'numpy'.  \n",
    "- **approximation_order**: When approximating the Meijer $G$-function for intermediate computation, we use a Taylor series approximation. This argument determines the number of terms in a Taylor expansion to include in the approximation, and is set to a default value of 15.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned earlier, the main strength of the Meijer $G$-function is that they reduce to familiar functional forms for different selections of poles and zeros. Notable examples include the following:\n",
    "\n",
    "- $G^{0,1}_{3, 1}\\left(^{2, 2, 2}_{1}\\,\\Big|\\, x\\right) = x$ (which we have demonstrated above).\n",
    "\n",
    "\n",
    "- $G^{1,0}_{0, 1}\\left(^{-}_{0}\\,\\Big|\\, x\\right) = e^{-x}$.\n",
    "\n",
    "\n",
    "- $G^{1,2}_{2, 2}\\left(^{1,1}_{1,0}\\,\\Big|\\, x\\right) = \\log(x+1)$.\n",
    "\n",
    "\n",
    "- $G^{1,2}_{2, 2}\\left(^{1,1}_{1,1}\\,\\Big|\\, x\\right) = \\frac{x}{(x+1)}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In what follows, we show how can we recover all the identities above using the MeijerG object. To enable printing in LaTex format inside the notebook, we use the \"init_print\" command in sympy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "identity_function    = MeijerG(theta=[2, 2, 2, 1, 1], order=[0, 1, 3, 1])\n",
    "exponential_function = MeijerG(theta=[0, 1], order=[1, 0, 0, 1])\n",
    "logarithmic_function = MeijerG(theta=[1,1,1,0,1], order=[1, 2, 2, 2])\n",
    "rational_function    = MeijerG(theta=[1,1,1,1,1], order=[1, 2, 2, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us check if the magic works!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "identity_function.expression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exponential_function.expression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logarithmic_function.expression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rational_function.expression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "[1] Richard Beals and Jacek Szmigielski. Meijer g-functions: a gentle introduction. *Notices of the AMS*, 60(7):866–872, 2013.\n",
    "\n",
    "[2] CS Meijer. On the G-function. *North-Holland*, 1946.\n",
    "\n",
    "[3] CS Meijer. Uber whittakersche bezw. besselsche funktionen und deren produkte (english translation: About whittaker and bessel functions and their products). *Nieuw Archief voor Wiskunde*, 18(2):10–29, 1936.\n",
    "\n",
    "[4] I. Gradshteyn and I. Ryzhik. Table of integrals, series, and products. *Academic press*, 2014."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}