{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"A co-visitation matrix is essentially an \"analog\" approximation to matrix factorization! I talk a bit more about this idea here: [💡 What is the co-visitation matrix, really?](https://www.kaggle.com/competitions/otto-recommender-system/discussion/365358).\n\nBut matrix factorization has a lot of advantages as compared to co-visitation matrices. First of all, it can make better use of data -- it operates on the notion of similarity between categories. We can construct a more powerful representation if our model understands that aid `1` is similar to aid `142` as opposed to it treating each aid as an atomic entity (this is the jump from unigram/bigram/trigram models to word2vec in NLP).\n\nLet us thus train a matrix factorization model and replace the co-visitation matrices with it!\n\nNow, I don't expect that the first version of the model will be particularly well tuned. There has already been a lot of work put into co-visitation matrices and in the later versions we work off 3 different matrices, one for each category of actions! A similar progression can and will happen with matrix factorization 🙂 This notebook hopefully will enable us to jumpstart this type of exploration 🙂\n\nTo streamline the work, we will use data in `parquet` format. (Here is the notebook [💡 [Howto] Full dataset as parquet/csv files](https://www.kaggle.com/code/radek1/howto-full-dataset-as-parquet-csv-files) and here is [the most up-to-date version of the dataset](https://www.kaggle.com/datasets/radek1/otto-full-optimized-memory-footprint), no need for dealing with `jasonl` files and the associated mess any longer! Please upvote if you find this useful!)\n\nFor data processing we will use [polars](https://www.pola.rs/). `Polars` has a much smaller memory footprint than `pandas` and is quite fast. Plus it has really clean, intuitive API.\n\nLet's get to work! 🙂\n\n**UPDATE:** [@CPMP](https://www.kaggle.com/cpmpml) ported this notebook to run fully on the GPU! 🥳 This is awesome as it can allow you to experiment and ensemble models much faster 🔥 The other notebook also demonstrates how to accelerate a workflow very important to RecSys (and at the heart of the predictions in the notebook you are reading now) -- the nearest neighbor search algorithm. By running it on the GPU not only do you get significantly better results (you don't have to rely on approximate nearest neighbor search anymore, there is enough compute to develiver __true nearest neighbor search faster than you could run ANN on the CPU!__)\n\nPlease find the GPU accelerated notebook here: [💡Matrix Factorization with GPU](https://www.kaggle.com/code/cpmpml/matrix-factorization-with-gpu)\n\nAnd if you would be interested in a further discussion of NN (nearest neighbor) vs ANN (approximate nearest neighbor) and running them the CPU/GPU, please see my post [here](https://www.kaggle.com/competitions/otto-recommender-system/discussion/371111) (though it is a bit outdated as `cuml` has a nicer API).\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"# Data Preprocessing","metadata":{}},{"cell_type":"code","source":"!pip install polars\n\nimport polars as pl\n\ntrain = pl.read_parquet('../input/otto-full-optimized-memory-footprint/train.parquet')\ntest = pl.read_parquet('../input/otto-full-optimized-memory-footprint/test.parquet')","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:30:53.870635Z","iopub.execute_input":"2022-12-11T04:30:53.870974Z","iopub.status.idle":"2022-12-11T04:31:32.556351Z","shell.execute_reply.started":"2022-12-11T04:30:53.870945Z","shell.execute_reply":"2022-12-11T04:31:32.555570Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We need to create `aid-aid` pairs to train our matrix factorization model!\n\nLet's us grab the pairs both from the train and test set.","metadata":{}},{"cell_type":"code","source":"%%time\n\ntrain_pairs = (pl.concat([train, test])\n    .groupby('session').agg([\n        pl.col('aid'),\n        pl.col('aid').shift(-1).alias('aid_next')\n    ])\n    .explode(['aid', 'aid_next'])\n    .drop_nulls()\n)[['aid', 'aid_next']]","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:31:32.558023Z","iopub.execute_input":"2022-12-11T04:31:32.558306Z","iopub.status.idle":"2022-12-11T04:32:04.292930Z","shell.execute_reply.started":"2022-12-11T04:31:32.558275Z","shell.execute_reply":"2022-12-11T04:32:04.290830Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_pairs.shape[0] / 1_000_000","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:32:04.296363Z","iopub.execute_input":"2022-12-11T04:32:04.296764Z","iopub.status.idle":"2022-12-11T04:32:04.308494Z","shell.execute_reply.started":"2022-12-11T04:32:04.296725Z","shell.execute_reply":"2022-12-11T04:32:04.307430Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That is 209 million pairs created in 40 seconds without running out of RAM! 🙂 Not too bad","metadata":{}},{"cell_type":"code","source":"train_pairs.head()","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:32:04.311737Z","iopub.execute_input":"2022-12-11T04:32:04.312934Z","iopub.status.idle":"2022-12-11T04:32:04.334366Z","shell.execute_reply.started":"2022-12-11T04:32:04.312897Z","shell.execute_reply":"2022-12-11T04:32:04.333618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see what is the cardinality of our aids -- we will need this to create the embedding layer.","metadata":{}},{"cell_type":"code","source":"cardinality_aids = max(train_pairs['aid'].max(), train_pairs['aid_next'].max())\ncardinality_aids","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:32:04.335521Z","iopub.execute_input":"2022-12-11T04:32:04.336372Z","iopub.status.idle":"2022-12-11T04:32:04.450165Z","shell.execute_reply.started":"2022-12-11T04:32:04.336337Z","shell.execute_reply":"2022-12-11T04:32:04.449125Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will have up to `1855602` -- that is a lot! But our matrix factorization model will be able to handle this.\n\nLet's construct a `PyTorch` dataset and `dataloader`.","metadata":{}},{"cell_type":"code","source":"import torch\nfrom torch import nn\nfrom torch.utils.data import Dataset, DataLoader\n\nclass ClicksDataset(Dataset):\n    def __init__(self, pairs):\n        self.aid1 = pairs['aid'].to_numpy()\n        self.aid2 = pairs['aid_next'].to_numpy()\n    def __getitem__(self, idx):\n        aid1 = self.aid1[idx]\n        aid2 = self.aid2[idx]\n        return [aid1, aid2]\n    def __len__(self):\n        return len(self.aid1)\n\ntrain_ds = ClicksDataset(train_pairs[:-10_000_000])\nvalid_ds = ClicksDataset(train_pairs[10_000_000:])","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:32:04.451851Z","iopub.execute_input":"2022-12-11T04:32:04.452985Z","iopub.status.idle":"2022-12-11T04:32:04.469608Z","shell.execute_reply.started":"2022-12-11T04:32:04.452935Z","shell.execute_reply":"2022-12-11T04:32:04.468174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let us see how quickly we can iterate over a single epoch with a batch size of `65536`.","metadata":{}},{"cell_type":"code","source":"train_ds = ClicksDataset(train_pairs)\ntrain_dl_pytorch = DataLoader(train_ds, 65536, True, num_workers=2)","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:32:04.470913Z","iopub.execute_input":"2022-12-11T04:32:04.471255Z","iopub.status.idle":"2022-12-11T04:32:04.478209Z","shell.execute_reply.started":"2022-12-11T04:32:04.471226Z","shell.execute_reply":"2022-12-11T04:32:04.477438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nfor batch in train_dl_pytorch:\n    aid1, aid2 = batch[0], batch[1]","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:32:04.479488Z","iopub.execute_input":"2022-12-11T04:32:04.480293Z","iopub.status.idle":"2022-12-11T04:34:20.731301Z","shell.execute_reply.started":"2022-12-11T04:32:04.480235Z","shell.execute_reply":"2022-12-11T04:34:20.730468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Oh dear, that took forever! Mind you, were are not doing anything here, apart from iterating over the dataset for a single epoch (and that is without validation!).\n\nThe reason this is taking so long is that indexing into the the arrays and collating results into batches is very computationally expensive.\n\nThere are ways to work around this but they require writing a lot of code (you could use the iterable-style dataset). And still our solution wouldn't be particularly well optimized.\n\nLet us do something else instead!\n\nWe will use a brand new [Merlin Dataloader](https://github.com/NVIDIA-Merlin/dataloader). It is a library that my team launched just a couple of days ago 🙂\n\nNow this library shines when you have a GPU, which is what you generally want when training DL models. But, alas, Kaggle gives you only 13 GB of RAM on a kernel with a GPU, and that wouldn't allow us to process our dataset!\n\nLet's see how far we can get with CPU only.","metadata":{}},{"cell_type":"code","source":"!pip install merlin-dataloader==0.0.2","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:28:46.018619Z","iopub.execute_input":"2022-12-11T04:28:46.019308Z","iopub.status.idle":"2022-12-11T04:30:10.703915Z","shell.execute_reply.started":"2022-12-11T04:28:46.019203Z","shell.execute_reply":"2022-12-11T04:30:10.702958Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from merlin.loader.torch import Loader ","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:34:20.732806Z","iopub.execute_input":"2022-12-11T04:34:20.734834Z","iopub.status.idle":"2022-12-11T04:34:20.740084Z","shell.execute_reply.started":"2022-12-11T04:34:20.734784Z","shell.execute_reply":"2022-12-11T04:34:20.738968Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can read data directly from the disk -- even better!\n\nLet's write our datasets to disk.","metadata":{}},{"cell_type":"code","source":"train_pairs[:-10_000_000].to_pandas().to_parquet('train_pairs.parquet')\ntrain_pairs[-10_000_000:].to_pandas().to_parquet('valid_pairs.parquet')","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:34:20.743416Z","iopub.execute_input":"2022-12-11T04:34:20.744262Z","iopub.status.idle":"2022-12-11T04:34:29.225191Z","shell.execute_reply.started":"2022-12-11T04:34:20.744221Z","shell.execute_reply":"2022-12-11T04:34:29.224082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from merlin.loader.torch import Loader \nfrom merlin.io import Dataset\n\ntrain_ds = Dataset('train_pairs.parquet')\ntrain_dl_merlin = Loader(train_ds, 65536, True)","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:34:29.226438Z","iopub.execute_input":"2022-12-11T04:34:29.226794Z","iopub.status.idle":"2022-12-11T04:34:30.266400Z","shell.execute_reply.started":"2022-12-11T04:34:29.226760Z","shell.execute_reply":"2022-12-11T04:34:30.265394Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nfor batch, _ in train_dl_merlin:\n    aid1, aid2 = batch['aid'], batch['aid_next']","metadata":{"execution":{"iopub.status.busy":"2022-12-11T04:36:19.765430Z","iopub.execute_input":"2022-12-11T04:36:19.766328Z","iopub.status.idle":"2022-12-11T04:36:51.008486Z","shell.execute_reply.started":"2022-12-11T04:36:19.766288Z","shell.execute_reply":"2022-12-11T04:36:51.007632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That is much better 🙂. Let's train our matrix factorization model!","metadata":{}},{"cell_type":"code","source":"class MatrixFactorization(nn.Module):\n    def __init__(self, n_aids, n_factors):\n        super().__init__()\n        self.aid_factors = nn.Embedding(n_aids, n_factors, sparse=True)\n        \n    def forward(self, aid1, aid2):\n        aid1 = self.aid_factors(aid1)\n        aid2 = self.aid_factors(aid2)\n        \n        return (aid1 * aid2).sum(dim=1)\n    \nclass AverageMeter(object):\n    \"\"\"Computes and stores the average and current value\"\"\"\n    def __init__(self, name, fmt=':f'):\n        self.name = name\n        self.fmt = fmt\n        self.reset()\n\n    def reset(self):\n        self.val = 0\n        self.avg = 0\n        self.sum = 0\n        self.count = 0\n\n    def update(self, val, n=1):\n        self.val = val\n        self.sum += val * n\n        self.count += n\n        self.avg = self.sum / self.count\n\n    def __str__(self):\n        fmtstr = '{name} {val' + self.fmt + '} ({avg' + self.fmt + '})'\n        return fmtstr.format(**self.__dict__)\n\nvalid_ds = Dataset('valid_pairs.parquet')\nvalid_dl_merlin = Loader(valid_ds, 65536, True)","metadata":{"execution":{"iopub.status.busy":"2022-11-18T04:04:04.413613Z","iopub.execute_input":"2022-11-18T04:04:04.414111Z","iopub.status.idle":"2022-11-18T04:04:04.830825Z","shell.execute_reply.started":"2022-11-18T04:04:04.414075Z","shell.execute_reply":"2022-11-18T04:04:04.829486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from torch.optim import SparseAdam\n\nnum_epochs=1\nlr=0.1\n\nmodel = MatrixFactorization(cardinality_aids+1, 32)\noptimizer = SparseAdam(model.parameters(), lr=lr)\ncriterion = nn.BCEWithLogitsLoss()","metadata":{"execution":{"iopub.status.busy":"2022-11-18T04:04:14.877773Z","iopub.execute_input":"2022-11-18T04:04:14.878705Z","iopub.status.idle":"2022-11-18T04:04:15.530201Z","shell.execute_reply.started":"2022-11-18T04:04:14.878655Z","shell.execute_reply":"2022-11-18T04:04:15.528978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nfor epoch in range(num_epochs):\n    for batch, _ in train_dl_merlin:\n        model.train()\n        losses = AverageMeter('Loss', ':.4e')\n            \n        aid1, aid2 = batch['aid'], batch['aid_next']\n        output_pos = model(aid1, aid2)\n        output_neg = model(aid1, aid2[torch.randperm(aid2.shape[0])])\n        \n        output = torch.cat([output_pos, output_neg])\n        targets = torch.cat([torch.ones_like(output_pos), torch.zeros_like(output_pos)])\n        loss = criterion(output, targets)\n        losses.update(loss.item())\n        \n        optimizer.zero_grad()\n        loss.backward()\n        optimizer.step()\n        \n    model.eval()\n    \n    with torch.no_grad():\n        accuracy = AverageMeter('accuracy')\n        for batch, _ in valid_dl_merlin:\n            aid1, aid2 = batch['aid'], batch['aid_next']\n            output_pos = model(aid1, aid2)\n            output_neg = model(aid1, aid2[torch.randperm(aid2.shape[0])])\n            accuracy_batch = torch.cat([output_pos.sigmoid() > 0.5, output_neg.sigmoid() < 0.5]).float().mean()\n            accuracy.update(accuracy_batch, aid1.shape[0])\n            \n    print(f'{epoch+1:02d}: * TrainLoss {losses.avg:.3f}  * Accuracy {accuracy.avg:.3f}')","metadata":{"execution":{"iopub.status.busy":"2022-11-18T04:04:16.046547Z","iopub.execute_input":"2022-11-18T04:04:16.047340Z","iopub.status.idle":"2022-11-18T04:36:13.371266Z","shell.execute_reply.started":"2022-11-18T04:04:16.047285Z","shell.execute_reply":"2022-11-18T04:36:13.369786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's grab the embeddings!","metadata":{}},{"cell_type":"code","source":"embeddings = model.aid_factors.weight.detach().numpy()","metadata":{"execution":{"iopub.status.busy":"2022-11-18T04:39:45.835130Z","iopub.execute_input":"2022-11-18T04:39:45.835623Z","iopub.status.idle":"2022-11-18T04:39:45.842282Z","shell.execute_reply.started":"2022-11-18T04:39:45.835584Z","shell.execute_reply":"2022-11-18T04:39:45.841126Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And construct create the index for approximate nearest neighbor search.","metadata":{}},{"cell_type":"code","source":"%%time\n\nfrom annoy import AnnoyIndex\n\nindex = AnnoyIndex(32, 'euclidean')\nfor i, v in enumerate(embeddings):\n    index.add_item(i, v)\n    \nindex.build(10)","metadata":{"execution":{"iopub.status.busy":"2022-11-18T04:58:54.560248Z","iopub.execute_input":"2022-11-18T04:58:54.560673Z","iopub.status.idle":"2022-11-18T04:59:24.407218Z","shell.execute_reply.started":"2022-11-18T04:58:54.560640Z","shell.execute_reply":"2022-11-18T04:59:24.405675Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now for any `aid`, we can find its nearest neighbor!","metadata":{}},{"cell_type":"code","source":"index.get_nns_by_item(123, 10)","metadata":{"execution":{"iopub.status.busy":"2022-11-18T04:59:24.409450Z","iopub.execute_input":"2022-11-18T04:59:24.409829Z","iopub.status.idle":"2022-11-18T04:59:24.419200Z","shell.execute_reply.started":"2022-11-18T04:59:24.409798Z","shell.execute_reply":"2022-11-18T04:59:24.417561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's create a submission! 🙂","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\n\nfrom collections import defaultdict\n\nsample_sub = pd.read_csv('../input/otto-recommender-system//sample_submission.csv')\n\nsession_types = ['clicks', 'carts', 'orders']\ntest_session_AIDs = test.to_pandas().reset_index(drop=True).groupby('session')['aid'].apply(list)\ntest_session_types = test.to_pandas().reset_index(drop=True).groupby('session')['type'].apply(list)\n\nlabels = []\n\ntype_weight_multipliers = {0: 1, 1: 6, 2: 3}\nfor AIDs, types in zip(test_session_AIDs, test_session_types):\n    if len(AIDs) >= 20:\n        # if we have enough aids (over equals 20) we don't need to look for candidates! we just use the old logic\n        weights=np.logspace(0.1,1,len(AIDs),base=2, endpoint=True)-1\n        aids_temp=defaultdict(lambda: 0)\n        for aid,w,t in zip(AIDs,weights,types): \n            aids_temp[aid]+= w * type_weight_multipliers[t]\n            \n        sorted_aids=[k for k, v in sorted(aids_temp.items(), key=lambda item: -item[1])]\n        labels.append(sorted_aids[:20])\n    else:\n        # here we don't have 20 aids to output -- we will use approximate nearest neighbor search and our embeddings\n        # to generate candidates!\n        AIDs = list(dict.fromkeys(AIDs[::-1]))\n        \n        # let's grab the most recent aid\n        most_recent_aid = AIDs[0]\n        \n        # and look for some neighbors!\n        nns = index.get_nns_by_item(most_recent_aid, 21)[1:]\n                        \n        labels.append((AIDs+nns)[:20])","metadata":{"execution":{"iopub.status.busy":"2022-11-18T04:59:31.828840Z","iopub.execute_input":"2022-11-18T04:59:31.829303Z","iopub.status.idle":"2022-11-18T05:02:14.006532Z","shell.execute_reply.started":"2022-11-18T04:59:31.829269Z","shell.execute_reply":"2022-11-18T05:02:14.004874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's now pull it all together and write to a file,","metadata":{}},{"cell_type":"code","source":"labels_as_strings = [' '.join([str(l) for l in lls]) for lls in labels]\n\npredictions = pd.DataFrame(data={'session_type': test_session_AIDs.index, 'labels': labels_as_strings})\n\nprediction_dfs = []\n\nfor st in session_types:\n    modified_predictions = predictions.copy()\n    modified_predictions.session_type = modified_predictions.session_type.astype('str') + f'_{st}'\n    prediction_dfs.append(modified_predictions)\n\nsubmission = pd.concat(prediction_dfs).reset_index(drop=True)\nsubmission.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2022-11-18T05:02:29.973292Z","iopub.execute_input":"2022-11-18T05:02:29.974763Z","iopub.status.idle":"2022-11-18T05:03:10.037558Z","shell.execute_reply.started":"2022-11-18T05:02:29.974713Z","shell.execute_reply":"2022-11-18T05:03:10.036204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And we are done!\n\n\n**If you like this notebook, please smash the upvote button! Thank you! 😊**\n\nThere are many ways in which this can be expanded:\n* we can train on the GPU\n* we can train for longer\n* maybe we would get better results if we were to filter our train data by type?\n* should we train only on adjacent aids? maybe we should expand the neighborhood we train on\n\nWe can keep asking ourselves many questions like this 🙂 Now we have a framework to start answering them!\n\nThank you for reading! Happy Kaggling! 🙌","metadata":{"execution":{"iopub.status.busy":"2022-11-18T02:49:02.940358Z","iopub.execute_input":"2022-11-18T02:49:02.940858Z","iopub.status.idle":"2022-11-18T02:49:02.973867Z","shell.execute_reply.started":"2022-11-18T02:49:02.940760Z","shell.execute_reply":"2022-11-18T02:49:02.972223Z"}}}]}