nb-fs-topn.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter-Based Feature Selection, select top N-features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we will use a public dataset to perform a selection of features using the top N-features that are best informative for the class vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The website https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/ is a collection of datasets for classification and regression. We will use some of them to test our feature selection algorithms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import urllib\n",
    "\n",
    "filename = \"german.numer_scale\"  # 1000 x 24\n",
    "url = \"https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/\" + filename\n",
    "f = urllib.urlretrieve(url, filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing of the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_MLlib_ relies on _LabeledPoint_ as data structure to that stores a numerical vector (dense or sparse) and a numerical label. An RDD of LabeledPoint represents the dataset given as input to train or test supervised machine learning models.\n",
    "\n",
    "Spark provides a built-in function to tranforms a libsvm dataset into a RDD[LabeledPoint]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pyspark.mllib.util import MLUtils\n",
    "from pyspark.mllib.regression import LabeledPoint\n",
    "\n",
    "rdd = MLUtils.loadLibSVMFile(sc, filename)\n",
    "ncols = rdd.first().features.size  # number of columns (no class) of the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LabeledPoint(-1.0, (24,[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23],[-1.0,-0.941176,1.0,-0.89011,1.0,1.0,0.333333,1.0,-1.0,0.714286,1.0,-0.333333,-1.0,1.0,-1.0,-1.0,-1.0,1.0,-1.0,-1.0,1.0,-1.0,-1.0,1.0]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rdd.first()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pearson correlation coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we create first the Pearson correlation coefficients (PCCs) between the class and each features (_scoreClass_), then the PCCs between every pair of feature (_scoreMatrix_). Once these intermediate results are completed, we proceed into performing the feature selection."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here below we define the functions used during the _map_ and _reduce_ stage of the distributed calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.stats.stats import pearsonr\n",
    "from sklearn.metrics import normalized_mutual_info_score\n",
    "from sklearn.feature_selection import mutual_info_regression\n",
    "import numpy as np\n",
    "\n",
    "def meltLPclass(lp):\n",
    "    '''\n",
    "    This function creates a list of k,v tuples, one per each\n",
    "    label-feature combination. 'k' corresponds to the index\n",
    "    of the feature and 'v' corresponds to a tuple of two\n",
    "    elements: value of the label, value of the feature\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lp : LabeledPoint\n",
    "        a point in the feature space with label\n",
    "    '''\n",
    "    label = lp.label\n",
    "    features = lp.features\n",
    "    r = range(features.size)\n",
    "    return [(i, (label, features[i])) for i in r]\n",
    "\n",
    "def corr(x):\n",
    "    '''\n",
    "    This function calculates the Pearson correlation coefficient\n",
    "    among two variables. It returns the index of the feature and\n",
    "    its correlation coefficient\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x : tuple\n",
    "        x[0] is a scalar value (or a tuple), representing the index(es) of the feature(s)\n",
    "        x[1] is a pyspark.resultiterable.ResultIterable object\n",
    "    '''\n",
    "    idx = x[0]\n",
    "    values = list(x[1])\n",
    "    \n",
    "    l = list(values)\n",
    "    v1, v2 = zip(*values)\n",
    "    p = pearsonr(v1, v2)[0]\n",
    "    \n",
    "    return (idx, p)\n",
    "\n",
    "def mi(x):\n",
    "    '''\n",
    "    This function calculates the Mutual information between two discrete variables.\n",
    "    It relies on the 'adjusted_mutual_info_score' function available in the\n",
    "    sklearn package.\n",
    "    \n",
    "    Doc: http://scikit-learn.org/stable/modules/generated/sklearn.metrics.mutual_info_score.html\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x : tuple\n",
    "        x[0] is a scalar value (or a tuple), representing the index(es) of the feature(s)\n",
    "        x[1] is a pyspark.resultiterable.ResultIterable object\n",
    "    '''\n",
    "    idx = x[0]\n",
    "    values = list(x[1])\n",
    "    \n",
    "    l = list(values)\n",
    "    v1, v2 = zip(*values)\n",
    "    res = normalized_mutual_info_score(v1, v2)\n",
    "    \n",
    "    return (idx, res)\n",
    "\n",
    "def miCont(x):\n",
    "    '''\n",
    "    This function calculates the Mutual information between two continuous variables.\n",
    "    It relies on the 'mutual_info_regression' function available in the\n",
    "    sklearn package.\n",
    "    \n",
    "    Doc: http://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.mutual_info_regression.html\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x : tuple\n",
    "        x[0] is a scalar value (or a tuple), representing the index(es) of the feature(s)\n",
    "        x[1] is a pyspark.resultiterable.ResultIterable object\n",
    "    '''\n",
    "    idx = x[0]\n",
    "    values = list(x[1])\n",
    "    \n",
    "    l = list(values)\n",
    "    v1, v2 = zip(*values)\n",
    "    V1 = np.array(v1).reshape(len(v1),1)  # 1-column matrix layout transformation\n",
    "    res = mutual_info_regression(V1, v2, discrete_features=False, random_state=42)[0]  # random_state is set to provide deterministic results\n",
    "    \n",
    "    return (idx, res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From our RDD (_rdd_) that represents our dataset,\n",
    "- the _flatMap_ operation iterates over every single instance and produce intermediate tuples containing the values of the pair label-feature and the feature index. We need the feature index to be the key of the tuple, so we can group every tuple regarding such feature in the Reducers.\n",
    "- the _groupByKey_ operation sort and gather tuples having the same key in the Reducers (one Reducer per each key)\n",
    "- the _map_ operation of feature _j_ has been provided with all the data label-feature of the feature _j_. That is the vectors of the label and the feature _j_. Having such data in one place, the _corr_ function can calculate the Pearson correlation coefficient. You can plug in other (custom) functions to assess the association of the two variables, such as _mi_ and _miCont_.\n",
    "- the _collect_ operation with collect all the resulting data in the spark driver process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fscores = rdd.flatMap(meltLPclass).groupByKey().map(corr).collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, -0.35084747767184321),\n",
       " (2, -0.22878473305454464),\n",
       " (4, -0.1789427359379214),\n",
       " (6, -0.088184281454553884),\n",
       " (8, 0.14261199150183776)]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fscores[0:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of the distributed computation, the order of the scores in _fscores_ can be different with respect to the order of the features. This is the reason for which _corr_ function returns the feature index along with the Pearson correlation coefficient. We therefore need to sort the data according to the feature index. The result is the _scoreClass_ vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fscoresIdx, fscoresScore = zip(*fscores)\n",
    "scoreClass = [fscoresScore[fscoresIdx.index(i)] for i in range(ncols)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 2, 4)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fscoresIdx[0:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.35084747767184321, -0.22878473305454464, -0.1789427359379214)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fscoresScore[0:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-0.35084747767184321, 0.21492668774990711, -0.22878473305454464]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scoreClass[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Feature Selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Top n-features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the scores label-feature, we select the top _nfs_ features that best correlate with the label. _fsIdx_ stores the indexes of the selected features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nfs = 5  # number of feature to select"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df = zip([abs(x) for x in scoreClass], range(len(scoreClass)))\n",
    "df.sort(key=lambda tup: tup[0], reverse=True)\n",
    "fsIdx = [x[1] for x in df[0:nfs]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0.35084747767184321, 0),\n",
       " (0.22878473305454464, 2),\n",
       " (0.21492668774990711, 1),\n",
       " (0.1789427359379214, 4),\n",
       " (0.15406676409013534, 3)]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 2, 1, 4, 3]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fsIdx[0:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final step is to reduce the dimensionality of _rdd_ according to the selected features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pyspark.mllib.linalg import Vectors\n",
    "\n",
    "def reduceLP(lp, fsIdx):\n",
    "    label = lp.label\n",
    "    features = lp.features\n",
    "    v = [features[i] for i in fsIdx]\n",
    "    return LabeledPoint(label, Vectors.dense(v))\n",
    "\n",
    "rddFS = rdd.map(lambda x: reduceLP(x, fsIdx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LabeledPoint(-1.0, [-1.0,1.0,-0.941176,1.0,-0.89011])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rddFS.first()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Food for thought"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. does this algorithm scale with the number of instances? That is, given 100M instances instead of the current 1K, will this code work or will it crash?\n",
    "2. what about the same question above concerning the number of feature instead?\n",
    "3. for this feature selection, do we need to rely on the LabeledPoint data structure? Can we use another (simpler) data structure?\n",
    "4. can I directly calculate the correlation in _map_ phase instead of going through the _map_ and _reduce_ phases? Why?\n",
    "5. this dataset has 1000 instances, 24 features and 1 class feature. Can you calcuate the number of tuples produced at the _flatMap_ operation? Can you estimate the ratio between the size of such intermediate results and the original dataset?"
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}