13_Feature_Engineering_II_Automated_code.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup\n",
    "from mlwpy import *\n",
    "%matplotlib inline\n",
    "\n",
    "# sklearn got chirpy again\n",
    "warnings.filterwarnings('ignore',\n",
    "                        category=sklearn.exceptions.ConvergenceWarning,\n",
    "                        module='sklearn')\n",
    "\n",
    "kwargs = {'test_size':.25, 'random_state':42}\n",
    "\n",
    "iris = datasets.load_iris()\n",
    "tts = skms.train_test_split(iris.data, iris.target, **kwargs)\n",
    "(iris_train,     iris_test, \n",
    " iris_train_tgt, iris_test_tgt) = tts\n",
    "\n",
    "wine = datasets.load_wine()\n",
    "tts = skms.train_test_split(wine.data, wine.target, **kwargs)\n",
    "(wine_train,     wine_test, \n",
    " wine_train_tgt, wine_test_tgt) = tts\n",
    "\n",
    "diabetes = datasets.load_diabetes()\n",
    "tts = skms.train_test_split(diabetes.data, diabetes.target, **kwargs)\n",
    "(diabetes_train_ftrs, diabetes_test_ftrs, \n",
    " diabetes_train_tgt,  diabetes_test_tgt) = tts\n",
    "\n",
    "# these are entire datasets\n",
    "iris_df = pd.DataFrame(iris.data, columns=iris.feature_names)\n",
    "wine_df = pd.DataFrame(wine.data, columns=wine.feature_names)\n",
    "diabetes_df = pd.DataFrame(diabetes.data, columns=diabetes.feature_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(wine_df.var())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(wine_df['hue'].max() - wine_df['hue'].min())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(wine_df['proline'].max() - wine_df['proline'].min())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# variance selection example without scaling\n",
    "varsel = ftr_sel.VarianceThreshold(threshold=1.0)\n",
    "varsel.fit_transform(wine_train)\n",
    "\n",
    "print(\"first example\")\n",
    "print(varsel.fit_transform(wine_train)[0],\n",
    "      wine_train[0, wine_train.var(axis=0) > 1.0], sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(varsel.get_support())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "keepers_idx = varsel.get_support()\n",
    "keepers = np.array(wine.feature_names)[keepers_idx]\n",
    "print(keepers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minmax = skpre.MinMaxScaler().fit_transform(wine_train)\n",
    "print(np.sort(minmax.var(axis=0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# scaled variance selection example\n",
    "pipe = pipeline.make_pipeline(skpre.MinMaxScaler(),\n",
    "                              ftr_sel.VarianceThreshold(threshold=0.05))\n",
    "pipe.fit_transform(wine_train).shape\n",
    "\n",
    "# pipe.steps is list of (name, step_object)\n",
    "keepers_idx = pipe.steps[1][1].get_support()\n",
    "print(np.array(wine.feature_names)[keepers_idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cov(X,Y) = np.dot(X-E(X), Y-E(Y)) / n\n",
    "n = len(diabetes_train_ftrs)\n",
    "\n",
    "# abbreviate names\n",
    "x = diabetes_train_tgt[np.newaxis,:]\n",
    "y = diabetes_train_ftrs\n",
    "cov_via_dot = np.dot(x-x.mean(), y-y.mean()) / n\n",
    "\n",
    "# compute all covariances, extract the ones between ftr and target\n",
    "# bias=True to divide by n instead of n-1; np.cov defaults to bias=False\n",
    "cov_via_np  = np.cov(diabetes_train_ftrs, diabetes_train_tgt, \n",
    "                     rowvar=False, bias=True)[-1, :-1]\n",
    "print(np.allclose(cov_via_dot, cov_via_np))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.var default ddof=0 equates to bias=True\n",
    "# np.corrcoef is a bit of a hot mess to extract values from\n",
    "\n",
    "# cov()/sqrt(var() var())\n",
    "cor_via_cov = cov_via_np / np.sqrt(np.var(diabetes_train_tgt) * \n",
    "                                   np.var(diabetes_train_ftrs, axis=0))\n",
    "cor_via_np = np.corrcoef(diabetes_train_ftrs, diabetes_train_tgt, \n",
    "                         rowvar=False)[-1, :-1]\n",
    "print(np.allclose(cor_via_cov, cor_via_np))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# note, we use the squres of the correlations ... r^2\n",
    "corrs = np.corrcoef(diabetes_train_ftrs, \n",
    "                    diabetes_train_tgt, rowvar=False)[-1, :-1]\n",
    "cor_order = np.argsort(corrs**2) # r^2 (!)\n",
    "cor_names = np.array(diabetes.feature_names)[cor_order[::-1]]\n",
    "\n",
    "# and sklearn's f_regression calculation\n",
    "f_scores = ftr_sel.f_regression(diabetes_train_ftrs, \n",
    "                                diabetes_train_tgt)[0]\n",
    "freg_order = np.argsort(f_scores)\n",
    "freg_names = np.array(diabetes.feature_names)[freg_order[::-1]]\n",
    "\n",
    "# numpy arrays don't like comparing strings\n",
    "print(tuple(cor_names) == tuple(freg_names))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = np.linspace(-10,10,1000).reshape(-1,1)\n",
    "data = np.c_[xs, np.random.uniform(-10,10,xs.shape)]\n",
    "tgt = (np.cos(xs) > 0).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4,3))\n",
    "plt.scatter(data[:,0], tgt);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mi = ftr_sel.mutual_info_classif(data, tgt, \n",
    "                                 discrete_features=False)\n",
    "print(mi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = np.linspace(-10,10,1000).reshape(-1,1)\n",
    "data = np.c_[xs, np.random.uniform(-10,10,xs.shape)]\n",
    "tgt = np.cos(xs).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4,3))\n",
    "plt.plot(data[:,0], tgt);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ftr_sel.f_regression(data, tgt)[0],\n",
    "      ftr_sel.mutual_info_regression(data, tgt), \n",
    "      sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs, ys = np.mgrid[-2:2:.2, -2:2:.2]\n",
    "c_tgt = (ys > xs**2).flatten()\n",
    "\n",
    "# basically, turn r_tgt off if y<x**2\n",
    "r_tgt = ((xs**2 + ys**2)*(ys>xs**2))\n",
    "\n",
    "data = np.c_[xs.flat, ys.flat]\n",
    "\n",
    "# print out a few examples\n",
    "combined = np.c_[data, c_tgt, r_tgt.flat]\n",
    "combined[np.arange(0,401,66)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axes = plt.subplots(1,2, figsize=(6,3), sharey=True)\n",
    "axes[0].scatter(xs, ys, c=np.where(c_tgt, 'r', 'b'), marker='.')\n",
    "axes[0].set_aspect('equal');\n",
    "\n",
    "bound_xs = np.linspace(-np.sqrt(2), np.sqrt(2), 100)\n",
    "bound_ys = bound_xs**2\n",
    "axes[0].plot(bound_xs, bound_ys, 'k')\n",
    "axes[0].set_title('Classification')\n",
    "\n",
    "axes[1].pcolormesh(xs, ys, r_tgt, shading='auto', cmap='binary')\n",
    "axes[1].set_aspect('equal')\n",
    "axes[1].set_title('Regression');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ftr_sel.mutual_info_classif(data, c_tgt, \n",
    "                                  discrete_features=False, random_state=42))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ftr_sel.mutual_info_regression(data, r_tgt.flat, discrete_features=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftrsel = ftr_sel.SelectKBest(ftr_sel.mutual_info_classif, k=5)\n",
    "ftrsel.fit_transform(wine_train, wine_train_tgt)\n",
    "\n",
    "# extract names\n",
    "keepers_idx = ftrsel.get_support()\n",
    "print(np.array(wine.feature_names)[keepers_idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftrsel = ftr_sel.SelectKBest(ftr_sel.f_classif, k=5)\n",
    "ftrsel.fit_transform(wine_train, wine_train_tgt)\n",
    "\n",
    "keepers_idx = ftrsel.get_support()\n",
    "print(np.array(wine.feature_names)[keepers_idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftrsel = ftr_sel.SelectPercentile(ftr_sel.f_regression, \n",
    "                                  percentile=25)\n",
    "ftrsel.fit_transform(diabetes_train_ftrs, \n",
    "                     diabetes_train_tgt)\n",
    "\n",
    "print(np.array(diabetes.feature_names)[ftrsel.get_support()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftrsel = ftr_sel.SelectPercentile(ftr_sel.mutual_info_regression, \n",
    "                                  percentile=25)\n",
    "ftrsel.fit_transform(diabetes_train_ftrs, \n",
    "                     diabetes_train_tgt)\n",
    "\n",
    "print(np.array(diabetes.feature_names)[ftrsel.get_support()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftrsel = ftr_sel.SelectFromModel(ensemble.RandomForestClassifier(), \n",
    "                                 threshold='mean') # default\n",
    "ftrsel.fit_transform(wine_train, wine_train_tgt)\n",
    "print(np.array(wine.feature_names)[ftrsel.get_support()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lmlr = linear_model.LogisticRegression\n",
    "# default solver if lbfgs, which complains about L1 now\n",
    "# so switching to liblinear solver\n",
    "ftrsel = ftr_sel.SelectFromModel(lmlr(penalty='l1',\n",
    "                                 solver='liblinear')) # thesh is \"small\" coeffs\n",
    "ftrsel.fit_transform(wine_train, wine_train_tgt)\n",
    "\n",
    "print(np.array(wine.feature_names)[ftrsel.get_support()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftrsel = ftr_sel.RFE(ensemble.RandomForestClassifier(),\n",
    "                     n_features_to_select=5)\n",
    "\n",
    "res = ftrsel.fit_transform(wine_train, wine_train_tgt)\n",
    "print(np.array(wine.feature_names)[ftrsel.get_support()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# statisticians be warned (see end-of-chapter)\n",
    "# this picks based on feature weights (coefficients)\n",
    "# not on significance of coeffs nor on r^2/anova/F of (whole) model\n",
    "ftrsel = ftr_sel.RFE(linear_model.LinearRegression(),\n",
    "                     n_features_to_select=5)\n",
    "ftrsel.fit_transform(wine_train, wine_train_tgt)\n",
    "print(np.array(wine.feature_names)[ftrsel.get_support()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all the 1s are selected; non-1s are when they were dropped\n",
    "# go to the estimator and ask its coefficients\n",
    "print(ftrsel.ranking_,\n",
    "      ftrsel.estimator_.coef_, sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the order for the five 1s\n",
    "keepers_idx = np.argsort(np.abs(ftrsel.estimator_.coef_))\n",
    "# find the indexes of the 1s and get their ordering\n",
    "keepers_order_idx = np.where(ftrsel.ranking_ == 1)[0][keepers_idx]\n",
    "print(np.array(wine.feature_names)[keepers_order_idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skms.cross_val_score(linear_model.LogisticRegression(), \n",
    "                     wine.data, wine.target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# do it\n",
    "# theshold is \"small\" coeffs\n",
    "lmlr = linear_model.LogisticRegression\n",
    "ftrsel = ftr_sel.SelectFromModel(lmlr(penalty='l1',\n",
    "                                      solver='liblinear'))\n",
    "\n",
    "pipe = pipeline.make_pipeline(ftrsel, linear_model.LogisticRegression())\n",
    "skms.cross_val_score(pipe, wine.data, wine.target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftrsel = ftr_sel.RFE(ensemble.RandomForestClassifier(), \n",
    "                     n_features_to_select=5)\n",
    "pipe = pipeline.make_pipeline(ftrsel, linear_model.LogisticRegression())\n",
    "skms.cross_val_score(pipe, wine.data, wine.target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftrsel = ftr_sel.SelectPercentile(ftr_sel.mutual_info_classif, \n",
    "                                  percentile=25)\n",
    "pipe = pipeline.make_pipeline(ftrsel, linear_model.LogisticRegression())\n",
    "skms.cross_val_score(pipe, wine.data, wine.target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftrsel = ftr_sel.SelectPercentile(ftr_sel.mutual_info_classif, percentile=25)\n",
    "pipe = pipeline.make_pipeline(ftrsel, linear_model.LogisticRegression())\n",
    "\n",
    "param_grid = {'selectpercentile__percentile':[5,10,15,20,25]}\n",
    "grid = skms.GridSearchCV(pipe, param_grid=param_grid, cv=3)\n",
    "grid.fit(wine.data, wine.target)\n",
    "\n",
    "print(grid.best_params_)\n",
    "print(grid.best_score_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs, ys = np.mgrid[-2:2:.2, -2:2:.2]\n",
    "tgt = (xs**2 + ys**2 > 1).flatten()\n",
    "data = np.c_[xs.flat, ys.flat]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(4,3))\n",
    "\n",
    "# draw the points\n",
    "ax.scatter(xs, ys, c=np.where(tgt, 'r', 'b'), marker='.')\n",
    "ax.set_aspect('equal');\n",
    "\n",
    "# draw the circle boundary\n",
    "circ = plt.Circle((0,0), 1, color='k', fill=False)\n",
    "ax.add_patch(circ);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shootout_linear = [svm.SVC(kernel='linear'),\n",
    "                   linear_model.LogisticRegression()]\n",
    "\n",
    "fig, axes = plt.subplots(1,2,figsize=(4,2), sharey=True)\n",
    "for mod, ax in zip(shootout_linear, axes):\n",
    "    plot_boundary(ax, data, tgt, mod, [0,1])\n",
    "    ax.set_title(get_model_name(mod))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create some non-linear learning models\n",
    "knc_p, dtc_p = [1,20], [1,3,5,10]\n",
    "KNC = neighbors.KNeighborsClassifier\n",
    "DTC = tree.DecisionTreeClassifier\n",
    "shootout_nonlin = ([(KNC(n_neighbors=p), p) for p in knc_p] +\n",
    "                   [(DTC(max_depth=p), p)   for p in dtc_p ])\n",
    "\n",
    "# plot 'em\n",
    "fig, axes = plt.subplots(2,3,figsize=(9, 6),\n",
    "                        sharex=True, sharey=True)\n",
    "for (mod, param), ax in zip(shootout_nonlin, axes.flat):\n",
    "    plot_boundary(ax, data, tgt, mod, [0,1])\n",
    "    ax.set_title(get_model_name(mod) + \"({})\".format(param))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_data = np.concatenate([data, data**2], axis=1)\n",
    "print(\"augmented data shape:\", new_data.shape)\n",
    "print(\"first row:\", new_data[0])\n",
    "\n",
    "fig, axes = plt.subplots(1,2,figsize=(5,2.5))\n",
    "for mod, ax in zip(shootout_linear, axes):\n",
    "    # using the squares for prediction and graphing\n",
    "    plot_boundary(ax, new_data, tgt, mod, [2,3]) \n",
    "    ax.set_title(get_model_name(mod))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pick a few points to show off the before/after differences\n",
    "test_points = np.array([[.5,.5],\n",
    "                       [-1, -1.25]])\n",
    "\n",
    "# wonderful trick from trig class:  \n",
    "# if we walk around circle (fractions of pi), \n",
    "# sin/cos give x and y value from pi\n",
    "circle_pts = np.linspace(0,2*np.pi,100)\n",
    "circle_xs, circle_ys = np.sin(circle_pts), np.cos(circle_pts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,2, figsize=(6,3))\n",
    "\n",
    "labels = [('x',     'y',     'Original Space'),\n",
    "          ('$x^2$', '$y^2$', 'Squares Space')]\n",
    "\n",
    "funcs = [lambda x:x,    # for original ftrs\n",
    "         lambda x:x**2] # for squared ftrs\n",
    "\n",
    "for ax, func, lbls in zip(axes, funcs, labels):\n",
    "    ax.plot(func(circle_xs), func(circle_ys), 'k')\n",
    "    ax.scatter(*func(data.T), c=np.where(tgt, 'r', 'b'), marker='.')\n",
    "    ax.scatter(*func(test_points.T), c=['k', 'y'], s=100, marker='^')\n",
    "    \n",
    "    ax.axis('equal')\n",
    "    ax.set_xlabel(lbls[0])\n",
    "    ax.set_ylabel(lbls[1])\n",
    "    ax.set_title(lbls[2])\n",
    "    \n",
    "axes[1].yaxis.tick_right()\n",
    "axes[1].yaxis.set_label_position(\"right\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_data = metrics.pairwise.polynomial_kernel(data, data, \n",
    "                                            degree=2) # squares\n",
    "print('first example: ', data[0])\n",
    "\n",
    "print('example in   original features:', data[0].shape)\n",
    "print('example in kernelized features:', k_data[0].shape)\n",
    "\n",
    "print('# examples in both:', len(data), len(k_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(4,3))\n",
    "\n",
    "# learn from k_data instead of original data\n",
    "preds  = (linear_model.LogisticRegression()\n",
    "                      .fit(k_data, tgt)\n",
    "                      .predict(k_data))\n",
    "\n",
    "ax.scatter(xs, ys, c=np.where(preds, 'r', 'b'), marker='.')\n",
    "ax.set_aspect('equal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.base import TransformerMixin\n",
    "\n",
    "class PolyKernel(TransformerMixin):\n",
    "    def __init__(self, degree):\n",
    "        self.degree = degree\n",
    "\n",
    "    def transform(self, ftrs):\n",
    "        pk = metrics.pairwise.pairwise_kernels\n",
    "        return pk(ftrs, self.ftrs, metric='poly', degree=self.degree)\n",
    "\n",
    "    def fit(self, ftrs, tgt=None):\n",
    "        self.ftrs = ftrs\n",
    "        return self"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import kernel_approximation\n",
    "\n",
    "kn = kernel_approximation.Nystroem(kernel='polynomial', \n",
    "                                   degree=2, n_components=6)\n",
    "LMLR = linear_model.LogisticRegression()\n",
    "k_logreg1 = pipeline.make_pipeline(kn, LMLR)\n",
    "k_logreg2 = pipeline.make_pipeline(PolyKernel(2), LMLR)\n",
    "\n",
    "shootout_fancy = [(k_logreg1, 'Nystroem'),\n",
    "                  (k_logreg2, 'PolyKernel')]\n",
    "\n",
    "fig, axes = plt.subplots(1,2,figsize=(6,3), sharey=True)\n",
    "for (mod, kernel_name), ax in zip(shootout_fancy, axes):\n",
    "    plot_boundary(ax, data, tgt, mod, [0,1])\n",
    "    ax.set_title(get_model_name(mod)+\"({})\".format(kernel_name))\n",
    "    ax.set_aspect('equal')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_logreg = pipeline.make_pipeline(PolyKernel(2),\n",
    "                                  linear_model.LogisticRegression())\n",
    "\n",
    "shootout_fancy = [svm.SVC(kernel='poly', degree=2),\n",
    "                  k_logreg]\n",
    "\n",
    "fig, axes = plt.subplots(1,2,figsize=(6,3), sharey=True)\n",
    "for mod, ax in zip(shootout_fancy, axes):\n",
    "    plot_boundary(ax, data, tgt, mod, [0,1])\n",
    "    ax.set_title(get_model_name(mod))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first three are linear (but different)\n",
    "sv_classifiers = {\"LinearSVC\"   : svm.LinearSVC(),\n",
    "                  \"SVC(linear)\" : svm.SVC(kernel='linear'),\n",
    "                  \n",
    "                  \"SVC(poly=1)\" : svm.SVC(kernel='poly', degree=1), \n",
    "                  \"SVC(poly=2)\" : svm.SVC(kernel='poly', degree=2),\n",
    "                  \"SVC(poly=3)\" : svm.SVC(kernel='poly', degree=3),\n",
    "                  \n",
    "                  \"SVC(rbf,.5)\" : svm.SVC(kernel='rbf', gamma=0.5),\n",
    "                  \"SVC(rbf,1.0)\": svm.SVC(kernel='rbf', gamma=1.0),\n",
    "                  \"SVC(rbf,2.0)\": svm.SVC(kernel='rbf', gamma=2.0)}\n",
    "\n",
    "fig, axes = plt.subplots(4,2,figsize=(8,8),sharex=True, sharey=True)\n",
    "for ax, (name, mod) in zip(axes.flat, sv_classifiers.items()):\n",
    "    plot_boundary(ax, iris.data, iris.target, mod, [0,1])\n",
    "    ax.set_title(name)\n",
    "plt.tight_layout() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "digits = datasets.load_digits()\n",
    "\n",
    "param_grid = {\"gamma\" : np.logspace(-10, 1, 11, base=2),\n",
    "              \"C\"     : [0.5, 1.0, 2.0]}\n",
    "\n",
    "svc = svm.SVC(kernel='rbf')\n",
    "\n",
    "grid_model = skms.GridSearchCV(svc, param_grid = param_grid, \n",
    "                               cv=10)\n",
    "grid_model.fit(digits.data, digits.target);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid_model.best_params_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_gamma = grid_model.best_params_['gamma']\n",
    "my_svc = svm.SVC(kernel='rbf', **grid_model.best_params_)\n",
    "scores = skms.cross_val_score(my_svc,\n",
    "                              digits.data, digits.target, \n",
    "                              cv=10,\n",
    "                              scoring='f1_macro')\n",
    "print(scores)\n",
    "print(\"{:5.3f}\".format(scores.mean()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.array([[1, 2, 4, 5], \n",
    "                   [2.5,.75,5.25,3.5]]).T\n",
    "mean = data.mean(axis=0)\n",
    "centered_data = data - mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(3,3))\n",
    "\n",
    "# original data in red; mean is larger dot\n",
    "ax.plot(*data.T, 'r.')\n",
    "ax.plot(*mean, 'ro')\n",
    "\n",
    "# centered data in blue, mean is larger dot (and at (0,0))\n",
    "ax.plot(*centered_data.T, 'b.')\n",
    "ax.plot(*centered_data.mean(axis=0), 'bo')\n",
    "\n",
    "#ax.set_aspect('equal');\n",
    "high_school_style(ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can reproduce the original data\n",
    "fig,ax = plt.subplots(figsize=(3,3))\n",
    "orig_data = centered_data + mean\n",
    "plt.plot(*orig_data.T, 'r.')\n",
    "plt.plot(*orig_data.mean(axis=0), 'ro')\n",
    "\n",
    "ax.set_xlim((0,6))\n",
    "ax.set_ylim((0,6))\n",
    "high_school_style(ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,3,figsize=(9,3), sharex=True,sharey=True)\n",
    "xs = np.linspace(0,6,30)\n",
    "lines = [(1,0), (0,3), (-1,6)]\n",
    "data = np.array([[1, 2, 4, 5], \n",
    "                   [2.5,.75,5.25,3.5]]).T\n",
    "plot_lines_and_projections(axes, lines, data, xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# draw data\n",
    "ax = plt.gca()\n",
    "ax.scatter(data[:,0], data[:,1], c='b', marker='.')\n",
    "\n",
    "# draw mean\n",
    "mean = np.mean(data, axis=0, keepdims=True)\n",
    "centered_data = data - mean\n",
    "ax.scatter(*mean.T, c='k')\n",
    "\n",
    "# compute PCA\n",
    "pca = decomposition.PCA()\n",
    "P = pca.fit_transform(centered_data)\n",
    "\n",
    "# extract useful bits for drawing\n",
    "directions = pca.components_\n",
    "lengths = pca.explained_variance_\n",
    "print(\"Lengths:\", lengths)\n",
    "var_wgt_prindirs = -np.diag(lengths).dot(directions) # negate so point up/right\n",
    "\n",
    "# draw principal axes\n",
    "sane_quiver(var_wgt_prindirs, ax, \n",
    "            origin=np.mean(data, axis=0), \n",
    "            colors='r')\n",
    "ax.set_xlim(0,10)\n",
    "ax.set_ylim(0,10)\n",
    "ax.set_aspect('equal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.allclose(np.cov(data, rowvar=False),\n",
    "                  np.cov(centered_data, rowvar=False)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "orig_cov = np.cov(centered_data, rowvar=False)\n",
    "print(orig_cov)\n",
    "print(np.diag(orig_cov).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EPS = 2.2e-16 # EPS (epsilon) is computer sciency for \"really small\"\n",
    "p_cov = np.cov(P, rowvar=False)\n",
    "p_cov[p_cov<EPS] = 0.0  # throw out \"really small\" values\n",
    "print(p_cov)\n",
    "print(p_cov.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_rotation(theta):\n",
    "    ''' ccw rotation of theta when it post-multiplies \n",
    "        a row vector (an example) '''\n",
    "    return np.array([[np.cos(theta), -np.sin(theta)],\n",
    "                     [np.sin(theta),  np.cos(theta)]]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spacing = np.linspace(0,2*np.pi,17)\n",
    "points  = np.c_[np.sin(spacing), np.cos(spacing)] # sin/cos walk around circle\n",
    "two_points = points[[0,3]]\n",
    "\n",
    "rot = make_rotation(np.pi/8) # 1/16th turn degrees ccw\n",
    "scale = np.diag([2,.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,4,figsize=(8,2), sharey=True)\n",
    "\n",
    "# original vanilla circle\n",
    "axes[0].plot(*two_points.T, 'k^')\n",
    "axes[0].plot(*points.T, 'b.')\n",
    "\n",
    "# make a rotation\n",
    "axes[1].plot(*np.dot(two_points, rot).T, 'k^')\n",
    "axes[1].plot(*np.dot(points, rot).T, 'r.')\n",
    "\n",
    "# stretch along x and y axes\n",
    "axes[2].plot(*two_points.dot(rot).dot(scale).T, 'k^')\n",
    "axes[2].plot(*points.dot(rot).dot(scale).T, 'r.')\n",
    "\n",
    "# undo initial rotation\n",
    "axes[3].plot(*two_points.dot(rot).dot(scale).dot(rot.T).T, 'k^')\n",
    "axes[3].plot(*points.dot(rot).dot(scale).dot(rot.T).T, 'b.')\n",
    "\n",
    "names = ['circle', 'rotate', 'scale', 'unrotate']\n",
    "for ax,name in zip(axes,names):\n",
    "    ax.set_aspect('equal')\n",
    "    ax.set_title(name)\n",
    "    ax.set_xlim(-2,2)\n",
    "    ax.set_ylim(-2,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def signs_like(A, B):\n",
    "    ' produce new A,B with signs of A '\n",
    "    assert np.allclose(np.abs(A), np.abs(B))\n",
    "    signs = np.sign(A) * np.sign(B)\n",
    "    return A, B * signs\n",
    "signs_like([1,-1], [1,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(centered_data)\n",
    "\n",
    "U, s, Vt = np.linalg.svd(centered_data, full_matrices=False)\n",
    "S = np.diag(s)\n",
    "print(np.allclose(centered_data, U.dot(S).dot(Vt)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# align signs - \n",
    "# note, U.S.Vt won't give centered now b/c \n",
    "# U,Vt work together and we monkeyed with Vt\n",
    "_, Vt = signs_like(pca.components_, Vt)   \n",
    "V = Vt.T\n",
    "\n",
    "# directions come from Vt;  amounts come from S\n",
    "# div. by n-1 to get unbiased ... see EOC\n",
    "print(all((np.allclose(pca.components_,         Vt),\n",
    "           np.allclose(pca.explained_variance_, s**2/(N-1)),\n",
    "           np.allclose(P,                       centered_data.dot(V)))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('original covariance:\\n', orig_cov)\n",
    "print('centered covariance:\\n',\n",
    "      centered_data.T.dot(centered_data) / (N-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigval, eigvec = np.linalg.eig(orig_cov)\n",
    "\n",
    "# two differences:\n",
    "# 1.  eigvals aren't ordered high to low (S from svd is ordered)\n",
    "# 2.  final signs with SVD go through V *and* U, eigvec bit different\n",
    "order = np.argsort(eigval)[::-1]\n",
    "print(np.allclose(eigval[order], s**2/(N-1)))\n",
    "\n",
    "_,ev = signs_like(Vt,eigvec[:,order])\n",
    "print(np.allclose(ev, Vt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import sklearn.manifold as manifold\n",
    "data_3d, color = datasets.make_s_curve(n_samples=750, \n",
    "                                       random_state=42)\n",
    "\n",
    "cmap = plt.cm.Spectral\n",
    "fig = plt.figure(figsize=(4,4))\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "ax.scatter(*data_3d.T, c=color, cmap=cmap)\n",
    "ax.view_init(20, -50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,2,figsize=(8,4))\n",
    "n_components = 2\n",
    "\n",
    "# method 1:  map to 2D using isomap\n",
    "isomap = manifold.Isomap(n_neighbors=10, n_components=n_components)\n",
    "data_2d = isomap.fit_transform(data_3d)\n",
    "axes[0].scatter(*data_2d.T, c=color, cmap=cmap)\n",
    "\n",
    "# method 2:  map to 2D using TSNE\n",
    "tsne = manifold.TSNE(n_components=n_components,                     \n",
    "                     init='pca', \n",
    "                     random_state=42)\n",
    "data_2d = tsne.fit_transform(data_3d)\n",
    "axes[1].scatter(*data_2d.T, c=color, cmap=cmap);"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}