notebooks/chapter05_book.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lasso selection for linear regression/classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def towdir(s):\n",
    "    return (str('./datasets_book/'+s))\n",
    "\n",
    "import deepglmlib.utils as utils\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib\n",
    "importlib.reload(utils)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "for filenamefix in {'','_train','_test'}:\n",
    "    with open(\"\".join([towdir('x_y'),filenamefix,'_450d_lasso.npz']), 'rb') as f:\n",
    "        xy = np.load(f)\n",
    "        namex_ = \"\".join(['x',filenamefix]); x_ = xy[namex_]\n",
    "        namey_ = \"\".join(['y',filenamefix]); y_ = xy[namey_]\n",
    "        print(namex_, x_.shape, namey_, y_.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple regression without lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import deepglmlib.utils as utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.utils.data import DataLoader, TensorDataset\n",
    "import torch\n",
    "\n",
    "x_y_train_450d_lasso=np.load(towdir('./x_y_train_450d_lasso.npz'))\n",
    "x_y_test_450d_lasso=np.load(towdir('./x_y_test_450d_lasso.npz'))\n",
    "\n",
    "x_train = x_y_train_450d_lasso['x_train']\n",
    "y_train = x_y_train_450d_lasso['y_train']\n",
    "\n",
    "x_test = x_y_test_450d_lasso['x_test']\n",
    "y_test = x_y_test_450d_lasso['y_test']\n",
    "\n",
    "\n",
    "dataset_train = TensorDataset( torch.from_numpy(x_train.astype(np.float32)), \n",
    "                               torch.from_numpy(y_train.astype(np.float32)) )\n",
    "\n",
    "dataset_test  = TensorDataset( torch.from_numpy(x_test.astype(np.float32)), \n",
    "                               torch.from_numpy(y_test.astype(np.float32)) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_train, p_train = x_train.shape\n",
    "n_test, p_test   = x_test.shape\n",
    "\n",
    "n_train, p_train, n_test, p_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "# import glmlib.utils as utils\n",
    "from torch.utils.data import TensorDataset, DataLoader\n",
    "dl_train            = DataLoader(dataset_train,shuffle=False,batch_size=10)#,num_workers=1,pin_memory=True)\n",
    "dl_test             = DataLoader(dataset_test,shuffle=False,batch_size=10)#,num_workers=1,pin_memory=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"cuda.is_available()     = \", torch.cuda.is_available())\n",
    "# print(torch.cuda.device(0)\n",
    "# print(torch.cuda.device_count())\n",
    "print(\"cuda.get_device_name(0) = \",torch.cuda.get_device_name(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The training function is now implementing with the gpu device for faster training, when available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = torch.device(\"cuda:0\" \\\n",
    "  if torch.cuda.is_available() else \"cpu\")\n",
    "\n",
    "device"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch.nn as nn\n",
    "import copy\n",
    "\n",
    "#px = dataset_train.x.shape[1]\n",
    "px = p_train\n",
    "\n",
    "layers_regress = []\n",
    "layers_regress.append(nn.Linear(px,1,bias=True))\n",
    "\n",
    "model =  utils.GNLMRegression(\"LinearRegression\",\n",
    "                        copy.deepcopy(layers_regress))\n",
    "\n",
    "model.to(device)\n",
    "print(model)\n",
    "\n",
    "\n",
    "nbmax_epoqs = 500\n",
    "alpha_t     = 1e-4\n",
    "debug_out   = 5\n",
    "model.train()\n",
    "\n",
    "loss      = torch.nn.MSELoss(reduction='sum')\n",
    "optimizer = torch.optim.SGD(model.parameters(), lr=alpha_t, momentum=0.0)\n",
    "monitor   = utils.MyMonitorTest(model,loss,dl_train,dl_test,nbmax_epoqs,\n",
    "                                debug_out,device=device)\n",
    "\n",
    "# loss_s,tmax,monistopc  = utils.f_train_glmr(dl_train,model,optimizer,loss,monitor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_train_s,tmax,monistopc  = utils.f_train_glmr(dl_train,model,optimizer,\n",
    "                                       loss,monitor,device=device,printed=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lasso selection for neural networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda_l1 = 0.020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lambda_l1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_yy_model(lossb,model):\n",
    "    lossb_b_rg = lossb\n",
    "    lossb_b_l1 = (torch.abs(list(model.parameters())[0])+0.000001).sum()\n",
    "    loss_b = lossb_b_rg + lambda_l1 * lossb_b_l1\n",
    "    return loss_b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbmax_epoqs = 350\n",
    "alpha_t     = 0.001\n",
    "\n",
    "model          = utils.GNLMRegression(\"LinearRegression\",copy.deepcopy(layers_regress))\n",
    "loss           = torch.nn.MSELoss(reduction='mean')\n",
    "optimizer      = torch.optim.SGD(model.parameters(), lr=alpha_t, momentum=0.0)\n",
    "monitor        = utils.MyMonitorTest(model,loss,dl_train,dl_test,nbmax_epoqs,debug_out,device)\n",
    "\n",
    "loss_train_s,tmax,monistopc = \\\n",
    "    utils.f_train_glmr(dl_train,model,optimizer,loss,monitor,device=device,\n",
    "                       loss_yy_model=loss_yy_model,printed=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Post-processing and mean square error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat_train_l1, y_train_l1 = utils.f_get_yhat(model.cpu(),dl_train)\n",
    "yhat_test_l1, y_test_l1   = utils.f_get_yhat(model.cpu(),dl_test)\n",
    "\n",
    "y_train_l1     = y_train_l1.squeeze()\n",
    "y_test_l1      = y_test_l1.squeeze()\n",
    "yhat_train_l1  = yhat_train_l1.squeeze()\n",
    "yhat_test_l1   = yhat_test_l1.squeeze()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_,_, = utils.f_metrics_regression(y_test_l1,yhat_test_l1,True)\n",
    "_,_, = utils.f_metrics_regression(y_train_l1,yhat_train_l1,True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretation of the obtained model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- <table><tr>\n",
    "<td> <img src=\"./images/pytorch_nn_450d_lasso_without_l1.png\" alt=\"Drawing\" style=\"width: 350px;\"/> </td>\n",
    "<td> <img src=\"./images/pytorch_nn_450d_lasso_with_l1.png\" alt=\"Drawing\" style=\"width: 350px;\"/> </td>\n",
    "</tr></table> -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search for optimal learning rate via a grid from parameters ranges with hdf5 files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.utils.data import DataLoader, TensorDataset\n",
    "import torch\n",
    "\n",
    "x_y_all_450d_lasso=np.load(towdir('./x_y_450d_lasso.npz'))\n",
    "\n",
    "x = x_y_all_450d_lasso['x']\n",
    "y = x_y_all_450d_lasso['y']\n",
    "\n",
    "dataset = TensorDataset( torch.from_numpy(x.astype(np.float32)), \n",
    "                         torch.from_numpy(y.astype(np.float32)) )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_all, d_all = x.shape\n",
    "print( n_all, d_all )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_fold = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.random as rd\n",
    "\n",
    "def f_idx_traintest_kfolds(n,k_fold=5,shuffle = True):\n",
    "    if not shuffle : idx_all = range(n)\n",
    "    if shuffle     : idx_all = rd.permutation(range(n))\n",
    "    idx_s = dict()\n",
    "    for k,idx_test in enumerate(np.array_split(idx_all,k_fold)):\n",
    "        idx_train = [e for e in idx_all if e not in idx_test]\n",
    "        idx_s[str(k)] = dict({\"train\":np.asarray(idx_train),\n",
    "                              \"test\":np.asarray(idx_test)})\n",
    "    return idx_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_s = f_idx_traintest_kfolds(n_all,k_fold=k_fold)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_t     = 0.001\n",
    "lambda_l1   = 0.02\n",
    "nbmax_epoqs = 350\n",
    "batch_size  = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The call to the function is thus: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch.nn as nn\n",
    "import copy\n",
    "layers_regress = []\n",
    "layers_regress.append(nn.Linear(px,1,bias=True))\n",
    "\n",
    "model_          = utils.GNLMRegression(\"LinearRegression\",copy.deepcopy(layers_regress))\n",
    "loss_           = torch.nn.MSELoss(reduction='sum')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_train_s_s, loss_test_s_s, yhat_train_s, \\\n",
    "y_train_s, yhat_test_s, y_test_s = \\\n",
    "    utils.f_reg_l1_nn_cv(idx_s,None,dataset,model_,loss_,batch_size,alpha_t,\n",
    "                   nbmax_epoqs,debug_out,device=device,\n",
    "                  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_train_s_s, loss_test_s_s, yhat_train_s, \\\n",
    "y_train_s, yhat_test_s, y_test_s = \\\n",
    "    utils.f_reg_l1_nn_cv(idx_s,None,dataset,model_,loss_,batch_size,alpha_t,\n",
    "                   nbmax_epoqs,debug_out,device=device,\n",
    "                  loss_yy_model=loss_yy_model,printed=2,\n",
    "                  hyperparameter_to_print=str(f\"lambda_l1={lambda_l1}\")\n",
    "                  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_t_s = np.array([0.1,0.075,0.05,0.01,0.0075,0.005,0.003,0.001,0.00075,0.0005,0.0001])[::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "namefile_s = utils.f_save_cvdatasets_to_h5py(idx_s,dataset,\n",
    "                                       towdir(\"x_y__450d_lasso__\"),\n",
    "                                       x.shape[1],\n",
    "                                       transformx=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "namefile_s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " \n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_t_s   = np.array([0.00005,0.000075,0.001,0.00125,0.00150])\n",
    "lambda_l1_s = np.array([0.001,0.005,0.01,0.015,0.02,0.025,0.03,0.04,0.05,0.1,0.2,0.3]) #[::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "nbmax_epoqs = 100\n",
    "\n",
    "resu_s = []\n",
    "para_s = []\n",
    "\n",
    "for lambda_l1 in iter(lambda_l1_s):\n",
    "    for alpha_t in iter(alpha_t_s):\n",
    "\n",
    "        print(f\"lambda_l1={lambda_l1} alpha_t={alpha_t}\", end='')\n",
    "\n",
    "        loss_train_s_s, loss_test_s_s, yhat_train_s, \\\n",
    "        y_train_s, yhat_test_s, y_test_s = \\\n",
    "        utils.f_reg_l1_nn_cv(idx_s,None,dataset,model_,loss_,batch_size,alpha_t,\n",
    "                       nbmax_epoqs,debug_out,device=device,\n",
    "                       loss_yy_model=loss_yy_model,printed=0,\n",
    "                       hyperparameter_to_print=str(f\"lambda_l1={lambda_l1}\")\n",
    "                       )\n",
    "        \n",
    "        mse_train_s = []\n",
    "        for y_train_, yhat_train_ in iter(zip(y_train_s,yhat_train_s)):\n",
    "            mse_,r2_ = utils.f_metrics_regression(y_train_,yhat_train_,False)\n",
    "            mse_train_s.append(mse_)\n",
    "        \n",
    "        mse_test_s = []\n",
    "        for y_test_, yhat_test_ in iter(zip(y_test_s,yhat_test_s)):\n",
    "            mse_,r2_ = utils.f_metrics_regression(y_test_,yhat_test_,False)\n",
    "            mse_test_s.append(mse_)\n",
    "            \n",
    "        resu_s.append([mse_train_s,mse_test_s,yhat_train_s,y_train_s,yhat_test_s,y_test_s])\n",
    "        para_s.append([lambda_l1,alpha_t])\n",
    "        \n",
    "        msetr = round(np.mean(mse_train_s),4)\n",
    "        msete = round(np.mean(mse_test_s),4)\n",
    "        \n",
    "        print(\" mse_tr_mean=\", msetr, end='')\n",
    "        print(\" mse_te_mean=\", msete, end='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "meanmsetrain = np.asarray([np.mean(resu_[0]) for resu_ in resu_s])\n",
    "stdmsetrain  = np.asarray([np.std(resu_[0]) for resu_ in resu_s])\n",
    "meanmsetest  = np.asarray([np.mean(resu_[1]) for resu_ in resu_s])\n",
    "stdmsetest   = np.asarray([np.std(resu_[1]) for resu_ in resu_s])\n",
    "\n",
    "results = [ [t[0] for t in para_s],[t[1] for t in para_s],\n",
    "           meanmsetrain, meanmsetest, stdmsetrain, stdmsetest, \n",
    "           stdmsetrain/meanmsetrain, stdmsetest/meanmsetest]\n",
    "results_pd = pd.DataFrame(results).transpose()\n",
    "results_pd.columns = [\"lambda_l1\",\"alpha_t\",\n",
    "                      \"mse_tr_mean\", \"mse_te_mean\", \n",
    "                      \"mse_tr_std\", \"mse_te_std\", \"mse_tr_coeffvar\", \n",
    "                      \"mse_te_coeffvar\",]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "g = sns.lineplot(data=results_pd, x='lambda_l1', y='mse_te_mean', style='alpha_t')\n",
    "\n",
    "# g.set(xscale='log')\n",
    "g.set(xticks=results_pd['lambda_l1'])\n",
    "g.set(xticklabels=results_pd['lambda_l1'])\n",
    "plt.xticks(rotation=90)\n",
    "plt.legend(bbox_to_anchor=(1.01, 1), borderaxespad=0, title = \"alpha_t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "myenv",
   "language": "python",
   "name": "myenv"
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}