Select_K_best-updated version.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import multiscale_graphcorr\n",
    "from scipy._lib._util import MapWrapper\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {},
   "outputs": [],
   "source": [
    "def k_sample_test(X, y,score_func=\"mgc\"):\n",
    "    \"\"\"Compute the multivariate independence test statisitic for the provided sample.\n",
    "    Read more in the :ref:`User Guide <multivariate_feature_selection>`.\n",
    "    Parameters\n",
    "    ----------\n",
    "    X : {array-like, sparse matrix} of shape (n_samples, n_features)\n",
    "        The set of regressors that will be tested sequentially.\n",
    "    y : ndarray of shape (n_samples,)\n",
    "        The target vector.\n",
    "    score_func : string that refers to a k_sample multivariate independence test from scipy\n",
    "                The default and only existing test is multiscale graph correlation.\n",
    "    Returns\n",
    "    -------\n",
    "    stat : float that refers to the computed k_sample test statistic\n",
    "    \"\"\"\n",
    "    k_array = np.unique(y)\n",
    "    matrices = []\n",
    "    for i in k_array:\n",
    "        indices = np.where(y == i)[0] \n",
    "        if len(X.shape) == 1:\n",
    "            xi = X[indices]\n",
    "        else:\n",
    "            xi = X[indices,:]\n",
    "        matrices.append(xi)\n",
    "    X = np.concatenate(matrices)\n",
    "    vs = []\n",
    "    for i in range(len(np.unique(y))):\n",
    "        n = matrices[i].shape[0]\n",
    "        encode = np.zeros(shape=(n, len(matrices)))\n",
    "        encode[:, i] = np.ones(shape=n)\n",
    "        vs.append(encode)\n",
    "    y = np.concatenate(vs)\n",
    "    warnings.filterwarnings(\"ignore\")\n",
    "    mgc = multiscale_graphcorr(X,y,reps = 0)\n",
    "    stat = mgc.stat #mgc case\n",
    "    return(stat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_transform(X,y,k,workers = -1): \n",
    "    features = np.arange(X.shape[1])#features of matrix X(n samples by p features, creates an array of feature indexes\n",
    "    best_features = []\n",
    "    while (len(best_features) < k): # outerloop to obtain up to k features \n",
    "        X_new = np.array(X)# establishes matrix of data as np array\n",
    "        parallel = _Parallel(X_new=X_new, y=y,best_features = best_features)#establishes parallel operation helper class object\n",
    "        with MapWrapper(workers) as mapwrapper:\n",
    "            scores = list(mapwrapper(parallel, features)) #maps in parallel the parallel operation that calcs score with the iterable features list to test, with best_features already obtained\n",
    "        scores_index = np.zeros((len(features),2)) #temp array\n",
    "        scores_index[:,0] = features #input features as first column, all features tested(exclude best_features)\n",
    "        scores_index[:,1] = scores #input scores in second column\n",
    "        sorted_index = scores_index[scores_index[:, 1].argsort()] #sort by scores column\n",
    "        best = sorted_index[len(scores)-1,0] #find best of the scores\n",
    "        best_features.append(int(best)) #append new best feature column index \n",
    "        features = np.delete(features,np.where(features == best))\n",
    "    return X_new[:,best_features] # obtain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 215,
   "metadata": {},
   "outputs": [],
   "source": [
    "class _Parallel:\n",
    "    \"\"\"Helper function to calculate parallel test value.\"\"\"\n",
    "\n",
    "    def __init__(self, X_new, y,best_features):\n",
    "        self.X_new = X_new\n",
    "        self.y = y\n",
    "        self.best_features = best_features\n",
    "\n",
    "    def __call__(self, index):\n",
    "        if np.var(self.X_new[:,index]) == 0:\n",
    "            stat = -1000.0\n",
    "        else:   \n",
    "            if len(self.best_features)==0:\n",
    "                X_j =  self.X_new[:,index] #each feature from j to last feature \n",
    "                stat= k_sample_test(X_j,self.y)\n",
    "            else:\n",
    "                columns = self.best_features #construct array for indexing \n",
    "                columns.append(index)\n",
    "                X_j = self.X_new[:,columns]#perform test with obtained features against every feature after to then obtain the best group of features with one additional feature \n",
    "                stat= k_sample_test(X_j,self.y)\n",
    "        return stat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  \n",
      "/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:11: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[13.5745, 12.5243, 14.6717],\n",
       "       [13.2476, 11.8767, 14.715 ],\n",
       "       [12.8502, 11.3969, 14.1845],\n",
       "       [13.0014, 11.621 , 13.6435],\n",
       "       [12.9124, 11.8306, 15.0428],\n",
       "       [13.0333, 11.9655, 14.4929],\n",
       "       [12.2114, 11.3213, 14.0294],\n",
       "       [12.6004, 10.997 , 14.8682],\n",
       "       [12.4893, 11.6024, 14.455 ],\n",
       "       [12.3688, 11.4932, 14.6224],\n",
       "       [12.8225, 11.8076, 14.5764],\n",
       "       [11.7899, 11.4094, 14.3924],\n",
       "       [12.2337, 11.2851, 14.2587],\n",
       "       [11.7723, 11.5346, 14.0707],\n",
       "       [13.0272, 11.7531, 14.5536],\n",
       "       [12.8708, 11.7705, 14.2633],\n",
       "       [12.4727, 11.1859, 14.0748],\n",
       "       [12.1196, 12.1531, 13.7019],\n",
       "       [13.0948, 11.9245, 14.9422],\n",
       "       [12.2289, 10.9051, 14.3723],\n",
       "       [12.4988, 11.6073, 14.7345],\n",
       "       [11.858 , 10.7416, 14.3773],\n",
       "       [13.047 , 11.3527, 14.1223],\n",
       "       [12.8655, 11.9285, 13.8704],\n",
       "       [13.041 , 11.35  , 15.072 ],\n",
       "       [12.9453, 11.5694, 14.1182],\n",
       "       [13.0688, 10.9557, 13.8358],\n",
       "       [12.9483, 11.5178, 14.4335],\n",
       "       [13.2116, 11.1974, 14.0294],\n",
       "       [13.3686, 11.8125, 14.8615],\n",
       "       [13.7047, 11.5599, 15.4929],\n",
       "       [12.6136, 11.337 , 13.9539],\n",
       "       [12.511 , 11.5683, 13.6258],\n",
       "       [12.3943, 11.4127, 14.0707],\n",
       "       [12.3429, 10.5155, 13.9094],\n",
       "       [12.702 , 11.3413, 14.5039],\n",
       "       [12.5139, 11.7966, 14.5093],\n",
       "       [13.072 , 11.8469, 13.9228],\n",
       "       [13.3953, 11.6531, 14.8477],\n",
       "       [12.2041, 11.3449, 14.7025],\n",
       "       [13.858 , 11.9151, 14.6407],\n",
       "       [12.5656, 11.6717, 14.787 ],\n",
       "       [13.098 , 11.4604, 14.6467],\n",
       "       [13.4444, 11.7214, 14.4982],\n",
       "       [12.4435, 11.2566, 14.2813],\n",
       "       [13.2697, 11.5224, 14.4234],\n",
       "       [12.3753, 11.0585, 14.3481],\n",
       "       [12.4206, 13.1244, 14.6224],\n",
       "       [12.4338, 11.1026, 14.4075],\n",
       "       [12.3763, 11.5493, 14.4982],\n",
       "       [12.1211, 11.0454, 14.477 ],\n",
       "       [12.6059, 11.2178, 14.9422],\n",
       "       [12.8816, 11.3778, 14.4335],\n",
       "       [14.1679, 10.7825, 15.8757],\n",
       "       [12.783 , 10.8412, 14.3675],\n",
       "       [12.3534, 11.3025, 14.4819],\n",
       "       [12.7806, 11.686 , 14.9921],\n",
       "       [12.9165, 11.5922, 14.3097],\n",
       "       [11.9438, 10.9903, 14.2813],\n",
       "       [12.4435, 11.0932, 14.888 ],\n",
       "       [12.7223, 11.2515, 13.9063],\n",
       "       [13.1112, 11.0348, 15.1972],\n",
       "       [12.2975, 11.3744, 14.6656],\n",
       "       [12.4271, 10.7479, 14.3723],\n",
       "       [12.5182, 11.0524, 14.3005],\n",
       "       [12.4038, 10.869 , 14.8066],\n",
       "       [12.5253, 10.8467, 14.418 ],\n",
       "       [11.4372, 11.0688, 14.3577],\n",
       "       [14.1508, 10.7526, 15.2341],\n",
       "       [12.5698, 11.2583, 15.1306],\n",
       "       [11.9438, 11.5123, 14.7671],\n",
       "       [13.3001, 11.8903, 14.793 ],\n",
       "       [12.0547, 10.9848, 14.6349],\n",
       "       [12.4387, 11.2375, 14.5536],\n",
       "       [12.3135, 10.9398, 14.6407],\n",
       "       [12.2727, 10.7812, 14.5258],\n",
       "       [12.1813, 11.5577, 14.787 ],\n",
       "       [13.0766, 11.3178, 14.6349],\n",
       "       [12.3366, 11.8522, 13.7074],\n",
       "       [12.3242, 11.263 , 14.2549],\n",
       "       [13.3041, 11.8651, 14.0443],\n",
       "       [12.51  , 11.3308, 14.4283],\n",
       "       [12.5982, 11.0083, 14.7798],\n",
       "       [13.2061, 12.103 , 14.5423],\n",
       "       [13.2623, 11.5411, 14.3097],\n",
       "       [12.0096, 11.5783, 14.3577],\n",
       "       [12.4122, 12.0116, 15.0062],\n",
       "       [12.0715, 11.2465, 14.787 ],\n",
       "       [12.4788, 11.2477, 14.329 ],\n",
       "       [12.4056, 11.4677, 14.4929],\n",
       "       [12.4455, 10.9076, 14.1022],\n",
       "       [12.7931, 11.5736, 13.7669],\n",
       "       [12.4904, 11.2011, 14.3625],\n",
       "       [12.3543, 11.0143, 14.1974],\n",
       "       [11.9204, 11.3993, 14.5093],\n",
       "       [12.3999, 11.6254, 13.9126],\n",
       "       [13.0848, 11.6364, 14.5483],\n",
       "       [12.5872, 11.6717, 13.2257],\n",
       "       [12.4967, 11.0612, 14.0905],\n",
       "       [12.3771, 12.3437, 14.2319],\n",
       "       [12.9732, 11.2187, 14.5996],\n",
       "       [12.853 , 11.96  , 14.793 ],\n",
       "       [12.5505, 11.1192, 14.5423],\n",
       "       [13.1704, 12.3063, 14.5093],\n",
       "       [14.1508, 11.5868, 14.3531],\n",
       "       [12.3269, 11.406 , 13.9857],\n",
       "       [13.1499, 11.745 , 14.0109],\n",
       "       [13.5841, 11.9567, 14.6467],\n",
       "       [13.8515, 12.1673, 14.7025],\n",
       "       [13.0932, 11.7366, 13.9857],\n",
       "       [12.5808, 10.9391, 14.5708],\n",
       "       [12.4893, 11.2746, 14.135 ],\n",
       "       [12.7008, 11.4227, 14.4234],\n",
       "       [12.0665, 11.7081, 13.9892],\n",
       "       [12.1144, 11.8325, 13.4015],\n",
       "       [12.4474, 11.3581, 13.9576],\n",
       "       [12.4999, 11.5376, 13.8048],\n",
       "       [11.9472, 11.7838, 13.8387],\n",
       "       [12.1554, 12.0827, 13.9363],\n",
       "       [12.4884, 12.2913, 14.2232],\n",
       "       [12.0752, 11.1715, 14.0707],\n",
       "       [12.4669, 11.7717, 13.3911],\n",
       "       [11.7777, 11.6014, 13.5768],\n",
       "       [11.7826, 11.1939, 13.6833],\n",
       "       [12.7494, 11.2724, 15.0501],\n",
       "       [12.8475, 11.1156, 14.0294],\n",
       "       [12.2057, 11.7989, 14.3244],\n",
       "       [13.4913, 11.204 , 14.4391],\n",
       "       [12.2794, 10.0261, 14.6967],\n",
       "       [12.8829, 11.2938, 13.9502],\n",
       "       [12.326 , 11.2244, 13.2221],\n",
       "       [11.7285, 10.9667, 14.1223],\n",
       "       [12.07  , 11.3057, 14.3244],\n",
       "       [11.7674,  9.7232, 14.793 ],\n",
       "       [12.4084, 11.5301, 13.8643],\n",
       "       [11.8426, 11.891 , 13.6966],\n",
       "       [11.4893, 11.1885, 14.1928],\n",
       "       [12.4875, 11.296 , 14.3335],\n",
       "       [11.6046, 11.8344, 13.8802],\n",
       "       [12.3743, 11.8281, 14.6467],\n",
       "       [11.8754, 11.191 , 14.0143],\n",
       "       [12.0629, 11.1518, 13.5214],\n",
       "       [12.6953, 11.756 , 13.8802],\n",
       "       [12.5029, 10.9249, 14.7671],\n",
       "       [12.4367, 11.5572, 14.2232],\n",
       "       [12.2557, 11.0184, 14.4234],\n",
       "       [12.5656, 11.2855, 13.877 ],\n",
       "       [11.9996, 11.7692, 14.2861],\n",
       "       [12.7582, 12.2566, 13.8207],\n",
       "       [10.7956, 11.6447, 13.7269],\n",
       "       [12.2984, 11.4898, 13.9927],\n",
       "       [12.3153, 12.1991, 14.1551],\n",
       "       [11.8044, 12.3839, 13.2532],\n",
       "       [12.5232, 11.5704, 14.3531],\n",
       "       [12.5253, 11.7127, 14.418 ],\n",
       "       [12.5893, 11.5689, 14.7537],\n",
       "       [12.3669, 11.59  , 14.0593],\n",
       "       [12.6629, 11.2738, 13.6646],\n",
       "       [11.9743, 10.9851, 13.2459],\n",
       "       [11.8735, 11.7966, 13.613 ],\n",
       "       [12.5982, 11.9533, 13.4937],\n",
       "       [12.8557, 11.9326, 13.4555],\n",
       "       [12.1767, 10.6152, 13.3973],\n",
       "       [12.2213, 11.2877, 13.4555],\n",
       "       [12.1243, 11.1719, 13.9431],\n",
       "       [12.7346, 11.3472, 13.6725],\n",
       "       [11.9696, 12.0504, 13.2642],\n",
       "       [13.3399, 11.1578, 15.3372],\n",
       "       [12.7957, 11.6391, 13.9716],\n",
       "       [12.3437, 11.8767, 13.5431],\n",
       "       [12.0563, 12.0789, 13.4466],\n",
       "       [12.6884, 11.9098, 13.4665],\n",
       "       [12.2098, 12.4126, 13.5987],\n",
       "       [12.6345, 11.3778, 14.8202],\n",
       "       [12.5172, 11.9587, 13.4735],\n",
       "       [12.9568, 11.5646, 14.1803],\n",
       "       [12.5808, 12.7894, 14.477 ],\n",
       "       [13.133 , 12.2821, 13.1046],\n",
       "       [12.242 , 11.4581, 12.8829],\n",
       "       [13.0241, 11.8272, 13.8079],\n",
       "       [13.7241, 12.2804, 14.6833],\n",
       "       [11.3647, 11.6436, 12.5558],\n",
       "       [12.2305, 11.7705, 13.6409],\n",
       "       [12.2993, 11.378 , 13.1261],\n",
       "       [11.8804, 11.122 , 13.8643],\n",
       "       [12.2057, 10.8651, 14.7605],\n",
       "       [11.6062, 11.4591, 14.4929],\n",
       "       [12.6539, 11.3917, 14.4127],\n",
       "       [13.3563, 11.1209, 14.6656],\n",
       "       [12.6238, 11.821 , 13.9609],\n",
       "       [12.1383, 12.0592, 14.1143],\n",
       "       [12.7756, 11.1851, 14.3723],\n",
       "       [12.5452, 11.0318, 14.2364],\n",
       "       [12.4415, 11.3025, 14.6349],\n",
       "       [13.1095, 12.1286, 14.6168],\n",
       "       [13.0932, 11.955 , 13.71  ],\n",
       "       [12.5317, 11.2224, 14.1761],\n",
       "       [13.4015, 12.4038, 13.8515],\n",
       "       [12.7284, 11.4382, 14.3383],\n",
       "       [13.4983, 11.0947, 14.2102],\n",
       "       [12.8867, 11.0505, 13.6081],\n",
       "       [12.6617, 12.3743, 13.9822],\n",
       "       [12.1775, 12.0887, 15.2573],\n",
       "       [12.7743, 11.9703, 12.9568],\n",
       "       [11.8406, 11.5356, 13.5308],\n",
       "       [12.0468, 11.3685, 14.0371],\n",
       "       [12.7743, 11.3518, 13.6966],\n",
       "       [12.4387, 11.0424, 14.6288],\n",
       "       [13.0424, 12.2777, 13.0736],\n",
       "       [13.3563, 12.2913, 13.2679],\n",
       "       [12.8072, 11.0621, 14.5594],\n",
       "       [12.3046, 11.5189, 14.3335]])"
      ]
     },
     "execution_count": 137,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test on dataset from last sem in CMM, high amount of dimensions, low samples select very few features, time about 5 mins\n",
    "import pandas as pd\n",
    "expr_data = pd.read_csv('expr.txt', sep= \"\\t\", header='infer')\n",
    "pheno_data = pd.read_csv('pheno.txt', sep= \"\\t\", header='infer')\n",
    "pheno_data['RelapseGroup_new'] = np.zeros(pheno_data.shape[0])\n",
    "i = 0\n",
    "while (i < pheno_data.shape[0]):\n",
    "    if pheno_data['RelapseGroup'][i] == 'Relapse':\n",
    "        pheno_data['RelapseGroup_new'][i] = 1\n",
    "        i = i + 1\n",
    "    else:\n",
    "        pheno_data['RelapseGroup_new'][i] = 0\n",
    "        i = i + 1\n",
    "expr_data = expr_data.T\n",
    "expr_data.drop('GENE', axis = 0, inplace=True)\n",
    "expr_data = expr_data.astype(float)\n",
    "expr_data['Phenotype'] = np.zeros(pheno_data.shape[0])\n",
    "pheno_data.set_index('GEO_ACCESSION',drop = True, inplace = True)\n",
    "phenotype = pheno_data.loc[:,'RelapseGroup_new']\n",
    "expr_data.index.name = 'GEO_ACCESSION'\n",
    "expr_data['Phenotype'] = phenotype\n",
    "data = np.array(expr_data)\n",
    "X = data[:,0:22215]\n",
    "y = data[:,22215]\n",
    "fit_transform(X,y,3)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5.,  1.,  8., ..., 10.,  0.,  0.],\n",
       "       [ 0., 16.,  0., ..., 16., 16., 16.],\n",
       "       [ 1.,  5.,  0., ..., 11., 13., 15.],\n",
       "       ...,\n",
       "       [ 0., 16.,  5., ..., 13., 16., 16.],\n",
       "       [ 0.,  0.,  7., ..., 16., 14., 14.],\n",
       "       [ 0.,  4.,  0., ..., 14., 15., 15.]])"
      ]
     },
     "execution_count": 216,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test on sklearn dataset, medium-high samples, medium-low features, select like a third of features, 5 mins to run\n",
    "from sklearn.datasets import load_digits\n",
    "X, y = load_digits(return_X_y=True)\n",
    "fit_transform(X,y,20)"
   ]
  },
  {
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}