01-feature-selection.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Feature Selection with the D-Wave System\n",
    "This notebook demonstrates the formulation of an $n \\choose k$ optimization problem for solution on a D-Wave quantum computer. The method used in the example problem of this notebook&mdash;feature-selection for machine learning&mdash; is applicable to problems from a wide range of domains; for example financial portfolio optimization. \n",
    "    \n",
    "1. [What is Feature Selection?](#What-is-Feature-Selection?) defines and explains the feature-selection problem.\n",
    "2. [Feature Selection by Mutual Information](#Feature-Selection-by-Mutual-Information) describes a particular method of feature selection that is demonstrated in this notebook.\n",
    "3. [Solving Feature Selection on a Quantum Computer](#Solving-Feature-Selection-on-a-Quantum-Computer) shows how such optimization problems can be formulated for solution on a D-Wave quantum computer. \n",
    "4. [Example Application: Predicting Survival of Titanic Passengers](#Example-Application:-Predicting-Survival-of-Titanic-Passengers) demonstrates the use of *Kerberos*, an out-of-the-box classical-quantum [hybrid](https://docs.ocean.dwavesys.com/en/stable/docs_hybrid/sdk_index.html) sampler, to select optimal features for a public-domain dataset. \n",
    "\n",
    "This notebook should help you understand both the techniques and [Ocean software](https://docs.ocean.dwavesys.com) tools used for solving optimization problems on D-Wave quantum computers."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**New to Jupyter Notebooks?** JNs are divided into text or code cells. Pressing the **Run** button in the menu bar moves to the next cell. Code cells are marked by an \"In: \\[\\]\" to the left; when run, an asterisk displays until code completion: \"In: \\[\\*\\]\"."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What is Feature Selection?\n",
    "Statistical and machine-learning models use sets of input variables (\"features\") to predict output variables of interest. Feature selection can be part of the model design process: selecting from a large set of potential features a highly informative subset simplifies the model and reduces dimensionality. \n",
    "\n",
    "For example, to build a model that predicts the ripening of hothouse tomatoes, Farmer MacDonald daily records the date, noontime temperature, daylight hours, degree of cloudiness, rationed water and fertilizer, soil humidity, electric-light hours, etc. These measurements constitute a list of potential features. After a growth cycle or two, her analysis reveals some correlations between these features and crop yields:\n",
    "\n",
    "* fertilizer seems a strong predictor of fruit size\n",
    "* cloudiness and daylight hours seem poor predictors of growth \n",
    "* water rations and soil humidity seem a highly correlated pair of strong predictors of crop rot   \n",
    "\n",
    "Farmer MacDonald suspects that her hothouse's use of electric light reduces dependency on seasons and sunlight. She can simplify her model by discarding date, daylight hours, and cloudiness. She can record just water ration or just soil humidity rather than both.\n",
    "\n",
    "For systems with large numbers of potential input information&mdash;for example, weather forecasting or image recognition&mdash;model complexity and required compute resources can be daunting. Feature selection can help make such models tractable. \n",
    "\n",
    "However, optimal feature selection can itself be a hard problem. This example introduces a powerful method of optimizing feature selection based on a complex probability calculation. This calculation is submitted for solution to a quantum computer. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Illustrative Toy Problem\n",
    "This subsection illustrates the use of feature selection with a simple example: a toy system with a single output generated from three inputs. \n",
    "\n",
    "The model built to predict the system's output is even simpler: it uses just two of three possible features (inputs). You can expect it to perform better when the selected two features are more independent, assuming all three contribute somewhat commensurately to the system output (if an independent feature contributes less than the difference between two dependent ones, this might not be true). In the case of Farmer MacDonalds's tomatoes, a model using rationed water and fertilizer should perform better than one using rationed water and soil humidity. \n",
    "\n",
    "The code cell below uses the NumPy library to define three inputs, the first two of which are very similar: a sine, a noisy sine, and a linear function with added random noise. It defines an output that is a simple linear combination of these three inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "sig_len = 100\n",
    "# Three inputs: in1 & in2 are similar \n",
    "in1 = np.sin(np.linspace(-np.pi, np.pi, sig_len)).reshape(sig_len, 1)\n",
    "in2 = np.sin(np.linspace(-np.pi+0.1, np.pi+0.2, sig_len)).reshape(sig_len, 1) + 0.3*np.random.rand(sig_len, 1)\n",
    "in3 = np.linspace(-1, 1, sig_len).reshape(sig_len,1) + 2*np.random.rand(sig_len, 1)\n",
    "\n",
    "out = 2*in1 + 3*in2 + 6*in3"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the features and variable of interest (the output). In this and other cells below, graphics code is imported from a `helpers` module. To see this code, navigate to the `helpers` folder of this repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers.plots import plot_toy_signals\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "\n",
    "# Store problem in a pandas DataFrame for later use\n",
    "toy = pd.DataFrame(np.hstack((in1, in2, in3, out)), columns=[\"in1\", \"in2\", \"in3\", \"out\"])\n",
    "\n",
    "plot_toy_signals(toy)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `two_var_model` function in the next cell defines a linear model of two variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def two_var_model(in_tuple, a, b):\n",
    "    ina, inb = in_tuple\n",
    "    return a*ina + b*inb"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the SciPy library's `curve_fit` function to try to predict the variable of interest from two of three features. The model with less-correlated features performs better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers.plots import plot_two_var_model \n",
    "import itertools\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "model = []\n",
    "two_vars = []\n",
    "for f0, f1 in itertools.combinations(['in1', 'in2', 'in3'], 2):  \n",
    "    two_vars.append((f0, f1))\n",
    "    popt, pcov = curve_fit(two_var_model, (toy[f0].values, toy[f1].values), toy['out'].values)\n",
    "    model.append(two_var_model((toy[f0].values, toy[f1].values), popt[0], popt[1]).reshape(len(toy), 1))\n",
    "    print(\"Standard deviation for selection of features {} and {} is {:.4f}.\".format(f0, f1, max(np.sqrt(np.diag(pcov)))))\n",
    "model_df = pd.DataFrame(np.hstack(model), columns=two_vars)\n",
    "\n",
    "plot_two_var_model(model_df, toy)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following section introduces a method for optimizing feature selection that can be useful in modeling complex systems. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Feature Selection by Mutual Information\n",
    "There are various methods for [feature selection](https://en.wikipedia.org/wiki/Feature_selection). If you are building a machine-learning model, for example, and have six potential features, you might naively consider training it first on each of the features by itself, then on all 15 combinations of subsets of two features, then 20 combinations of subsets of three features, and so on. However, statistical methods are more efficient.  \n",
    "\n",
    "One statistical criterion that can guide this selection is mutual information (MI). The following subsections explain information and MI with some simple examples.\n",
    "\n",
    "If you already understand MI and Shannon entropy, please skip ahead to section [Solving Feature Selection on a Quantum Computer](#Solving-Feature-Selection-on-a-Quantum-Computer) (and then from the menu bar, click the *Run All Above* option). "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantifying Information: Shannon Entropy\n",
    "[Shannon entropy](https://en.wiktionary.org/wiki/Shannon_entropy), $H(X)$,  mathematically quantifies the information in a signal:\n",
    "\n",
    "$H(X) = - \\sum_{x \\in X} p(x) \\log p(x)$\n",
    "\n",
    "where $p(x)$ represents the probability of an event's occurrence. The Shannon Entropy (SE) formula can be understood as weighing by an event's probability a value of $\\log \\frac{1}{p(x)}$ for the event, where the reciprocal is due to the minus sign. This value means that the less likely the occurrence of an event, the more information is attributed to it (intuitively, when a man bites a dog it's news). \n",
    "\n",
    "To calculate SE, the `prob` function defined below calculates probability for a dataset representing some variables (a training set in a machine learning context) by dividing it into bins as a histogram using the NumPy library's `histogramdd` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prob(dataset, max_bins=10):\n",
    "    \"\"\"Joint probability distribution P(X) for the given data.\"\"\"\n",
    "\n",
    "    # bin by the number of different values per feature\n",
    "    num_rows, num_columns = dataset.shape\n",
    "    bins = [min(len(np.unique(dataset[:, ci])), max_bins) for ci in range(num_columns)]\n",
    "\n",
    "    freq, _ = np.histogramdd(dataset, bins)\n",
    "    p = freq / np.sum(freq)\n",
    "    return p\n",
    "\n",
    "def shannon_entropy(p):\n",
    "    \"\"\"Shannon entropy H(X) is the sum of P(X)log(P(X)) for probabilty distribution P(X).\"\"\"\n",
    "    p = p.flatten()\n",
    "    return -sum(pi*np.log2(pi) for pi in p if pi)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Illustration of Shannon Entropy\n",
    "For an intuitive example of measuring SE, this subsection applies the `shannon_entropy` function to three signals with defined distributions:\n",
    "\n",
    "* Uniform: this distribution maximizes values of SE because all outcomes are equally likely, meaning every outcome is equally unpredictable. $H(X) = log(N)$ for uniform distribution, where $N$ is the number of possible outcomes, ${x_1, x_2, ...x_N}$. \n",
    "* Exponential: the steeper the curve, the more outcomes are in the \"tail\" part (have higher probability) with lower information value. \n",
    "* Binomial: the stronger this signal is biased to one outcome, the more predictable its values, the lower its information value. $H(X) = -p \\log(p) - (1-p) \\log(1-p)$ for binomial distribution; for $p = 0.1$, for example, $H(X) = 0.468$.  \n",
    "\n",
    "Define the three signals and plot the SE. The red dots show the maximal values of SE for different numbers of bits (Shannon developed the formula to calculate channel bandwidth, which for digital communications is measured in bits) or, as here, the bins into which the signals' possible values are divided.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers.plots import plot_se \n",
    "\n",
    "max_bins = 10\n",
    "\n",
    "# Signals with defined distributions\n",
    "x_uniform = np.random.uniform(0, 1, (1000, 1))\n",
    "x_exp = np.exp(-np.linspace(0, 10, 1000)/2).reshape(1000, 1)\n",
    "x_vals = np.random.choice([0, 1],(1000, 1), p=[0.1, 0.9])\n",
    "\n",
    "data = list()\n",
    "for bins in range(1, max_bins):\n",
    "    uniform_se = shannon_entropy(prob(x_uniform, bins))\n",
    "    exp_se = shannon_entropy(prob(x_exp, bins))\n",
    "    vals_se = shannon_entropy(prob(x_vals, bins))                               \n",
    "    data.append({'Bins': bins, 'Uniform': uniform_se, 'Maximum': np.log2(bins), 'Exp': exp_se, 'Vals': vals_se})\n",
    "\n",
    "plot_se(data)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conditional Shannon Entropy\n",
    "\n",
    "Conditional SE (CSE) measures the information in one signal, $X$, when the value of another signal, $Y$, is known: \n",
    "\n",
    "$\\begin{aligned} H(X|Y) \n",
    "& = H(X,Y)-H(Y) \\\\\n",
    "& = - \\sum_{x \\in X} p(x, y) \\log p(x, y) - H(Y) \\end{aligned}$\n",
    "\n",
    "where joint SE, $H(X,Y)$, measures the information in both signals together, with $p(x,y)$ being their joint probability. For example, knowing that it's winter reduces the information value of news that it is raining.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conditional_shannon_entropy(p, *conditional_indices):\n",
    "    \"\"\"Shannon entropy of P(X) conditional on variable j\"\"\"\n",
    "\n",
    "    axis = tuple(i for i in np.arange(len(p.shape)) if i not in conditional_indices)\n",
    "\n",
    "    return shannon_entropy(p) - shannon_entropy(np.sum(p, axis=axis))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Illustration of CSE \n",
    "Apply CSE to the toy problem. Because signals `in1` and `in2` are similar, knowing the value of one provides a good estimate of the other; in contrast, the value of signal `in3` is less good for estimating the first two.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"H(in1) = {:.2f}\".format(shannon_entropy(prob(toy[[\"in1\"]].values))))\n",
    "print(\"H(in1|in3) = {:.2f}\".format(conditional_shannon_entropy(prob(toy[[\"in1\", \"in3\"]].values), 1)))\n",
    "print(\"H(in1|in2) = {:.2f}\".format(conditional_shannon_entropy(prob(toy[[\"in1\", \"in2\"]].values), 1)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mutual Information\n",
    "[Mutual information](https://en.wikipedia.org/wiki/Mutual_information) between variables $X$ and $Y$ is defined as \n",
    "\n",
    "$I(X;Y)  = \\sum_{y \\in Y} \\sum_{x \\in X} p(x, y) \\log \\frac{p(x,y)}{p(x)p(y)}$\n",
    "\n",
    "where $p(x)$ and $p(y)$ are marginal probabilities of $X$ and $Y$, and $p(x,y)$ the joint probability. Equivalently, \n",
    "\n",
    "$I(X;Y)  = H(Y) - H(Y|X)$\n",
    "\n",
    "where $H(Y)$ is the SE of $Y$ and $H(Y|X)$ is the CSE of $Y$ conditional on $X$.\n",
    "\n",
    "Mutual information (MI) quantifies how much one knows about one random variable from observations of another. Intuitively, a model based on just one of a pair of features (e.g., farmer MacDonald's water rations and soil humidity) will better reproduce their combined contribution when MI between them is high. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mutual_information(p, j):\n",
    "    \"\"\"Mutual information between all variables and variable j\"\"\"\n",
    "    return shannon_entropy(np.sum(p, axis=j)) - conditional_shannon_entropy(p, j)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mutual Information on the Toy Problem\n",
    "Calculate and plot MI between the output of the toy problem and its three input signals. This measures the suitability of each on its own as a feature in a model of the system, or how much each shapes the output.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from helpers.plots import plot_mi \n",
    "\n",
    "mi = {}\n",
    "for column in toy.columns:\n",
    "    if column == 'out':\n",
    "        continue\n",
    "    mi[column] = mutual_information(prob(toy[['out', column]].values), 1)\n",
    "\n",
    "plot_mi(mi)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot of input and output signals in the first section might give an impression that the toy model's output is closer to the two sine signals than to `in3`, but the linear regression below confirms the MI result. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import linregress\n",
    "from helpers.plots import plot_lingress \n",
    "\n",
    "model = []\n",
    "var_rval = []\n",
    "for column in toy.columns:\n",
    "    if column == 'out':\n",
    "        continue\n",
    "    slope, intercept, rvalue, pvalue, stderr = linregress(toy[column].values, toy['out'].values)  \n",
    "    model.append((slope*toy[column].values + intercept).reshape(len(toy), 1))\n",
    "    var_rval.append((column, rvalue))\n",
    "\n",
    "plot_lingress(pd.DataFrame(np.hstack(model), columns=var_rval), toy)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result should in fact be expected given an output defined as $out = 2 \\times in_1 + 3 \\times in_2 + 6 \\times in_3$ with the sixfold multiplier on $in_3$ greater than the sum of multipliers on the other signals, all three of which have an amplitude of 1.  "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conditional Mutual Information\n",
    "\n",
    "Conditional mutual information (CMI) between a variable of interest, $X$, and a feature, $Y$, given the selection of another feature, $Z$, is given by\n",
    "\n",
    "$I(X;Y|Z) = H(X|Z)-H(X|Y,Z)$\n",
    "\n",
    "where $H(X|Z)$ is the CSE of $X$ conditional on $Z$ and $H(X|Y, Z)$ is the CSE of $X$ conditional on both $Y$ and $Z$.\n",
    "\n",
    "In this code cell, because marginalization over $j$ removes a dimension, any conditional indices pointing to subsequent dimensions are decremented by 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conditional_mutual_information(p, j, *conditional_indices):\n",
    "    \"\"\"Mutual information between variables X and variable Y conditional on variable Z.\"\"\"\n",
    "\n",
    "    marginal_conditional_indices = [i-1 if i > j else i for i in conditional_indices]\n",
    "\n",
    "    return (conditional_shannon_entropy(np.sum(p, axis=j), *marginal_conditional_indices)\n",
    "            - conditional_shannon_entropy(p, j, *conditional_indices))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply `conditional_mutual_information` to the toy problem to find CMI between `out` and either `in2` or `in3`, conditional on `in1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"I(out;in2|in1) = {:.2f}\".format(conditional_mutual_information(prob(toy[['out', 'in2', 'in1']].values), 1, 2)))\n",
    "print(\"I(out;in3|in1) = {:.2f}\".format(conditional_mutual_information(prob(toy[['out', 'in3', 'in1']].values), 1, 2)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the sine signal, `in1`, if you try to predict the output from one of the remaining signals, you find that the linear function, `in3` contributes more information than the noisy sine, `in2`.\n",
    "\n",
    "Ideally, to select a model's $k$ most relevant of $n$ features, you could maximize $I({X_k}; Y)$, the MI between a set of $k$ features, $X_k$, and variable of interest, $Y$. This is a hard calculation because $n \\choose k$ grows rapidly in real-world problems."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving Feature Selection on a Quantum Computer\n",
    "There are different methods of approximating the hard calculation of optimally selecting $n \\choose k$ features to maximize MI. The approach followed here assumes conditional independence of features and limits CMI calculations to permutations of three features. The optimal set of features, $S$, is then approximated by:\n",
    "\n",
    "$\\arg \\max_S \\sum_{i \\in S} \\left \\{ I(X_i;Y) + \\sum_{j \\in S, j \\ne i} I(X_j;Y |X_i) \\right \\}$\n",
    "\n",
    "\n",
    "The left-hand component, $I(X_i;Y)$, represents MI between the variable of interest and a particular feature; maximizing selects features that best predict the variable of interest. The right-hand component, $I(X_j;Y |X_i)$, represents conditional MI between the variable of interest and a feature given the prior selection of another feature; maximizing selects features that complement information about the variable of interest rather than provide redundant information.\n",
    "\n",
    "This approximation is still a hard calculation. The following subsection demonstrates a method for formulating it for solution on the D-Wave quantum computer. The method is based on the 2014 paper, [Effective Global Approaches for Mutual Information Based Feature Selection](https://dl.acm.org/citation.cfm?id=2623611), by Nguyen, Chan, Romano, and Bailey published in the Proceedings of the 20th ACM SIGKDD international conference on knowledge discovery and data mining."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MIQUBO: QUBO Representation of Feature Selection\n",
    "D-Wave systems solve binary quadratic models (BQM)&mdash;the Ising model traditionally used in statistical mechanics and its computer-science equivalent, the quadratic unconstrained binary optimization (QUBO) problem. Given $N$ variables $x_1,...,x_N$, where each variable $x_i$ can have binary values $0$ or $1$, the system finds assignments of values that minimize,\n",
    "    \n",
    "$\\sum_i^N q_ix_i + \\sum_{i<j}^N q_{i,j}x_i  x_j$,\n",
    "    \n",
    "where $q_i$ and $q_{i,j}$ are configurable (linear and quadratic) coefficients. To formulate a problem for the D-Wave system is to program $q_i$ and $q_{i,j}$ so that assignments of $x_1,...,x_N$ also represent solutions to the problem.\n",
    "\n",
    "For feature selection, the Mutual Information QUBO (MIQUBO) method formulates a QUBO based on the approximation above for $I({X_k}; Y)$, which can be submitted to the D-Wave quantum computer for solution.\n",
    "\n",
    "The reduction of scope to permutations of three variables in this approximate formulation for MI-based optimal feature selection makes it a natural fit for reformulation as a QUBO: \n",
    "\n",
    "|  | **Optimization** &nbsp; &nbsp; | **Linear <br>Terms** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; |\t**Quadratic <br>Terms** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;| **Formula** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; |\t\n",
    "|:-|:-----------------|:---------------------|:-------------------------|:------------|\n",
    "| **Feature <br>Selection** &nbsp; &nbsp; | Maximize &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; |$I(X_i;Y)$  | $I(X_j;Y \\|X_i)$  | $\\sum_{i \\in S} \\left \\{ I(X_i;Y) + \\sum_{j \\in S, j \\ne i} I(X_j;Y \\|X_i) \\right \\}$ |\n",
    "| **QUBO** | Minimize &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;| $q_ix_i$ | $q_{i,j}x_ix_j$ | $\\sum_i^N q_ix_i + \\sum_{i<j}^N q_{i,j}x_i  x_j$ | \n",
    "\n",
    "You can represent each choice of $n \\choose k$ features as the value of solution $x_1,...,x_N$ by encoding $x_i=1$ if feature $X_i$ should be selected and $x_i=0$ if not. With solutions encoded this way, you can represent the QUBO in matrix format, $\\mathbf{x}^T \\mathbf{Q x}$, where $\\mathbf Q$ is an $n$ x $n$ matrix and $\\mathbf{x}$ is an $n$ x $1$ matrix (a vector) that should have $k$ ones representing the selected features. \n",
    "\n",
    "To map the feature-selection formula to a QUBO, set the elements of $\\mathbf Q$ such that\n",
    "\n",
    " * diagonal elements (linear coefficients) represent MI: $Q_{ii} \\leftarrow -I(X_i;Y)$ \n",
    " * non-diagonal elements (quadratic elements) represent CMI: $Q_{ij} \\leftarrow -I(X_j;Y |X_i)$\n",
    "\n",
    "These QUBO terms are negative because the quantum computer seeks to minimize the programmed problem while the feature-selection formula maximizes. The following subsection codes this and then completes the formulation by adding the $n \\choose k$ constraint to the QUBO.  "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MIQUBO on the Toy Problem\n",
    "This subsection applies the MIQUBO formulation to the toy problem by configuring the QUBO in three parts: (1) linear biases that maximize MI between the variable of interest and each feature (2) quadratic biases that maximize CMI between the variable of interest and each feature, given the prior choice of another feature (3) selection of just $k$ features.\n",
    "\n",
    "Create a BQM and set the linear coefficients as the MI between `out` and each potential feature. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dimod\n",
    "\n",
    "bqm = dimod.BinaryQuadraticModel.empty(dimod.BINARY)\n",
    "\n",
    "for column in toy.columns:\n",
    "\n",
    "    if column == 'out':\n",
    "        continue\n",
    "\n",
    "    mi = mutual_information(prob(toy[['out', column]].values), 1)\n",
    "    bqm.add_variable(column, -mi)\n",
    "\n",
    "for item in bqm.linear.items():\n",
    "    print(\"{}: {:.3f}\".format(item[0], item[1]))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the quadratic coefficients as the MI between `out` and each potential feature conditional on the other features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for f0, f1 in itertools.combinations(['in1', 'in2', 'in3'], 2):\n",
    "    cmi_01 = conditional_mutual_information(prob(toy[['out', f0, f1]].values), 1, 2)\n",
    "    cmi_10 = conditional_mutual_information(prob(toy[['out', f1, f0]].values), 1, 2)\n",
    "    bqm.add_interaction(f0, f1, -cmi_01)\n",
    "    bqm.add_interaction(f1, f0, -cmi_10)\n",
    "\n",
    "bqm.normalize()     # scale the BQM to (-1, 1) biases\n",
    "\n",
    "for item in bqm.quadratic.items():\n",
    "    print(\"{}: {:.3f}\".format(item[0], item[1]))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use Ocean software's [dimod](https://docs.ocean.dwavesys.com/en/stable/docs_dimod/sdk_index.html) [ExactSolver()](https://docs.ocean.dwavesys.com/en/stable/docs_dimod/reference/sampler_composites/samplers.html) to find an exact solution of the BQM, which currently represents a minimization of MI and CMI, and plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers.plots import plot_solutions \n",
    "\n",
    "sampler = dimod.ExactSolver()\n",
    "\n",
    "result = sampler.sample(bqm)\n",
    "\n",
    "plot_solutions(result)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unsurprisingly, the best solution (lowest-energy solution for the minimized QUBO) employs all three input signals because the current QUBO mapping does not constrain the number of selected features. In those solutions where only two are selected, models that select `in3` are better than the one that selects just `in1` and `in2`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Penalizing Non-k Selections\n",
    "How do you program  on the quantum computer a constraint that exactly $k$ features be selected? By penalizing solutions that select greater or fewer than $k$ features. If you add \n",
    "\n",
    "$P = \\alpha (\\sum_{i=1}^n x_i - k)^2$ \n",
    "\n",
    "to the QUBO, where penalty $P$ is positive whenever the number of $1$s in solution $x_1,...,x_N$ is not $k$, a large enough $\\alpha$ can ensure that such solutions are no longer minima of the problem.  \n",
    "\n",
    "Set set a constraint that $k=2$ with a penalty amplitude of $4$ (you can rerun this cell with varying values of `strength` to see the penalty range from ineffective to overshadowing the problem) and plot the solution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 2\n",
    "bqm.update(dimod.generators.combinations(bqm.variables, k, strength=4))\n",
    "\n",
    "result = sampler.sample(bqm)\n",
    "\n",
    "plot_solutions(result)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example Application: Predicting Survival of Titanic Passengers\n",
    "This example illustrates the MIQUBO method by finding an optimal feature set for predicting survival of Titanic passengers. It uses records provided in file `formatted_titanic.csv`, which is a feature-engineered version of a public database of passenger information recorded by the ship's crew. In addition to a column showing survival for each passenger, its columns record gender, title, class, port of embarkation, etc. \n",
    "\n",
    "The next cell reads the file into a pandas DataFrame. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "titanic = pd.read_csv(\"data/formatted_titanic.csv\") # To see the data folder's contents, select Jupyter File Explorer View from the Online Learning page"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot a ranking of MI between each feature and the variable of interest (survival)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mi = {}\n",
    "features = list(set(titanic.columns).difference(('survived',)))\n",
    "\n",
    "for feature in features:\n",
    "    mi[feature] = mutual_information(prob(titanic[['survived', feature]].values), 1)\n",
    "\n",
    "plot_mi(mi)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exact Versus Good Solutions: a Note on this Dataset\n",
    "This example demonstrates a technique for solving a type of optimization problem on a quantum computer. The Titanic dataset provides a familiar, intuitive example available in the public domain. In itself, however, it is not a good fit for solving by sampling. \n",
    "\n",
    "The next cell plots the values of MI and CMI for all features and three-variable permutations. Notice the clustering of values in tight ranges of values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.plot(range(len(features)), [mutual_information(prob(titanic[['survived', feature]].values), 1) for feature in features], 'bo')\n",
    "\n",
    "plt.plot(range(len([x for x in itertools.combinations(features, 2)])), [conditional_mutual_information(prob(titanic[['survived', f0, f1]].values), 1, 2) for f0, f1 in itertools.combinations(features, 2)], 'go')\n",
    "plt.plot(range(len([x for x in itertools.combinations(features, 2)])), [conditional_mutual_information(prob(titanic[['survived', f1, f0]].values), 1, 2) for f0, f1 in itertools.combinations(features, 2)], 'go')\n",
    "\n",
    "plt.title(\"Titanic MI & CMI Values\")\n",
    "plt.ylabel(\"Shannon Entropy\")\n",
    "plt.xlabel(\"Variable\")\n",
    "plt.legend([\"MI\", \"CMI\"])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot below, obtained by exploiting the problem's small size and brute-force solving for all possible values, shows the solution space for a couple of choices of $k$. The left side shows the resulting energy for all possible assignments of values to $x_1...x_N$ (yellow) and those that satisfy the requirement of $n \\choose k$ (blue); the right side focuses on only those that satisfy $n \\choose k$ and highlights the optimal solution (red).\n",
    "\n",
    "<img src=\"images/k4_7_solution_space.png\" width=800x>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the high number of valid solutions that form a small cluster (the energy difference between the five best solutions in the depicted graph is in the fourth decimal place). The quantum computer's strength is in quickly finding diverse good solutions to hard problems, it is not best employed as a double-precision numerical calculator. Run naively on this dataset, it finds numerous good solutions but is unlikely to find the exact optimal solution.\n",
    "\n",
    "There are many techniques for reformulating problems for the D-Wave system that can improve performance on various metrics, some of which can help narrow down good solutions to closer approach an optimal solution. These are out of scope for this example. For more information, see Leap's other Jupyter Notebooks, the [D-Wave Problem-Solving Handbook](https://docs.dwavesys.com/docs/latest/doc_handbook.html), and examples in the [Ocean software documentation](https://docs.ocean.dwavesys.com/en/stable/index.html).\n",
    "\n",
    "The remainder of this section solves the problem for just the highest-scoring features."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building the MI-Based BQM\n",
    "Select 8 features with the top MI ranking found above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "keep = 8\n",
    "\n",
    "sorted_mi = sorted(mi.items(), key=lambda pair: pair[1], reverse=True)\n",
    "titanic = titanic[[column[0] for column in sorted_mi[0:keep]] + [\"survived\"]]\n",
    "features = list(set(titanic.columns).difference(('survived',)))\n",
    "\n",
    "print(\"Submitting for {} features: {}\".format(keep, features))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate a BQM based on the problem's MI and CMI as done previously for the toy problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers.draw import plot_bqm \n",
    "\n",
    "bqm = dimod.BinaryQuadraticModel.empty(dimod.BINARY)\n",
    "\n",
    "# add the features\n",
    "for feature in features:\n",
    "    mi = mutual_information(prob(titanic[['survived', feature]].values), 1)\n",
    "    bqm.add_variable(feature, -mi)\n",
    "\n",
    "for f0, f1 in itertools.combinations(features, 2):\n",
    "    cmi_01 = conditional_mutual_information(prob(titanic[['survived', f0, f1]].values), 1, 2)\n",
    "    cmi_10 = conditional_mutual_information(prob(titanic[['survived', f1, f0]].values), 1, 2)\n",
    "    bqm.add_interaction(f0, f1, -cmi_01)\n",
    "    bqm.add_interaction(f1, f0, -cmi_10)\n",
    "\n",
    "bqm.normalize()  \n",
    "\n",
    "plot_bqm(bqm)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting Up a QPU as a Solver\n",
    "Set up a D-Wave system as your solver in the standard way described in the Ocean documentation's [Configuring Access to D-Wave Solvers](https://docs.ocean.dwavesys.com/en/stable/overview/sapi.html). \n",
    "\n",
    "*minor-embedding*, the mapping between the problem's variables to the D-Wave QPU's numerically indexed qubits, can be handled in a variety of ways and this affects solution quality and performance. Ocean software provides tools suited for different types of problems; for example, [dwave-system](https://docs.ocean.dwavesys.com/en/stable/docs_system/sdk_index.html) [EmbeddingComposite()](https://docs.ocean.dwavesys.com/en/stable/docs_system/reference/composites.html) has a heuristic for automatic embedding. This example uses [FixedEmbeddingComposite()](https://docs.ocean.dwavesys.com/en/stable/docs_system/reference/composites.html) with the embedding found using an algorithm tuned for cliques (complete graphs).\n",
    "\n",
    "The code below creates a [*dwave-networkx*](https://docs.ocean.dwavesys.com/en/stable/docs_dnx/sdk_index.html) graph that represents the *working graph* of the QPU selected by [DWaveSampler](https://docs.ocean.dwavesys.com/en/stable/docs_system/reference/samplers.html), such as a Pegasus graph with the same sets of nodes (qubits) and edges (couplers) as the QPU. Ocean software's [*minorminer*](https://docs.ocean.dwavesys.com/en/stable/docs_minorminer/sdk_index.html) finds an embedding for the required clique size in the working graph. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dwave.system import DWaveSampler, FixedEmbeddingComposite\n",
    "from minorminer.busclique import find_clique_embedding\n",
    "\n",
    "qpu = DWaveSampler()\n",
    "\n",
    "qpu_working_graph = qpu.to_networkx_graph()\n",
    "embedding = find_clique_embedding(bqm.variables, qpu_working_graph)\n",
    "\n",
    "qpu_sampler = FixedEmbeddingComposite(qpu, embedding)\n",
    "\n",
    "print(\"Maximum chain length for minor embedding is {}.\".format(max(len(x) for x in embedding.values())))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This problem is small enough to be solved in its entirety on an Advantage QPU. For datasets with higher numbers of features, D-Wave Ocean's [dwave-hybrid](https://docs.ocean.dwavesys.com/en/stable/docs_hybrid/sdk_index.html) tool can be used to break the BQM into smaller pieces for serial submission to the QPU and/or parallel solution on classical resources. Here, an out-of-the-box hybrid sampler, [Kerberos](https://docs.ocean.dwavesys.com/en/stable/docs_hybrid/reference/reference.html) is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hybrid.reference.kerberos import KerberosSampler\n",
    "\n",
    "kerberos_sampler = KerberosSampler() "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Submit the Problem for All k Values\n",
    "For all numbers of selected features, $k$, set a $n \\choose k$ penalty, submit an updated BQM for solution, and at the end plot the selected features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers.draw import plot_feature_selection \n",
    "\n",
    "selected_features = np.zeros((len(features), len(features)))\n",
    "for k in range(1, len(features) + 1):\n",
    "    print(\"Submitting for k={}\".format(k))\n",
    "    kbqm = dimod.generators.combinations(features, k, strength=6)\n",
    "    kbqm.update(bqm)\n",
    "    kbqm.normalize()\n",
    "    \n",
    "    best = kerberos_sampler.sample(kbqm, \n",
    "                                   qpu_sampler=qpu_sampler, \n",
    "                                   qpu_reads=5000, \n",
    "                                   max_iter=1,\n",
    "                                   qpu_params={'label': 'Notebook - Feature Selection'}\n",
    "                                  ).first.sample\n",
    "    \n",
    "    for fi, f in enumerate(features):\n",
    "        selected_features[k-1, fi] = best[f]\n",
    "\n",
    "plot_feature_selection(features, selected_features)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright &copy; 2021 D-Wave Systems, Inc\n",
    "\n",
    "The software is licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "you may not use this file except in compliance with the License.\n",
    "You may obtain a copy of the License at\n",
    "\n",
    "    http://www.apache.org/licenses/LICENSE-2.0\n",
    "\n",
    "Unless required by applicable law or agreed to in writing, software\n",
    "distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "See the License for the specific language governing permissions and\n",
    "limitations under the License.\n",
    "\n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc-sa/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by-nc-sa/4.0/88x31.png\" /></a><br />This Jupyter Notebook is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc-sa/4.0/\">Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License</a>"
   ]
  }
 ],
 "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.5"
  },
  "toc": {
   "base_numbering": "1",
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}