2-ImplementationFactor.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation - Factors and operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Factor data structure\n",
    "\n",
    "Here factors are represented as multidimensional ```numpy``` arrays. Each axis corresponds to a particular random variable. Names of the variables are stored in ```self.__variables``` in the appropriate order.\n",
    "\n",
    "For example, if we have a factor ```psi```, the value of $(a_2, b_1, c_3)$ can be accessed as follows: ```psi.get_distribution()[1,0,2]```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class factor:\n",
    "    def __init__(self, variables = None, distribution = None):\n",
    "        if (distribution is None) and (variables is not None):\n",
    "            self.__set_data(np.array(variables), None, None)\n",
    "        elif (variables is None) or (len(variables) != len(distribution.shape)):\n",
    "            raise Exception('Data is incorrect')\n",
    "        else:\n",
    "            self.__set_data(np.array(variables),\n",
    "                            np.array(distribution),\n",
    "                            np.array(distribution.shape))\n",
    "    \n",
    "    def __set_data(self, variables, distribution, shape):\n",
    "        self.__variables    = variables\n",
    "        self.__distribution = distribution\n",
    "        self.__shape        = shape\n",
    "    \n",
    "    # ----------------------- Info --------------------------\n",
    "    def is_none(self):\n",
    "        return True if self.__distribution is None else False\n",
    "        \n",
    "    # ----------------------- Getters -----------------------\n",
    "    def get_variables(self):\n",
    "        return self.__variables\n",
    "    \n",
    "    def get_distribution(self):\n",
    "        return self.__distribution\n",
    "    \n",
    "    def get_shape(self):\n",
    "        return self.__shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Factor operations\n",
    "\n",
    "There are several important operations on factors:\n",
    "\n",
    "* factor product\n",
    "* factor marginalization\n",
    "* factor reduction\n",
    "* computing joint distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Factor product\n",
    "\n",
    "Representative example:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](./img/2/2.1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>Fig 2.1: Figure is from [1]</center>\n",
    "\n",
    "Here we are going to use ```numpy``` [broadcasting](https://docs.scipy.org/doc/numpy-1.17.0/user/basics.broadcasting.html) rules to perform tensor multiplication efficiently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def factor_product(x, y):\n",
    "    if x.is_none() or y.is_none():\n",
    "        raise Exception('One of the factors is None')\n",
    "    \n",
    "    xy, xy_in_x_ind, xy_in_y_ind = np.intersect1d(x.get_variables(), y.get_variables(), return_indices=True)\n",
    "    \n",
    "    if xy.size == 0:\n",
    "        raise Exception('Factors do not have common variables')\n",
    "    \n",
    "    if not np.all(x.get_shape()[xy_in_x_ind] == y.get_shape()[xy_in_y_ind]):\n",
    "        raise Exception('Common variables have different order')\n",
    "    \n",
    "    x_not_in_y = np.setdiff1d(x.get_variables(), y.get_variables(), assume_unique=True)\n",
    "    y_not_in_x = np.setdiff1d(y.get_variables(), x.get_variables(), assume_unique=True)\n",
    "    \n",
    "    x_mask = np.isin(x.get_variables(), xy, invert=True)\n",
    "    y_mask = np.isin(y.get_variables(), xy, invert=True)\n",
    "    \n",
    "    x_ind = np.array([-1]*len(x.get_variables()), dtype=int)\n",
    "    y_ind = np.array([-1]*len(y.get_variables()), dtype=int)\n",
    "    \n",
    "    x_ind[x_mask] = np.arange(np.sum(x_mask))\n",
    "    y_ind[y_mask] = np.arange(np.sum(y_mask)) + np.sum(np.invert(y_mask))\n",
    "    \n",
    "    x_ind[xy_in_x_ind] = np.arange(len(xy)) + np.sum(x_mask)\n",
    "    y_ind[xy_in_y_ind] = np.arange(len(xy))\n",
    "    \n",
    "    x_distribution = np.moveaxis(x.get_distribution(), range(len(x_ind)), x_ind)\n",
    "    y_distribution = np.moveaxis(y.get_distribution(), range(len(y_ind)), y_ind)\n",
    "                \n",
    "    res_distribution =   x_distribution[tuple([slice(None)]*len(x.get_variables())+[None]*len(y_not_in_x))] \\\n",
    "                       * y_distribution[tuple([None]*len(x_not_in_y)+[slice(None)])]\n",
    "    \n",
    "    return factor(list(x_not_in_y)+list(xy)+list(y_not_in_x), res_distribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example 1\n",
    "\n",
    "Let's reproduce the example above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_1 = factor(['a', 'b'], np.array([[0.5, 0.8], [0.1, 0.0], [0.3, 0.9]]))\n",
    "phi_2 = factor(['b', 'c'], np.array([[0.5, 0.7], [0.1, 0.2]]))\n",
    "\n",
    "phi_3 = factor_product(phi_1, phi_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['a', 'b', 'c'], dtype='<U1')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi_3.get_variables()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[0.25, 0.35],\n",
       "        [0.08, 0.16]],\n",
       "\n",
       "       [[0.05, 0.07],\n",
       "        [0.  , 0.  ]],\n",
       "\n",
       "       [[0.15, 0.21],\n",
       "        [0.09, 0.18]]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi_3.get_distribution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example 2\n",
    "\n",
    "Conditional probability distribution (CPD) and joint distribution can be represented as factors. Suppose we have $P(a|b)$ and $P(b)$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|        |  $$P(a_1| b)$$   |  $$P(a_2|b)$$   |  $$P(a_3|b)$$   |\n",
    "|---     |---               |---              |---              |\n",
    "| $b_1$  |     0.3          |      0.2        |      0.5        |\n",
    "| $b_2$  |     0.8          |      0.1        |      0.1        |\n",
    "\n",
    "<center>Fig 2.2: CPD $P(a|b)$</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|        |  $$P(b)$$   |  \n",
    "|---     |---          |\n",
    "| $b_1$  |     0.3     |\n",
    "| $b_2$  |     0.7     |\n",
    "\n",
    "<center>Fig 2.3: P(b)</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute $P(a,b) = P(a|b)P(b)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|        |  $a_1$    |  $a_2$    |  $a_3$   |\n",
    "|---     |---        |---        |---       |\n",
    "| $b_1$  |   0.09    |    0.06   |    0.15  |\n",
    "| $b_2$  |   0.56    |    0.07   |    0.7   |\n",
    "\n",
    "<center>Fig 2.4: Joint distribution P(a, b)</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['a', 'b', 'c'], dtype='<U1')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi_4 = factor(['a', 'b'], np.array([[0.3, 0.8], [0.2, 0.1], [0.5, 0.1]]))\n",
    "phi_5 = factor(['b'], np.array([0.3, 0.7]))\n",
    "\n",
    "phi_6 = factor_product(phi_1, phi_2)\n",
    "phi_6.get_variables()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[0.25, 0.35],\n",
       "        [0.08, 0.16]],\n",
       "\n",
       "       [[0.05, 0.07],\n",
       "        [0.  , 0.  ]],\n",
       "\n",
       "       [[0.15, 0.21],\n",
       "        [0.09, 0.18]]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi_6.get_distribution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Factor marginalization\n",
    "\n",
    "Representative example:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](./img/2/2.2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>Fig 2.5: Figure is from [1]</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def factor_marginalization(x, variables):\n",
    "    variables = np.array(variables)\n",
    "    \n",
    "    if x.is_none():\n",
    "        raise Exception('Factor is None')\n",
    "    \n",
    "    if not np.all(np.in1d(variables, x.get_variables())):\n",
    "        raise Exception('Factor do not contain given variables')\n",
    "    \n",
    "    res_variables    = np.setdiff1d(x.get_variables(), variables, assume_unique=True)\n",
    "    res_distribution = np.sum(x.get_distribution(),\n",
    "                              tuple(np.where(np.isin(x.get_variables(), variables))[0]))\n",
    "    \n",
    "    return factor(res_variables, res_distribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example\n",
    "\n",
    "Let's reproduce the example above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_4 = factor_marginalization(phi_3, ['b'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['a', 'c'], dtype='<U1')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi_4.get_variables()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Factor reduction\n",
    "\n",
    "Representative example:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](./img/2/2.3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>Fig 2.6: Figure is from [1]</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def factor_reduction(x, variable, value):\n",
    "    if x.is_none() or (variable is None) or (value is None):\n",
    "        raise Exception('Input is None')\n",
    "    \n",
    "    if not np.any(variable == x.get_variables()):\n",
    "        raise Exception('Factor do not contain given variable')\n",
    "    \n",
    "    if value >= x.get_shape()[np.where(variable==x.get_variables())[0]]:\n",
    "        raise Exception('Incorrect value of given variable')\n",
    "    \n",
    "    res_variables    = np.setdiff1d(x.get_variables(), variable, assume_unique=True)\n",
    "    res_distribution = np.take(x.get_distribution(),\n",
    "                               value,\n",
    "                               int(np.where(variable==x.get_variables())[0]))\n",
    "    \n",
    "    return factor(res_variables, res_distribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example\n",
    "\n",
    "Let's reproduce the example above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['a', 'b'], dtype='<U1')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi_7 = factor_reduction(phi_3, 'c', 0)\n",
    "phi_7.get_variables()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.25, 0.08],\n",
       "       [0.05, 0.  ],\n",
       "       [0.15, 0.09]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi_7.get_distribution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Joint distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def joint_distribution(ar):\n",
    "    for element in ar:\n",
    "        if element.is_none():\n",
    "            raise Exception('Factor is None')\n",
    "    \n",
    "    res = ar[0]\n",
    "    for element in ar[1:]:\n",
    "        res = factor_product(res, element)\n",
    "    \n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*References:*\n",
    "\n",
    "1. D Koller \"Probabilistic Graphical Models\" by Stanford at Coursera [[link](https://ru.coursera.org/lecture/probabilistic-graphical-models/factors-tEZ6S)]"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}