notebooks/chapter07_book.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Second order training of (deep) generalized linear model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to GLM\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model for count observations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example of linear poisson regression with artifical data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gc\n",
    "gc.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import psutil\n",
    "memory = psutil.virtual_memory()\n",
    "print(f\" Memory used      : {memory.percent} %\\n\",\n",
    "      f\"Memory available : { round(memory.free / (1024.0 ** 3),2)} GB\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!nvidia-smi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The datasets are loaded as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n  = 5000 #number of rows\n",
    "p1 = 7 #with intercept\n",
    "p  = p1-1 #number of vars\n",
    "\n",
    "print(n, p1, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train   = np.loadtxt(towdir(\"poisson_n5000_d7_Xtrain.txt\"))\n",
    "X_test    = np.loadtxt(towdir(\"poisson_n5000_d7_Xtest.txt\"))\n",
    "y_train   = np.loadtxt(towdir(\"poisson_n5000_d7_ytrain.txt\")).astype(np.int64)\n",
    "y_test    = np.loadtxt(towdir(\"poisson_n5000_d7_ytest.txt\")).astype(np.int64)\n",
    "beta0     = np.loadtxt(towdir(\"poisson_n5000_d7_beta0.txt\"))\n",
    "mu0       = np.loadtxt(towdir(\"poisson_n5000_d7_mu0.txt\"))\n",
    "idx_test  = np.loadtxt(towdir(\"poisson_n5000_d7_idxtest.txt\")).astype(np.int64)\n",
    "idx_train = np.loadtxt(towdir(\"poisson_n5000_d7_idxtrain.txt\")).astype(np.int64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_train, p_train = X_train.shape\n",
    "n_test, p_test   = X_test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson regression fitting with numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Checking the target variable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Training of the regression model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as stm\n",
    "ols= stm.Poisson(y_train, X_train)\n",
    "fit_ols_train = ols.fit()\n",
    "olssumy= fit_ols_train.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_stm = fit_ols_train.params\n",
    "beta_stm = beta_stm.reshape((len(beta_stm),1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Poisson fitting with numpy and full batches or mini-batches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import gammaln\n",
    "\n",
    "def f_poisson_logLik(beta,X,y,name=None):\n",
    "    beta = beta.reshape(len(beta),1)\n",
    "    y = y.reshape((len(y),1)).astype(np.float64)\n",
    "    mu_hat  = np.exp(X @ beta) #.ravel()\n",
    "    logL    = np.sum(y * np.log(mu_hat) - mu_hat - gammaln(y+1))\n",
    "    if name is not None: print(name+\"=\",np.round(logL,4))\n",
    "    return logL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_mu_mse_cor_poisson(X,y,fit,mu0,isprint=None):\n",
    "    beta = fit[\"beta\"]\n",
    "    algo = fit[\"algo\"]\n",
    "    mu_hat  = np.exp(X @ beta).ravel()\n",
    "    mse_mu_hat = ( (mu_hat-mu0.ravel())**2 ).mean()\n",
    "    cor_mu_hat = np.corrcoef(mu0.ravel(),mu_hat)[0,1]\n",
    "    logLik = f_poisson_logLik(beta,X,y)\n",
    "    \n",
    "    if isprint is not None:\n",
    "        print(str(\"mse_mu_\"+algo+\"=\"),np.round(mse_mu_hat,4), \n",
    "              #str(\"cor(cor_mu_\"+namethod+\",mu)=\"),np.round(cor_mu_hat,4),\n",
    "              str(\"logLik_\"+algo+\"=\"),np.round(logLik,4))\n",
    "    \n",
    "    return {\"mu\":mu_hat, \"msemu\":mse_mu_hat, \n",
    "            \"cormu\":cor_mu_hat, \"logL\":logLik,\n",
    "            \"fit\":fit}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Poisson training with pytorch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First order procedure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import deepglmlib.utils as utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch.nn as nn\n",
    "import torch\n",
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = X_train[0:500,:]\n",
    "y_train = y_train[0:500]\n",
    "X_test  = X_test[0:500,:]\n",
    "y_test  = y_test[0:500]\n",
    "\n",
    "idx_train = idx_train[0:500]\n",
    "idx_test = idx_test[0:500]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.utils.data import DataLoader, TensorDataset\n",
    "\n",
    "dt_train = TensorDataset( torch.from_numpy(X_train[:,1:].astype(np.float32)), \n",
    "                          torch.from_numpy(y_train.astype(np.float32)) )\n",
    "\n",
    "dt_test  = TensorDataset( torch.from_numpy(X_test[:,1:].astype(np.float32)), \n",
    "                          torch.from_numpy(y_test.astype(np.float32)) )\n",
    "\n",
    "batch_size= 8\n",
    "dl_train = DataLoader(dt_train, batch_size= batch_size, shuffle=False,num_workers=1)\n",
    "dl_test  = DataLoader(dt_test, batch_size= batch_size, shuffle=False,num_workers=1)\n",
    "n_train, p_train = dl_train.dataset.tensors[0].shape\n",
    "n_test, p_test   = dl_test.dataset.tensors[0].shape\n",
    "\n",
    "print(n_train, p_train, n_test, p_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_update_model(model,loss,optimizer,device,b=None,Xb=None,yb=None):\n",
    "    alpha_t = next(iter(optimizer.param_groups))['lr'] #here constant!\n",
    "    for p in iter(model.parameters()):\n",
    "        p.grad[p.grad>2] = 2\n",
    "        p.grad[p.grad<-2] = -2\n",
    "        p.data = p.data - alpha_t * p.grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
    "\n",
    "name_model = \"PoissonRegression\"\n",
    "nbmax_epoqs = 80\n",
    "debug_out   = 1\n",
    "alpha_t     = 0.0001\n",
    "\n",
    "layers = []\n",
    "layers.append(nn.Linear(p_train,1, bias=True))\n",
    "\n",
    "resus_gdth = utils.f_train_my_glm(dl_train, dl_test, layers, name_model,                   \n",
    "               nbmax_epoqs=nbmax_epoqs, debug_out=debug_out, \n",
    "               alpha_t=alpha_t, transform_yb = utils.transform_yb,\n",
    "               transform_yhatb = utils.transform_yhatb, device=device, \n",
    "               update_model=f_update_model,printed=2,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun_model2vector(model):\n",
    "    beta = [p.detach().numpy().ravel() \n",
    "               for p in model.parameters()]\n",
    "    beta = [beta[(i + 1) % len(beta)]\n",
    "               for i, x in enumerate(beta)]\n",
    "    beta = np.concatenate( beta, axis=0 )\n",
    "    return beta\n",
    "\n",
    "beta_gdth = fun_model2vector(resus_gdth[\"model\"].to(torch.device(\"cpu\"))).ravel()\n",
    "beta_gdth = beta_gdth.reshape(len(beta_gdth),1)\n",
    "\n",
    "fit_gdth         = {\"beta\":beta_gdth,\"algo\":\"gdth\"}\n",
    "quali_gdth_test  = f_mu_mse_cor_poisson(X_test,y_test,fit_gdth,mu0[idx_test])\n",
    "quali_gdth_train = f_mu_mse_cor_poisson(X_train,y_train,fit_gdth,mu0[idx_train])\n",
    "\n",
    "logL_gdth_train = quali_gdth_train[\"logL\"]\n",
    "logL_gdth_test = quali_gdth_test[\"logL\"]\n",
    "print(f\"logL_gdth_train= {logL_gdth_train:5.2f}\")\n",
    "print(f\"logL_gdth_test = {logL_gdth_test:5.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as stm\n",
    "ols= stm.Poisson(y_train, X_train)\n",
    "fit_ols_train = ols.fit()\n",
    "beta_stm = fit_ols_train.params\n",
    "beta_stm = beta_stm.reshape((len(beta_stm),1))\n",
    "\n",
    "fit_stm        = {\"beta\":beta_stm,\"algo\":\"stm\"}\n",
    "quali_stm_test = f_mu_mse_cor_poisson(X_test,y_test,fit_stm,mu0[idx_test])\n",
    "quali_stm_train = f_mu_mse_cor_poisson(X_train,y_train,fit_stm,mu0[idx_train])\n",
    "\n",
    "logL_stm_train = quali_stm_train[\"logL\"]\n",
    "logL_stm_test = quali_stm_test[\"logL\"]\n",
    "print(f\"logL_stm_train = {logL_stm_train:5.2f}\")\n",
    "print(f\"logL_stm_test  = {logL_stm_test:5.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example of training with a pytorch optimizer at second order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda_l1 = 0.01\n",
    "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",
    "    #lossb_b_l1 = (lossb_b_l1 + torch.abs(list(model.parameters())[1])+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": [
    "dl_train = DataLoader(dt_train, batch_size= 32, shuffle=False,num_workers=1)\n",
    "dl_test  = DataLoader(dt_test, batch_size= 32, shuffle=False,num_workers=1)\n",
    "\n",
    "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
    "\n",
    "name_model = \"PoissonRegression\"\n",
    "nbmax_epoqs = 8\n",
    "debug_out   = 1\n",
    "alpha_t     = 0.001\n",
    "\n",
    "layers_regress = []\n",
    "layers_regress.append(nn.Linear(p_train,1, bias=True))\n",
    "\n",
    "resus_lbfgsth = utils.f_train_my_glm(dl_train, dl_test, layers_regress, name_model,                   \n",
    "               nbmax_epoqs=nbmax_epoqs, debug_out=debug_out, \n",
    "               alpha_t=alpha_t, transform_yb = utils.transform_yb,\n",
    "               transform_yhatb = utils.transform_yhatb, device=device, \n",
    "               update_model=f_update_model,printed=1,\n",
    "               name_optimizer=\"LBFGS\", nbmax_iter_lbgs=30,\n",
    "               loss_yy_model=loss_yy_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_lbfgsth = fun_model2vector(resus_lbfgsth[\"model\"].to(torch.device(\"cpu\"))).ravel()\n",
    "beta_lbfgsth = beta_lbfgsth.reshape(len(beta_lbfgsth),1)\n",
    "\n",
    "fit_lbfgsth         = {\"beta\":beta_lbfgsth,\"algo\":\"lbfgsth\"}\n",
    "quali_lbfgsth_test  = f_mu_mse_cor_poisson(X_test,y_test,fit_lbfgsth,mu0[idx_test])\n",
    "quali_lbfgsth_train = f_mu_mse_cor_poisson(X_train,y_train,fit_lbfgsth,mu0[idx_train])\n",
    "\n",
    "logL_lbfgsth_train = quali_lbfgsth_train[\"logL\"]\n",
    "logL_lbfgsth_test = quali_lbfgsth_test[\"logL\"]\n",
    "print(f\"logL_lbfgsth_train = {logL_lbfgsth_train:5.2f}\")\n",
    "print(f\"logL_lbfgsth_test  = {logL_lbfgsth_test:5.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " \n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "method_s = [\"gd-mb-torch\", \n",
    "            \"lbfgs-mb-torch\", \n",
    "            \"stm (module)\"]\n",
    "\n",
    "logLik_s = [quali_gdth_test[\"logL\"],\n",
    "            quali_lbfgsth_test[\"logL\"],\n",
    "            quali_stm_test[\"logL\"]]\n",
    "\n",
    "mse_mu_s = [quali_gdth_test[\"msemu\"], \n",
    "            quali_lbfgsth_test[\"msemu\"],\n",
    "            quali_stm_test[\"msemu\"]]\n",
    "\n",
    "nbstep_s = [resus_gdth[\"tmax\"], \n",
    "            resus_lbfgsth[\"tmax\"], 4]\n",
    "\n",
    "n_train_s = [n_train,n_train,n_train]\n",
    "n_test_s  = [n_test,n_test,n_test]\n",
    "p_s       = [p_train,p_train,p_train]\n",
    "\n",
    "results = [method_s, logLik_s, mse_mu_s, nbstep_s,\n",
    "           n_train_s, n_test_s, p_s]\n",
    "\n",
    "results_pd = pd.DataFrame(results).transpose()\n",
    "results_pd.columns = [\"algo\", \"logL_te\", \n",
    "                      \"mse(mu_hat,mu)_te\", \"nb_steps_tr\",\n",
    "                      \"n_train\", \"n_test\", \"nb_vars\"]\n",
    "\n",
    "with pd.option_context('float_format', '{:.4f}'.format, \n",
    "                       'display.expand_frame_repr', False):\n",
    "    print(results_pd.to_string(index=False))#, header=False"
   ]
  },
  {
   "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
}