Practice 11 - MF with PyTorch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recommender Systems 2022/2023\n",
    "\n",
    "## Practice Session 11 - MF with PyTorch\n",
    "\n",
    "PyTorch, Tensorflow, Keras are useful framework that allow you to build machine learning models (from linear regression to complex deep learning methods) and hide almost all of the complexity related to the training. Usually, you only have to create an object that starting from the model parameters will be able to compute your prediction, then specify the loss and the framework automatically calculates the gradients.\n",
    "\n",
    "#### Performance warning!\n",
    "In image processing tasks one usually has an image, maybe a reasonably large one (1000x1000x3, hence 3\\*10^6 data points) on which a complex network is applied (multiple convolution operations, pooling etc). The computationally expensive part can be effectively parallelized by the framework and hence the speedup over using single-core operations is massive.\n",
    "\n",
    "Unfortunately here we are not dealing with images, but with user profiles. In terms of data this means each profile can be in the 10^5 - 10^6 items. The issue arises when considering the model. If you use a matrix factorization model, the core of the operation is a dot product between two embedding vectors, which is an extremely fast operation. There is hardly any speedup and the burden of the overall infrastructure is not offset by it. For this reason, if you use a profiler you will see that 80-90% of the time is spent in the data sampling phase (because it is done in python it can be quite slow) and the actual prediction computation is a tiny fraction of the time. Overall, a *simple* matrix factorization model may be 10x slower if implemented with pytorch without being careful. You must ensure to use a fast sampling script (Cython here again) and a reasonably powerful GPU.\n",
    "\n",
    "#### Prototyping\n",
    "Given how the complexity of gradients and such is hidden, pytorch becomes a great tool for prototyping. It is very easy to change someting in your model becayse you do not need to dig in Cython code. For example, you may implement a SLIM MSE method that uses as inital parameters the similarity computed with an item-based cosine method, or you may create a hybrid of multiple similarities and lear the weights to use for each similarity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What do we need\n",
    "\n",
    "* A Dataset object to load the data\n",
    "* Model object\n",
    "* Training loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Movielens10M: Verifying data consistency...\n",
      "Movielens10M: Verifying data consistency... Passed!\n",
      "DataReader: current dataset is: Movielens10M\n",
      "\tNumber of items: 10681\n",
      "\tNumber of users: 69878\n",
      "\tNumber of interactions in URM_all: 10000054\n",
      "\tValue range in URM_all: 0.50-5.00\n",
      "\tInteraction density: 1.34E-02\n",
      "\tInteractions per user:\n",
      "\t\t Min: 2.00E+01\n",
      "\t\t Avg: 1.43E+02\n",
      "\t\t Max: 7.36E+03\n",
      "\tInteractions per item:\n",
      "\t\t Min: 0.00E+00\n",
      "\t\t Avg: 9.36E+02\n",
      "\t\t Max: 3.49E+04\n",
      "\tGini Index: 0.57\n",
      "\n",
      "\tICM name: ICM_all, Value range: 1.00 / 69.00, Num features: 10126, feature occurrences: 128384, density 1.19E-03\n",
      "\tICM name: ICM_genres, Value range: 1.00 / 1.00, Num features: 20, feature occurrences: 21564, density 1.01E-01\n",
      "\tICM name: ICM_tags, Value range: 1.00 / 69.00, Num features: 10106, feature occurrences: 106820, density 9.90E-04\n",
      "\tICM name: ICM_year, Value range: 6.00E+00 / 2.01E+03, Num features: 1, feature occurrences: 10681, density 1.00E+00\n",
      "\n",
      "\n",
      "Warning: 82 (0.12 %) of 69878 users have no sampled items\n"
     ]
    }
   ],
   "source": [
    "from Data_manager.split_functions.split_train_validation_random_holdout import split_train_in_two_percentage_global_sample\n",
    "from Data_manager.Movielens.Movielens10MReader import Movielens10MReader\n",
    "\n",
    "data_reader = Movielens10MReader()\n",
    "data_loaded = data_reader.load_data()\n",
    "\n",
    "URM_all = data_loaded.get_URM_all()\n",
    "\n",
    "URM_train, URM_test = split_train_in_two_percentage_global_sample(URM_all, train_percentage = 0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MF models rely upon latent factors for users and items which are called 'embeddings'\n",
    "\n",
    "![latent factors](https://miro.medium.com/max/988/1*tiF4e4Y-wVH732_6TbJVmQ.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_factors = 10\n",
    "\n",
    "n_users, n_items = URM_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "\n",
    "# Creates U\n",
    "user_factors = torch.nn.Embedding(num_embeddings=n_users, embedding_dim=num_factors)\n",
    "\n",
    "# Creates V\n",
    "item_factors = torch.nn.Embedding(num_embeddings=n_items, embedding_dim=num_factors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Embedding(69878, 10)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Embedding(10681, 10)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_factors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In order to compute the prediction you have to:\n",
    "- Get a list of user and item indices (as tensors)\n",
    "- Get the user and item embedding\n",
    "- Compute the element-wise product of the embeddings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(tensor([42]), tensor([42]))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_index = torch.Tensor([42]).type(torch.LongTensor)\n",
    "item_index = torch.Tensor([42]).type(torch.LongTensor)\n",
    "\n",
    "user_index, item_index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notice that each object has a \"grad_fn=...\" attribute, which si going to be used for the automatic gradient compuation to go backwards in the operations required to compute the prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(tensor([[-0.4088,  0.8317,  1.7737, -0.0576,  1.3666, -0.5248, -0.8404, -1.2425,\n",
       "          -0.5362, -0.2207]], grad_fn=<EmbeddingBackward0>),\n",
       " tensor([[-0.3708,  0.5689, -0.5090, -0.4367, -1.1786,  0.8405, -1.5139,  1.2917,\n",
       "          -1.4407, -0.7400]], grad_fn=<EmbeddingBackward0>))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "current_user_factors = user_factors(user_index)\n",
    "current_item_factors = item_factors(item_index)\n",
    "\n",
    "current_user_factors, current_item_factors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the dot product is just a summation over the elementwise product. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(-1.7016, grad_fn=<SumBackward0>)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prediction = torch.mul(current_user_factors, current_item_factors).sum()\n",
    "prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the \"grad_fn\" states \"SubBackward\", the prediction was indeed due to a sum\n",
    "\n",
    "#### We can also use the einstein summation format, which is particularly useful when you have a more complex equation to compute the prediction\n",
    "\n",
    "The einstein summation allows you to write the prediction in terms of the indices of a summation. In this case we want to iterate both embedding vectors, perform an element-by-element product and then sum at the end. Be careful on the dimensions, in this case the factors have two dimensions (the row dimension is 1 so in practice it is useless). We use \"b\" to iterate over the rows (useful when we compute batches of predictions to parallelize) and \"i\" is the latent factor index.\n",
    "\n",
    "$$ r_b = \\sum_{i}U'_{bi}V'_{bi} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([-1.7016], grad_fn=<ViewBackward0>)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "torch.einsum(\"bi,bi->b\", current_user_factors, current_item_factors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To take the result of the prediction and transform it into a traditional numpy array you have to:\n",
    "- call .detach() to disconnect the tensor from the automatic gradient tracking\n",
    "- then .numpy()\n",
    "\n",
    "### The result is an array of 1 cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Prediction is -1.70\n"
     ]
    }
   ],
   "source": [
    "prediction_numpy = prediction.detach().numpy()\n",
    "print(\"Prediction is {:.2f}\".format(prediction_numpy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train a MF MSE model with PyTorch\n",
    "\n",
    "# Step 1 Create a Model python object\n",
    "\n",
    "### The model should implement the forward function which computes the prediction as we did before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MF_MSE_PyTorch_model(torch.nn.Module):\n",
    "    def __init__(self, n_users, n_items, n_factors):\n",
    "        super(MF_MSE_PyTorch_model, self).__init__()\n",
    "\n",
    "        self.n_users = n_users\n",
    "        self.n_items = n_items\n",
    "\n",
    "        self.user_factors = torch.nn.Embedding(num_embeddings=self.n_users, embedding_dim=n_factors)\n",
    "        self.item_factors = torch.nn.Embedding(num_embeddings=self.n_items, embedding_dim=n_factors)\n",
    "\n",
    "    def forward(self, user_batch, item_batch):\n",
    "        user_factors_batch = self.user_factors(user_batch)\n",
    "        item_factors_batch = self.item_factors(item_batch)\n",
    "\n",
    "        prediction_batch = torch.mul(user_factors_batch, item_factors_batch).sum()\n",
    "\n",
    "        return prediction_batch\n",
    "\n",
    "    def get_W(self):\n",
    "        return self.user_factors.weight.detach().cpu().numpy()\n",
    "\n",
    "    def get_H(self):\n",
    "        return self.item_factors.weight.detach().cpu().numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2 Setup PyTorch devices and Data Reader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MF_MSE_PyTorch: Using CUDA\n"
     ]
    }
   ],
   "source": [
    "if torch.cuda.is_available():\n",
    "    device = torch.device('cuda')\n",
    "    print(\"MF_MSE_PyTorch: Using CUDA\")\n",
    "else:\n",
    "    device = torch.device('cpu')\n",
    "    print(\"MF_MSE_PyTorch: Using CPU\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create an instance of the model and specify the device it should run on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = MF_MSE_PyTorch_model(n_users, n_items, num_factors).to(device)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Choose loss functions (Mean Squared Error in our case), there are quite a few to choose from"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "lossFunction = torch.nn.MSELoss(reduction=\"sum\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively one can implement it "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _my_MSE_loss(model, user, item):\n",
    "    \n",
    "    # Compute prediction for each element in batch\n",
    "    prediction = model.forward(user, item)\n",
    "\n",
    "    # Compute total loss for batch\n",
    "    loss = (prediction - rating).pow(2).mean()\n",
    "    \n",
    "    return loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select the optimizer to be used for the model parameters: Adam, AdaGrad, RMSProp etc... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "learning_rate = 1e-4\n",
    "l2_reg = 1e-3\n",
    "\n",
    "optimizer = torch.optim.Adagrad(model.parameters(), lr=learning_rate, weight_decay = l2_reg*learning_rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the DatasetInteraction, which will be used to load a specific data point\n",
    "\n",
    "A DatasetInteraction will implement the Dataset class and provide the \\_\\_getitem\\_\\_(self, index) method, which allows to get the data points indexed by that index.\n",
    "\n",
    "Since we need the data to be a tensor, we pre inizialize everything as a tensor. In practice we save the URM in coordinate format (user, item, rating)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.utils.data import Dataset\n",
    "import numpy as np\n",
    "\n",
    "class DatasetInteraction(Dataset):\n",
    "    def __init__(self, URM):\n",
    "\n",
    "        URM = URM.tocoo()\n",
    "        self.n_data_points = URM.nnz\n",
    "\n",
    "        self._row = torch.tensor(URM.row).type(torch.LongTensor)\n",
    "        self._col = torch.tensor(URM.col).type(torch.LongTensor)\n",
    "        self._data = torch.tensor(URM.data).type(torch.FloatTensor)\n",
    "       \n",
    "    def __getitem__(self, index):\n",
    "        return self._row[index], self._col[index], self._data[index]\n",
    "\n",
    "\n",
    "    def __len__(self):\n",
    "        return self.n_data_points\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We pass the DatasetIterator to a DataLoader object which manages the use of batches and so on..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.utils.data import DataLoader\n",
    "\n",
    "# A large batch_size (256, 512...) improves parallelization, but the gradient becomes more smoot\n",
    "# at some point the performance will increase but at the expense of the final prediction quality\n",
    "batch_size = 64\n",
    "\n",
    "dataset_iterator = DatasetInteraction(URM_train[0:1000])\n",
    "\n",
    "train_data_loader = DataLoader(dataset=dataset_iterator,\n",
    "                               batch_size=batch_size,\n",
    "                               shuffle=True,\n",
    "                               #num_workers = 2,\n",
    "                              )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## And now we ran the usual epoch steps\n",
    "* Data point sampling\n",
    "* Prediction computation\n",
    "* Loss function computation\n",
    "* Gradient computation\n",
    "* Update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[tensor([417, 173,  67, 442, 655, 821, 507, 486, 842, 507, 279, 591, 820, 243,\n",
       "         284, 405, 600, 855, 891, 139, 201, 233, 465, 233, 746, 604, 729, 111,\n",
       "         841, 551, 326, 717, 456, 422, 377, 451, 300, 929, 820, 379, 573, 405,\n",
       "         940, 605, 422, 407, 731, 280, 764, 828, 350, 758, 825, 402, 507, 377,\n",
       "         446, 823, 781, 562, 406, 770, 781, 862]),\n",
       " tensor([ 417, 1377,  331, 2901,  691, 1055, 3040, 1864,  528,   96, 1994,  691,\n",
       "           82,  213,  367, 1864, 1429, 1288,   16,  799,  466, 3718,   88,   36,\n",
       "         1306, 1040,  811, 1339, 1117,  563,  192, 1376, 1419, 1622, 1008,   31,\n",
       "           39,  302,   72, 1603,  368,  204,  613, 3071, 1139,  300,  113,  411,\n",
       "         1076, 1296,    4,  198, 2015, 1429, 2847, 2192, 1134,  437,   63,   41,\n",
       "          498,  177, 4109,  573]),\n",
       " tensor([4.0000, 4.0000, 4.0000, 3.5000, 4.0000, 2.0000, 4.0000, 3.5000, 1.0000,\n",
       "         5.0000, 4.0000, 2.0000, 2.0000, 4.5000, 4.0000, 2.0000, 5.0000, 4.0000,\n",
       "         4.0000, 4.0000, 5.0000, 2.0000, 3.0000, 4.0000, 4.0000, 5.0000, 4.0000,\n",
       "         3.0000, 1.0000, 4.0000, 5.0000, 3.0000, 3.0000, 3.0000, 4.0000, 5.0000,\n",
       "         3.0000, 4.0000, 4.0000, 4.0000, 2.0000, 2.0000, 4.0000, 4.0000, 3.0000,\n",
       "         4.0000, 4.0000, 4.0000, 3.0000, 2.0000, 2.0000, 3.0000, 4.0000, 5.0000,\n",
       "         4.0000, 4.5000, 5.0000, 3.5000, 4.0000, 3.0000, 3.5000, 5.0000, 3.5000,\n",
       "         4.5000])]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "batch = next(iter(train_data_loader))\n",
    "batch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cda7462426194e3c869b7017342f1bfb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1661 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: total: 4.41 s\n",
      "Wall time: 3.72 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "def run_epoch(data_iterator):\n",
    "    \n",
    "    epoch_loss = 0\n",
    "    for batch in tqdm(data_iterator):\n",
    "\n",
    "        # Clear previously computed gradients\n",
    "        optimizer.zero_grad()\n",
    "\n",
    "        user, item, rating = batch\n",
    "        user = user.to(device)\n",
    "        item = item.to(device)\n",
    "        rating = rating.to(device)\n",
    "\n",
    "        # Compute prediction for each element in batch\n",
    "        prediction = model.forward(user, item)\n",
    "\n",
    "        # Compute total loss for batch\n",
    "        loss = (prediction - rating).pow(2).mean()\n",
    "\n",
    "        # Compute gradients given current loss\n",
    "        loss.backward()\n",
    "\n",
    "        # Apply gradient using the selected optimizer\n",
    "        optimizer.step()\n",
    "\n",
    "        epoch_loss += loss.item()\n",
    "        \n",
    "        \n",
    "run_epoch(train_data_loader)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## After the train is complete (after many epochs), we can get the matrices in the usual numpy format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "user_factors = model.get_W()\n",
    "item_factors = model.get_H()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[ 0.95083886, -1.2302514 , -0.9758749 , ...,  0.43588272,\n",
       "         -1.2674538 , -1.889634  ],\n",
       "        [-0.9309265 ,  0.27153483,  0.50834876, ..., -1.0890548 ,\n",
       "         -0.06128926, -1.1721663 ],\n",
       "        [-0.29489905, -1.9374983 ,  1.8267319 , ..., -1.3253108 ,\n",
       "         -0.69661677, -0.9291013 ],\n",
       "        ...,\n",
       "        [ 0.00301835,  0.18262248, -0.03628764, ...,  0.21849915,\n",
       "          1.1100464 ,  0.07815299],\n",
       "        [-1.0627358 ,  0.2809174 , -0.8795823 , ..., -0.03074418,\n",
       "         -0.23348367, -2.0146549 ],\n",
       "        [ 0.8826899 , -1.3980063 ,  1.2482066 , ...,  1.8058219 ,\n",
       "         -0.33894625,  0.02159044]], dtype=float32),\n",
       " (69878, 10))"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_factors, user_factors.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[ 1.71100342e+00, -1.71951842e+00,  8.17174241e-02, ...,\n",
       "         -1.48534536e+00, -9.88895446e-02, -4.12131786e-01],\n",
       "        [ 1.06533170e+00,  8.93806458e-01,  2.44553253e-01, ...,\n",
       "          5.43021798e-01,  4.13350135e-01, -8.92155409e-01],\n",
       "        [ 9.41068828e-01, -5.04354715e-01,  5.77265263e-01, ...,\n",
       "          1.09992906e-01,  1.11340630e+00, -2.55454421e-01],\n",
       "        ...,\n",
       "        [ 1.23278880e+00,  1.24008727e+00,  1.24930787e+00, ...,\n",
       "         -2.73840427e-01,  7.57642210e-01,  1.29161060e+00],\n",
       "        [-1.34193742e+00,  1.46697676e+00, -1.98895597e+00, ...,\n",
       "          1.57310751e-05,  2.72938877e-01,  2.28593826e+00],\n",
       "        [-9.53395128e-01,  1.13207646e-01, -2.48263434e-01, ...,\n",
       "         -2.05180216e+00, -2.29296136e+00, -7.50487745e-01]], dtype=float32),\n",
       " (10681, 10))"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_factors, item_factors.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's use the profile to monitor which operation is taking the most time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following disables code output line wrapping to better read the profiler lor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "div.output_area pre {\n",
       "    white-space: pre;\n",
       "}\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%html\n",
    "<style>\n",
    "div.output_area pre {\n",
    "    white-space: pre;\n",
    "}\n",
    "</style>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f455a9512da745f4af1aef4987a7534d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1661 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------------------------------------------------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  \n",
      "                                                   Name    Self CPU %      Self CPU   CPU total %     CPU total  CPU time avg     Self CUDA   Self CUDA %    CUDA total  CUDA time avg    # of Calls  \n",
      "-------------------------------------------------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  \n",
      "enumerate(DataLoader)#_SingleProcessDataLoaderIter._...        21.53%        2.470s        65.54%        7.517s       4.523ms     468.773ms         3.78%        6.594s       3.967ms          1662  \n",
      "                                           aten::select        15.14%        1.736s        15.68%        1.798s       5.639us        2.092s        16.88%        2.504s       7.852us        318855  \n",
      "                                            aten::stack        13.42%        1.539s        28.23%        3.238s     649.719us        2.310s        18.64%        3.619s     726.290us          4983  \n",
      "                                        aten::unsqueeze        12.76%        1.464s        13.22%        1.516s       4.756us     849.078ms         6.85%        1.273s       3.992us        318855  \n",
      "                                            aten::copy_        10.79%        1.238s        10.79%        1.238s     124.175us        1.468s        11.85%        1.468s     147.340us          9966  \n",
      "                              aten::_local_scalar_dense         3.49%     400.552ms         3.49%     400.552ms      80.351us        2.298s        18.54%        2.298s     460.890us          4985  \n",
      "                            Optimizer.step#Adagrad.step         3.27%     375.270ms         6.27%     718.766ms     432.731us      41.382ms         0.33%     468.706ms     282.183us          1661  \n",
      "                                              aten::mul         1.69%     194.292ms         1.69%     194.292ms      23.395us        1.034s         8.34%        1.034s     124.454us          8305  \n",
      "                                              aten::cat         1.16%     133.329ms         1.58%     181.736ms      36.471us      10.865ms         0.09%      32.666ms       6.555us          4983  \n",
      "                                       aten::as_strided         1.10%     126.018ms         1.10%     126.018ms       0.194us     855.080ms         6.90%     855.080ms       1.317us        649335  \n",
      "-------------------------------------------------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  \n",
      "Self CPU time total: 11.469s\n",
      "Self CUDA time total: 12.392s\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from torch.profiler import profile, record_function, ProfilerActivity\n",
    "\n",
    "with profile(activities=[ProfilerActivity.CPU, ProfilerActivity.CUDA], \n",
    "             profile_memory=False, record_shapes=False) as prof:\n",
    "\n",
    "    run_epoch(train_data_loader)\n",
    "    \n",
    "print(prof.key_averages().table(sort_by=\"self_cpu_time_total\", row_limit=10))\n",
    "# prof.export_chrome_trace(\"trace.json\")    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how approximately 65\\% of the CPU time is spent by:\n",
    "* SingleProcessDataLoader\n",
    "* aten::select\n",
    "* aten::stack\n",
    "* aten::unsqueeze\n",
    "\n",
    "All of these are either exclusively or mostly related to the data loading and their type conversion to tensors.\n",
    "\n",
    "Only a few percent of the time is spent by the optmizer and actual embedding product. This indicates our model is terribly inefficient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A less terribly inefficient data sampling\n",
    "\n",
    "Despite following the strategy recommended by PyTorch, the previous data sampler is terribly slow likely resulting in 80-90% of the global training time of the model to be wasted waiting for the sampling process. This is partially because the *DataLoader* object calls your custom *Dataset* object to get each single data point and then concatenates them as tensors, making lots of intermediate calls and conversions.\n",
    "\n",
    "You can check this with the simple profiler PyTorch provides.\n",
    "\n",
    "I STRONGLY recommend to use a custom sampler that leverages some basic vectorization to load a full batch of data instead. For example as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.utils.data import Dataset\n",
    "import math\n",
    "import numpy as np\n",
    "import scipy.sparse as sps\n",
    "\n",
    "class InteractionIterator(object):\n",
    "    \"\"\"\n",
    "    This Sampler samples among all the existing user-item interactions *uniformly at random*:\n",
    "    - One of the interactions in the dataset is sampled\n",
    "\n",
    "    The sample is: user_id, item_id, rating\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, URM_train, batch_size = 1):\n",
    "        super(InteractionIterator, self).__init__()\n",
    "\n",
    "        self.URM_train = sps.coo_matrix(URM_train)\n",
    "        self.n_users, self.n_items = self.URM_train.shape\n",
    "        self.n_data_points = self.URM_train.nnz\n",
    "        self.n_sampled_points = 0\n",
    "        self.batch_size = batch_size\n",
    "\n",
    "    def __len__(self):\n",
    "        return math.ceil(self.n_data_points/self.batch_size)\n",
    "\n",
    "    def __iter__(self):\n",
    "        self.n_sampled_points = 0\n",
    "        return self\n",
    "\n",
    "    def __next__(self):\n",
    "\n",
    "        if self.n_sampled_points >= self.n_data_points:\n",
    "            raise StopIteration\n",
    "\n",
    "        this_batch_size = min(self.batch_size, self.n_data_points-self.n_sampled_points)\n",
    "        self.n_sampled_points += this_batch_size\n",
    "\n",
    "        index_batch = np.random.randint(self.n_data_points, size = this_batch_size)\n",
    "\n",
    "        return torch.from_numpy(self.URM_train.row[index_batch]).long(),\\\n",
    "               torch.from_numpy(self.URM_train.col[index_batch]).long(), \\\n",
    "               torch.from_numpy(self.URM_train.data[index_batch]).float()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cb1cd20eff1448569c278ffaafdd8d49",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1661 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------------------------------------------------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  \n",
      "                                                   Name    Self CPU %      Self CPU   CPU total %     CPU total  CPU time avg     Self CUDA   Self CUDA %    CUDA total  CUDA time avg    # of Calls  \n",
      "-------------------------------------------------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  \n",
      "                              aten::_local_scalar_dense        16.08%     541.972ms        16.08%     541.972ms     108.764us     800.528ms        18.51%     800.528ms     160.652us          4983  \n",
      "                                            aten::copy_        14.86%     500.583ms        14.86%     500.583ms      33.486us        1.355s        31.33%        1.355s      90.642us         14949  \n",
      "                            Optimizer.step#Adagrad.step        10.50%     353.954ms        20.85%     702.657ms     423.033us      40.606ms         0.94%     588.233ms     354.144us          1661  \n",
      "                                              aten::mul         5.48%     184.792ms         5.48%     184.792ms      22.251us     976.068ms        22.57%     976.068ms     117.528us          8305  \n",
      "                                         aten::_to_copy         5.04%     169.671ms        20.19%     680.329ms      51.199us      52.101ms         1.20%        1.419s     106.796us         13288  \n",
      "                                             aten::add_         2.93%      98.780ms         5.14%     173.243ms      17.383us     153.672ms         3.55%     182.830ms      18.345us          9966  \n",
      "                                              aten::sum         2.86%      96.220ms         3.00%     101.032ms      30.413us      22.840ms         0.53%      26.984ms       8.123us          3322  \n",
      "                                              aten::pow         2.63%      88.517ms         3.27%     110.113ms      33.147us      19.265ms         0.45%      32.970ms       9.925us          3322  \n",
      "                  Optimizer.zero_grad#Adagrad.zero_grad         2.37%      80.000ms         3.71%     124.943ms      75.222us       8.541ms         0.20%      56.373ms      33.939us          1661  \n",
      "                                               aten::to         2.34%      78.781ms        22.53%     759.110ms      45.702us      38.264ms         0.88%        1.457s      87.740us         16610  \n",
      "-------------------------------------------------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  \n",
      "Self CPU time total: 3.369s\n",
      "Self CUDA time total: 4.325s\n",
      "\n"
     ]
    }
   ],
   "source": [
    "batch_train_data_loader = InteractionIterator(URM_train[0:1000], batch_size=batch_size)\n",
    "\n",
    "with profile(activities=[ProfilerActivity.CPU, ProfilerActivity.CUDA], \n",
    "             profile_memory=False, record_shapes=False) as prof:\n",
    "\n",
    "    run_epoch(batch_train_data_loader)\n",
    "    \n",
    "print(prof.key_averages().table(sort_by=\"self_cpu_time_total\", row_limit=10))\n",
    "# prof.export_chrome_trace(\"trace.json\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much better! We have removed the data bottleneck and the training proceeds remarkably faster!\n",
    "\n",
    "Notice how we still have some data related processing \n",
    "* aten::copy_\n",
    "* aten::_local_scalar_dense \n",
    "* aten::_to_copy\n",
    "\n",
    "These are mostly related to the CPU-GPU data transfer and are still significant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What if I want to change the sampling?\n",
    "The DatasetInteraction can be modified to obtain the desired behaviour, for example adding some negative (zero-rated) items in the sampling. If we want our model to be able to distinguish between positive and negative items we need to let the model see negative data as well, in our case the negative data is the zero-rated items."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "class DatasetInteraction(Dataset):\n",
    "    def __init__(self, URM_train, positive_quota):\n",
    "        \n",
    "        self._URM_train = sps.csr_matrix(URM_train)\n",
    "        \n",
    "        URM_train = URM_train.tocoo()\n",
    "        self.n_data_points = URM.nnz\n",
    "\n",
    "        self._row = torch.tensor(URM_train.row).type(torch.LongTensor)\n",
    "        self._col = torch.tensor(URM_train.col).type(torch.LongTensor)\n",
    "        self._data = torch.tensor(URM_train.data).type(torch.FloatTensor)\n",
    "        self._positive_quota = positive_quota\n",
    "       \n",
    "    def __getitem__(self, index):\n",
    "        select_positive_flag = torch.rand(1, requires_grad=False) > self._positive_quota\n",
    "\n",
    "        if select_positive_flag[0]:\n",
    "            return self._row[index], self._col[index], self._data[index]\n",
    "        else:\n",
    "            user_id = self._row[index]\n",
    "            seen_items = self._URM_train.indices[self._URM_train.indptr[user_id]:self._URM_train.indptr[user_id+1]]\n",
    "            negative_selected = False\n",
    "\n",
    "            while not negative_selected:\n",
    "                negative_candidate = torch.randint(low=0, high=self.n_items, size=(1,))[0]\n",
    "\n",
    "                if negative_candidate not in seen_items:\n",
    "                    item_negative = negative_candidate\n",
    "                    negative_selected = True\n",
    "\n",
    "            return self._row[index], item_negative, torch.tensor(0.0)\n",
    "\n",
    "        \n",
    "        return self._row[index], self._col[index], self._data[index]\n",
    "\n",
    "\n",
    "    def __len__(self):\n",
    "        return self.n_data_points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may also change the dataset iterator to one that samples the user profile rather than the specific interaction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "class UserProfile_Dataset(Dataset):\n",
    "    def __init__(self, URM_train, device):\n",
    "        super().__init__()\n",
    "        URM_train = sps.csr_matrix(URM_train)\n",
    "        self.device = device\n",
    "\n",
    "        self.n_users, self.n_items = URM_train.shape\n",
    "        self._indptr = URM_train.indptr\n",
    "        self._indices = torch.tensor(URM_train.indices, dtype = torch.long, device=device)\n",
    "        self._data = torch.tensor(URM_train.data, dtype = torch.float, device=device)\n",
    "\n",
    "    def __len__(self):\n",
    "        return self.n_users\n",
    "\n",
    "    def __getitem__(self, user_id):\n",
    "        start_pos = self._indptr[user_id]\n",
    "        end_pos = self._indptr[user_id+1]\n",
    "\n",
    "        user_profile = torch.zeros(self.n_items, dtype=torch.float, requires_grad=False, device=self.device)\n",
    "        user_profile[self._indices[start_pos:end_pos]] = self._data[start_pos:end_pos]\n",
    "\n",
    "        return user_profile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What if I want to implement AsySVD? SLIM EN ... \n",
    "You just have to change the pytorch model with the desired one (easy to do)\n",
    "\n",
    "Note these two models work by sampling the whole user profile, they can be adapted to a sampler that provides single interactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch import nn\n",
    "\n",
    "class AsySVDModel(nn.Module):\n",
    "\n",
    "    def __init__(self, embedding_size = None, n_items = None, device = None):\n",
    "        super().__init__()\n",
    "\n",
    "        self._embedding_item_1 = torch.nn.Parameter(torch.randn((n_items, embedding_size)))\n",
    "        self._embedding_item_2 = torch.nn.Parameter(torch.randn((embedding_size, n_items)))\n",
    "\n",
    "    def forward(self, user_profile_batch):\n",
    "        # input shape is batch_size x n items\n",
    "        # r_hat_bi = SUM{e=0}{e=embedding_size} SUM{j=0}{j=n items} r_bj * V1_je * V2_ei\n",
    "        layer_output = torch.einsum(\"bj,je,ei->bi\", user_profile_batch, self._embedding_item_1, self._embedding_item_2)\n",
    "        return layer_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SDenseModel(nn.Module):\n",
    "\n",
    "    def __init__(self, n_items = None, device = None):\n",
    "        super().__init__()\n",
    "\n",
    "        self._S = torch.nn.Parameter(torch.zeros((n_items, n_items)))\n",
    "\n",
    "    def forward(self, user_profile_batch):\n",
    "        # input shape is batch_size x n items\n",
    "        # r_hat_bi = SUM{j=0}{j=n items} r_bj * S_ji\n",
    "        layer_output = torch.einsum(\"bj,ji->bi\", user_profile_batch, self._S)\n",
    "        return layer_output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What if I want to change the loss function?\n",
    "You can just implement the new one, BPR is quite simple. Make sure that the dataset iterator samples the right data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BPR_Dataset(Dataset):\n",
    "    def __init__(self, URM_train):\n",
    "        super().__init__()\n",
    "        self._URM_train = sps.csr_matrix(URM_train)\n",
    "        self.n_users, self.n_items = self._URM_train.shape\n",
    "\n",
    "    def __len__(self):\n",
    "        return self.n_users\n",
    "\n",
    "    def __getitem__(self, user_id):\n",
    "\n",
    "        seen_items = self._URM_train.indices[self._URM_train.indptr[user_id]:self._URM_train.indptr[user_id+1]]\n",
    "        item_positive = np.random.choice(seen_items)\n",
    "\n",
    "        negative_selected = False\n",
    "\n",
    "        while not negative_selected:\n",
    "            negative_candidate = np.random.randint(low=0, high=self.n_items, size=1)[0]\n",
    "\n",
    "            if negative_candidate not in seen_items:\n",
    "                item_negative = negative_candidate\n",
    "                negative_selected = True\n",
    "\n",
    "        return user_id, item_positive, item_negative\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_BPR(model, batch):\n",
    "    user, item_positive, item_negative = batch\n",
    "\n",
    "    # Compute prediction for each element in batch\n",
    "    x_ij = model.forward(user, item_positive) - model.forward(user, item_negative)\n",
    "\n",
    "    # Compute total loss for batch\n",
    "    loss = -x_ij.sigmoid().log().mean()\n",
    "\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}