MSD-low-memory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Markov Random Fields for Collaborative Filtering (Memory Efficient)\n",
    "\n",
    "This notebook provides a **memory efficient version** in Python 3.7 of the algorithm outlined in the paper \n",
    "\"[Markov Random Fields for Collaborative Filtering](https://arxiv.org/abs/1910.09645)\" \n",
    "at the 33rd Conference on Neural Information Processing Systems (NeurIPS 2019), Vancouver, Canada.\n",
    "\n",
    "For reproducibility, the experiments utilize publicly available [code](https://github.com/dawenl/vae_cf) for pre-processing three popular data-sets and for evaluating the learned model. That code accompanies the paper \"[Variational Autoencoders for Collaborative Filtering](https://arxiv.org/abs/1802.05814)\" by Dawen Liang et al. at The Web Conference 2018. While the code for the Movielens-20M data-set was made publicly available, the code for pre-processing the other two data-sets can easily be obtained by modifying their code as described in their paper.\n",
    "\n",
    "The experiments in the paper (where an AWS instance with 64 GB RAM and 16 vCPUs was used) may be re-run by following these three steps:\n",
    "- Step 1: Pre-processing the data (utilizing the publicly available [code](https://github.com/dawenl/vae_cf))\n",
    "- Step 2: Learning the MRF (this code implements the new algorithm)\n",
    "- Step 3: Evaluation (utilizing the publicly available [code](https://github.com/dawenl/vae_cf))\n",
    "\n",
    "This memory efficient version is modified by Yifei Shen @ Hong Kong University of Science and Technology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Pre-processing the data \n",
    "\n",
    "Utilizing the publicly available [code](https://github.com/dawenl/vae_cf), which is copied below (with kind permission of Dawen Liang):\n",
    " - run their cells 1-26 for data pre-processing \n",
    "      - note that importing matplotlib, seaborn, and tensorflow may not be necessary for our purposes here\n",
    " - run their cells 29-31 for loading the training data\n",
    " \n",
    "Note that the following code is modified as to pre-process the [MSD data-set](https://labrosa.ee.columbia.edu/millionsong/tasteprofile). For pre-processing the [MovieLens-20M data-set](https://grouplens.org/datasets/movielens/20m/), see their original publicly-available [code](https://github.com/dawenl/vae_cf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import sys\n",
    "\n",
    "import numpy as np\n",
    "from scipy import sparse\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "import bottleneck as bn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change to the location of the data\n",
    "DATA_DIR = 'MSD'\n",
    "\n",
    "itemId='songId'   # for MSD data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_data = pd.read_csv(os.path.join(DATA_DIR, 'train_triplets.txt'), sep='\\t', header=None, names=['userId', 'songId', 'playCount'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data splitting procedure\n",
    "- Select 50K users as heldout users, 50K users as validation users, and the rest of the users for training\n",
    "- Use all the items from the training users as item set\n",
    "- For each of both validation and test user, subsample 80% as fold-in data and the rest for prediction "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_count(tp, id):\n",
    "    playcount_groupbyid = tp[[id]].groupby(id, as_index=False)\n",
    "    count = playcount_groupbyid.size()\n",
    "    return count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter_triplets(tp, min_uc=5, min_sc=0):\n",
    "    # Only keep the triplets for items which were clicked on by at least min_sc users. \n",
    "    if min_sc > 0:\n",
    "        itemcount = get_count(tp, itemId)\n",
    "        tp = tp[tp[itemId].isin(itemcount.index[itemcount >= min_sc])]\n",
    "    \n",
    "    # Only keep the triplets for users who clicked on at least min_uc items\n",
    "    # After doing this, some of the items will have less than min_uc users, but should only be a small proportion\n",
    "    if min_uc > 0:\n",
    "        usercount = get_count(tp, 'userId')\n",
    "        tp = tp[tp['userId'].isin(usercount.index[usercount >= min_uc])]\n",
    "    \n",
    "    # Update both usercount and itemcount after filtering\n",
    "    usercount, itemcount = get_count(tp, 'userId'), get_count(tp, itemId) \n",
    "    return tp, usercount, itemcount"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_data, user_activity, item_popularity = filter_triplets(raw_data, min_uc=20, min_sc=200) # for MSD data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After filtering, there are 33633450 watching events from 571355 users and 41140 movies (sparsity: 0.143%)\n"
     ]
    }
   ],
   "source": [
    "sparsity = 1. * raw_data.shape[0] / (user_activity.shape[0] * item_popularity.shape[0])\n",
    "\n",
    "print(\"After filtering, there are %d watching events from %d users and %d movies (sparsity: %.3f%%)\" % \n",
    "      (raw_data.shape[0], user_activity.shape[0], item_popularity.shape[0], sparsity * 100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_uid = user_activity.index\n",
    "\n",
    "np.random.seed(98765)\n",
    "idx_perm = np.random.permutation(unique_uid.size)\n",
    "unique_uid = unique_uid[idx_perm]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create train/validation/test users\n",
    "n_users = unique_uid.size\n",
    "n_heldout_users = 50000 # for MSD data\n",
    "\n",
    "tr_users = unique_uid[:(n_users - n_heldout_users * 2)]\n",
    "vd_users = unique_uid[(n_users - n_heldout_users * 2): (n_users - n_heldout_users)]\n",
    "te_users = unique_uid[(n_users - n_heldout_users):]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_plays = raw_data.loc[raw_data['userId'].isin(tr_users)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_sid = pd.unique(train_plays[itemId])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "show2id = dict((sid, i) for (i, sid) in enumerate(unique_sid))\n",
    "profile2id = dict((pid, i) for (i, pid) in enumerate(unique_uid))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "pro_dir = os.path.join(DATA_DIR, 'pro_sg')\n",
    "\n",
    "if not os.path.exists(pro_dir):\n",
    "    os.makedirs(pro_dir)\n",
    "\n",
    "with open(os.path.join(pro_dir, 'unique_sid.txt'), 'w') as f:\n",
    "    for sid in unique_sid:\n",
    "        f.write('%s\\n' % sid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def split_train_test_proportion(data, test_prop=0.2):\n",
    "    data_grouped_by_user = data.groupby('userId')\n",
    "    tr_list, te_list = list(), list()\n",
    "\n",
    "    np.random.seed(98765)\n",
    "\n",
    "    for i, (_, group) in enumerate(data_grouped_by_user):\n",
    "        n_items_u = len(group)\n",
    "\n",
    "        if n_items_u >= 5:\n",
    "            idx = np.zeros(n_items_u, dtype='bool')\n",
    "            idx[np.random.choice(n_items_u, size=int(test_prop * n_items_u), replace=False).astype('int64')] = True\n",
    "\n",
    "            tr_list.append(group[np.logical_not(idx)])\n",
    "            te_list.append(group[idx])\n",
    "        else:\n",
    "            tr_list.append(group)\n",
    "\n",
    "        if i % 5000 == 0:\n",
    "            print(\"%d users sampled\" % i)\n",
    "            sys.stdout.flush()\n",
    "\n",
    "    data_tr = pd.concat(tr_list)\n",
    "    data_te = pd.concat(te_list)\n",
    "    \n",
    "    return data_tr, data_te"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "vad_plays = raw_data.loc[raw_data['userId'].isin(vd_users)]\n",
    "vad_plays = vad_plays.loc[vad_plays[itemId].isin(unique_sid)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 users sampled\n",
      "5000 users sampled\n",
      "10000 users sampled\n",
      "15000 users sampled\n",
      "20000 users sampled\n",
      "25000 users sampled\n",
      "30000 users sampled\n",
      "35000 users sampled\n",
      "40000 users sampled\n",
      "45000 users sampled\n"
     ]
    }
   ],
   "source": [
    "vad_plays_tr, vad_plays_te = split_train_test_proportion(vad_plays)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_plays = raw_data.loc[raw_data['userId'].isin(te_users)]\n",
    "test_plays = test_plays.loc[test_plays[itemId].isin(unique_sid)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 users sampled\n",
      "5000 users sampled\n",
      "10000 users sampled\n",
      "15000 users sampled\n",
      "20000 users sampled\n",
      "25000 users sampled\n",
      "30000 users sampled\n",
      "35000 users sampled\n",
      "40000 users sampled\n",
      "45000 users sampled\n"
     ]
    }
   ],
   "source": [
    "test_plays_tr, test_plays_te = split_train_test_proportion(test_plays)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save the data into (user_index, item_index) format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def numerize(tp):\n",
    "    uid = list(map(lambda x: profile2id[x], tp['userId']))\n",
    "    sid = list(map(lambda x: show2id[x], tp[itemId]))\n",
    "    return pd.DataFrame(data={'uid': uid, 'sid': sid}, columns=['uid', 'sid'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_data = numerize(train_plays)\n",
    "train_data.to_csv(os.path.join(pro_dir, 'train.csv'), index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "vad_data_tr = numerize(vad_plays_tr)\n",
    "vad_data_tr.to_csv(os.path.join(pro_dir, 'validation_tr.csv'), index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "vad_data_te = numerize(vad_plays_te)\n",
    "vad_data_te.to_csv(os.path.join(pro_dir, 'validation_te.csv'), index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_data_tr = numerize(test_plays_tr)\n",
    "test_data_tr.to_csv(os.path.join(pro_dir, 'test_tr.csv'), index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_data_te = numerize(test_plays_te)\n",
    "test_data_te.to_csv(os.path.join(pro_dir, 'test_te.csv'), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the pre-processed training and validation data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_sid = list()\n",
    "with open(os.path.join(pro_dir, 'unique_sid.txt'), 'r') as f:\n",
    "    for line in f:\n",
    "        unique_sid.append(line.strip())\n",
    "\n",
    "n_items = len(unique_sid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_train_data(csv_file):\n",
    "    tp = pd.read_csv(csv_file)\n",
    "    n_users = tp['uid'].max() + 1\n",
    "\n",
    "    rows, cols = tp['uid'], tp['sid']\n",
    "    data = sparse.csr_matrix((np.ones_like(rows),\n",
    "                             (rows, cols)), dtype='float64',\n",
    "                             shape=(n_users, n_items))\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_data = load_train_data(os.path.join(pro_dir, 'train.csv'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Learning the MRF model (implementation of the new algorithm)\n",
    "Now run the following code and choose to learn \n",
    "- either the dense MRF model \n",
    "- or the sparse MRF model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "from copy import deepcopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyClock:\n",
    "    startTime = time.time()\n",
    "    def tic(self):\n",
    "        self.startTime = time.time()\n",
    "    def toc(self):\n",
    "        secs = time.time() - self.startTime \n",
    "        print(\"... elapsed time: {} min {} sec\".format(int(secs//60), secs%60) )\n",
    "\n",
    "myClock = MyClock()\n",
    "totalClock = MyClock()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.75"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pre-computation of the training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter_XtX(train_data, block_size, thd4mem, thd4comp):\n",
    "    # To obtain and sparsify XtX at the same time to save memory\n",
    "    # block_size (2nd input) and threshold for memory (3rd input) controls the memory usage\n",
    "    # thd4comp is the threshold to control training efficiency\n",
    "    \n",
    "    XtXshape = train_data.shape[1]\n",
    "    userCount = train_data.shape[0]\n",
    "    bs = block_size\n",
    "    blocks = train_data.shape[1]// bs + 1\n",
    "    flag = False\n",
    "    thd = thd4mem\n",
    "    \n",
    "    #normalize data\n",
    "    mu = np.squeeze(np.array(np.sum(train_data, axis=0)))/ userCount \n",
    "    variance_times_userCount = (mu - mu * mu) * userCount\n",
    "    rescaling = np.power(variance_times_userCount, alpha / 2.0) \n",
    "    scaling = 1.0  / rescaling\n",
    "    \n",
    "    #block multiplication\n",
    "    for ii in range(blocks):\n",
    "        for jj in range(blocks):\n",
    "            XtX_tmp = np.asarray(train_data[:,bs*ii : bs*(ii+1)].T.dot(train_data[:,bs*jj : bs*(jj+1)]).todense(), dtype = np.float32)\n",
    "            XtX_tmp -= mu[bs*ii:bs*(ii+1),None] * (mu[bs*jj : bs*(jj+1)]* userCount)\n",
    "            XtX_tmp = scaling[bs*ii:bs*(ii+1),None] * XtX_tmp * scaling[bs*jj : bs*(jj+1)]\n",
    "\n",
    "            # sparsification filter 1 to control memory usage \n",
    "            ix = np.where(np.abs(XtX_tmp) > thd) \n",
    "            XtX_nz = XtX_tmp[ix]\n",
    "\n",
    "            ix = np.array(ix, dtype = 'int32')\n",
    "            ix[0,:] += bs*ii\n",
    "            ix[1,:] += bs*jj\n",
    "\n",
    "            if(flag):\n",
    "                ixs = np.concatenate((ixs, ix), axis = 1)\n",
    "                XtX_nzs = np.concatenate((XtX_nzs, XtX_nz), axis = 0)\n",
    "            else:\n",
    "                ixs = ix\n",
    "                XtX_nzs = XtX_nz\n",
    "                flag = True\n",
    "                \n",
    "    #sparsification filter 2 to control training time of the algorithm\n",
    "    ix2 = np.where(np.abs(XtX_nzs) >= thd4comp)\n",
    "    AA_nzs = XtX_nzs[ix2]\n",
    "    AA_ixs = np.squeeze(ixs[:,ix2])\n",
    "    \n",
    "    print(XtX_nzs.shape, AA_nzs.shape)\n",
    "    XtX = sparse.csc_matrix( (XtX_nzs, ixs), shape=(XtXshape,XtXshape), dtype=np.float32)\n",
    "    AA = sparse.csc_matrix( (AA_nzs, AA_ixs), shape=(XtXshape,XtXshape), dtype=np.float32)\n",
    "    return XtX, rescaling, XtX.diagonal(), AA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(53179992,) (8652486,)\n"
     ]
    }
   ],
   "source": [
    "XtX, rescaling, XtXdiag, AtA = filter_XtX(train_data, 10000, 0.04, 0.11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "ii_diag = np.diag_indices(XtX.shape[0])\n",
    "scaling = 1/rescaling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sparse MRF model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_sparsity_pattern(AtA, maxInColumn):\n",
    "    # this implements section 3.1 in the paper.\n",
    "    \n",
    "    print(\"sparsifying the data-matrix (section 3.1 in the paper) ...\")\n",
    "    myClock.tic()\n",
    "    # apply threshold\n",
    "    #ix = np.where( np.abs(XtX) > threshold)\n",
    "    #AA = sparse.csc_matrix( (XtX[ix], ix), shape=XtX.shape, dtype=np.float32)\n",
    "    AA = AtA\n",
    "    # enforce maxInColumn, see section 3.1 in paper\n",
    "    countInColumns=AA.getnnz(axis=0)\n",
    "    iiList = np.where(countInColumns > maxInColumn)[0]\n",
    "    print(\"    number of items with more than {} entries in column: {}\".format(maxInColumn, len(iiList)) )\n",
    "    for ii in iiList:\n",
    "        jj= AA[:,ii].nonzero()[0]\n",
    "        kk = bn.argpartition(-np.abs(np.asarray(AA[jj,ii].todense()).flatten()), maxInColumn)[maxInColumn:]\n",
    "        AA[  jj[kk], ii ] = 0.0\n",
    "    AA.eliminate_zeros()\n",
    "    print(\"    resulting sparsity of AA: {}\".format( AA.nnz*1.0 / AA.shape[0] / AA.shape[0]) )\n",
    "    myClock.toc()\n",
    "\n",
    "    return AA "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sparse_parameter_estimation(rr, XtX, AA, XtXdiag):\n",
    "    # this implements section 3.2 in the paper\n",
    "\n",
    "    # list L in the paper, sorted by item-counts per column, ties broken by item-popularities as reflected by np.diag(XtX)\n",
    "    AAcountInColumns = AA.getnnz(axis=0)\n",
    "    sortedList=np.argsort(AAcountInColumns+ XtXdiag /2.0/ np.max(XtXdiag)   )[::-1]  \n",
    "\n",
    "    print(\"iterating through steps 1,2, and 4 in section 3.2 of the paper ...\")\n",
    "    myClock.tic()\n",
    "    todoIndicators=np.ones(AAcountInColumns.shape[0])\n",
    "    blockList=[]   # list of blocks. Each block is a list of item-indices, to be processed in step 3 of the paper\n",
    "    for ii in sortedList:\n",
    "        if todoIndicators[ii]==1:\n",
    "            nn, _, vals=sparse.find(AA[:,ii])  # step 1 in paper: set nn contains item ii and its neighbors N\n",
    "            kk=np.argsort(np.abs(vals))[::-1]\n",
    "            nn=nn[kk]\n",
    "            blockList.append(nn) # list of items in the block, to be processed in step 3 below\n",
    "            # remove possibly several items from list L, as determined by parameter rr (r in the paper) \n",
    "            dd_count=max(1,int(np.ceil(len(nn)*rr)))\n",
    "            dd=nn[:dd_count] # set D, see step 2 in the paper\n",
    "            todoIndicators[dd]=0  # step 4 in the paper        \n",
    "    myClock.toc()\n",
    "\n",
    "    print(\"now step 3 in section 3.2 of the paper: iterating ...\")\n",
    "    # now the (possibly heavy) computations of step 3:\n",
    "    # given that steps 1,2,4 are already done, the following for-loop could be implemented in parallel.   \n",
    "    myClock.tic()\n",
    "    BBlist_ix1, BBlist_ix2, BBlist_val = [], [], []\n",
    "    for nn in blockList:\n",
    "        #calculate dense solution for the items in set nn\n",
    "        BBblock=np.linalg.inv( np.array(XtX[np.ix_(nn,nn)].todense()) )\n",
    "        #BBblock=np.linalg.inv( XtX[np.ix_(nn,nn)] )\n",
    "        BBblock/=-np.diag(BBblock)\n",
    "        # determine set D based on parameter rr (r in the paper) \n",
    "        dd_count=max(1,int(np.ceil(len(nn)*rr)))\n",
    "        dd=nn[:dd_count] # set D in paper\n",
    "        # store the solution regarding the items in D\n",
    "        blockix = np.meshgrid(dd,nn)\n",
    "        BBlist_ix1.extend(blockix[1].flatten().tolist())\n",
    "        BBlist_ix2.extend(blockix[0].flatten().tolist())\n",
    "        BBlist_val.extend(BBblock[:,:dd_count].flatten().tolist())\n",
    "    myClock.toc()\n",
    "\n",
    "    print(\"final step: obtaining the sparse matrix BB by averaging the solutions regarding the various sets D ...\")\n",
    "    myClock.tic()\n",
    "    BBsum = sparse.csc_matrix( (BBlist_val,  (BBlist_ix1, BBlist_ix2  )  ), shape=XtX.shape, dtype=np.float32) \n",
    "    BBcnt = sparse.csc_matrix( (np.ones(len(BBlist_ix1), dtype=np.float32),  (BBlist_ix1,BBlist_ix2  )  ), shape=XtX.shape, dtype=np.float32) \n",
    "    b_div= sparse.find(BBcnt)[2]\n",
    "    b_3= sparse.find(BBsum)\n",
    "    BBavg = sparse.csc_matrix( ( b_3[2] / b_div   ,  (b_3[0],b_3[1]  )  ), shape=XtX.shape, dtype=np.float32)\n",
    "    BBavg[ii_diag]=0.0\n",
    "    myClock.toc()\n",
    "\n",
    "    print(\"forcing the sparsity pattern of AA onto BB ...\")\n",
    "    myClock.tic()\n",
    "    BBavg = sparse.csr_matrix( ( np.asarray(BBavg[AA.nonzero()]).flatten(),  AA.nonzero() ), shape=BBavg.shape, dtype=np.float32)\n",
    "    \n",
    "    print(\"    resulting sparsity of learned BB: {}\".format( BBavg.nnz * 1.0 / AA.shape[0] / AA.shape[0]) )\n",
    "    myClock.toc()\n",
    "\n",
    "    return BBavg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sparse_solution(rr, maxInColumn, L2reg):\n",
    "    \n",
    "    # sparsity pattern, see section 3.1 in the paper\n",
    "    XtX[ii_diag] = XtXdiag  \n",
    "    AA = calculate_sparsity_pattern(AtA, maxInColumn)\n",
    "\n",
    "    # parameter-estimation, see section 3.2 in the paper \n",
    "    XtX[ii_diag] = XtXdiag+L2reg \n",
    "    BBsparse = sparse_parameter_estimation(rr, XtX, AA, XtXdiag+L2reg)\n",
    "    \n",
    "    return BBsparse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "training the sparse model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "training the sparse model:\n",
      "\n",
      "sparsifying the data-matrix (section 3.1 in the paper) ...\n",
      "    number of items with more than 1000 entries in column: 911\n",
      "    resulting sparsity of AA: 0.0049730209685130795\n",
      "... elapsed time: 0 min 1.0004370212554932 sec\n",
      "iterating through steps 1,2, and 4 in section 3.2 of the paper ...\n",
      "... elapsed time: 0 min 2.8682525157928467 sec\n",
      "now step 3 in section 3.2 of the paper: iterating ...\n",
      "... elapsed time: 0 min 33.74414944648743 sec\n",
      "final step: obtaining the sparse matrix BB by averaging the solutions regarding the various sets D ...\n",
      "... elapsed time: 0 min 29.50312352180481 sec\n",
      "forcing the sparsity pattern of AA onto BB ...\n",
      "    resulting sparsity of learned BB: 0.0049730209685130795\n",
      "... elapsed time: 0 min 2.2936534881591797 sec\n",
      "\n",
      "total training time (including the time for determining the sparsity-pattern):\n",
      "... elapsed time: 1 min 11.883310794830322 sec\n",
      "\n",
      "re-scaling BB back to the original item-popularities ...\n",
      "... elapsed time: 0 min 0.308868408203125 sec\n"
     ]
    }
   ],
   "source": [
    "maxInColumn = 1000\n",
    "# hyper-parameter r in the paper, which determines the trade-off between approximation-accuracy and training-time\n",
    "rr = 0.1\n",
    "# L2 norm regularization\n",
    "L2reg = 1.0 \n",
    "\n",
    "print(\"training the sparse model:\\n\")\n",
    "totalClock.tic()\n",
    "BBsparse = sparse_solution(rr, maxInColumn, L2reg)\n",
    "print(\"\\ntotal training time (including the time for determining the sparsity-pattern):\")\n",
    "totalClock.toc()\n",
    "\n",
    "print(\"\\nre-scaling BB back to the original item-popularities ...\")\n",
    "# assuming that mu.T.dot(BB) == mu, see Appendix in paper\n",
    "myClock.tic()\n",
    "BBsparse=sparse.diags(scaling).dot(BBsparse).dot(sparse.diags(rescaling))\n",
    "myClock.toc()\n",
    "\n",
    "#print(\"\\nfor the evaluation below: converting the sparse model into a dense-matrix-representation ...\")\n",
    "#myClock.tic()\n",
    "#BB = np.asarray(BBsparse.todense(), dtype=np.float32) \n",
    "#myClock.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Evaluating the MRF model \n",
    "\n",
    "Utilizing the publicly available [code](https://github.com/dawenl/vae_cf), which is copied below (with kind permission of Dawen Liang):\n",
    "\n",
    " - run their cell 32 for loading the test data\n",
    " - run their cells 35 and 36 for the ranking metrics (for later use in evaluation)\n",
    " - run their cells 45 and 46\n",
    " - modify and run their cell 50:\n",
    "    - remove 2 lines: the one that starts with ```with``` and the line below\n",
    "    - remove the indentation of the line that starts with ```for```\n",
    "    - modify the line that starts with ```pred_val``` as follows: ```pred_val = X.dot(BB)```\n",
    "        \n",
    " - run their cell 51\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_tr_te_data(csv_file_tr, csv_file_te):\n",
    "    tp_tr = pd.read_csv(csv_file_tr)\n",
    "    tp_te = pd.read_csv(csv_file_te)\n",
    "\n",
    "    start_idx = min(tp_tr['uid'].min(), tp_te['uid'].min())\n",
    "    end_idx = max(tp_tr['uid'].max(), tp_te['uid'].max())\n",
    "\n",
    "    rows_tr, cols_tr = tp_tr['uid'] - start_idx, tp_tr['sid']\n",
    "    rows_te, cols_te = tp_te['uid'] - start_idx, tp_te['sid']\n",
    "\n",
    "    data_tr = sparse.csr_matrix((np.ones_like(rows_tr),\n",
    "                             (rows_tr, cols_tr)), dtype='float64', shape=(end_idx - start_idx + 1, n_items))\n",
    "    data_te = sparse.csr_matrix((np.ones_like(rows_te),\n",
    "                             (rows_te, cols_te)), dtype='float64', shape=(end_idx - start_idx + 1, n_items))\n",
    "    return data_tr, data_te"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def NDCG_binary_at_k_batch(X_pred, heldout_batch, k=100):\n",
    "    '''\n",
    "    normalized discounted cumulative gain@k for binary relevance\n",
    "    ASSUMPTIONS: all the 0's in heldout_data indicate 0 relevance\n",
    "    '''\n",
    "    batch_users = X_pred.shape[0]\n",
    "    idx_topk_part = bn.argpartition(-X_pred, k, axis=1)\n",
    "    topk_part = X_pred[np.arange(batch_users)[:, np.newaxis],\n",
    "                       idx_topk_part[:, :k]]\n",
    "    idx_part = np.argsort(-topk_part, axis=1)\n",
    "    # X_pred[np.arange(batch_users)[:, np.newaxis], idx_topk] is the sorted\n",
    "    # topk predicted score\n",
    "    idx_topk = idx_topk_part[np.arange(batch_users)[:, np.newaxis], idx_part]\n",
    "    # build the discount template\n",
    "    tp = 1. / np.log2(np.arange(2, k + 2))\n",
    "\n",
    "    DCG = (heldout_batch[np.arange(batch_users)[:, np.newaxis],\n",
    "                         idx_topk].toarray() * tp).sum(axis=1)\n",
    "    IDCG = np.array([(tp[:min(n, k)]).sum()\n",
    "                     for n in heldout_batch.getnnz(axis=1)])\n",
    "    return DCG / IDCG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Recall_at_k_batch(X_pred, heldout_batch, k=100):\n",
    "    batch_users = X_pred.shape[0]\n",
    "\n",
    "    idx = bn.argpartition(-X_pred, k, axis=1)\n",
    "    X_pred_binary = np.zeros_like(X_pred, dtype=bool)\n",
    "    X_pred_binary[np.arange(batch_users)[:, np.newaxis], idx[:, :k]] = True\n",
    "\n",
    "    X_true_binary = (heldout_batch > 0).toarray()\n",
    "    tmp = (np.logical_and(X_true_binary, X_pred_binary).sum(axis=1)).astype(\n",
    "        np.float32)\n",
    "    recall = tmp / np.minimum(k, X_true_binary.sum(axis=1))\n",
    "    return recall"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the test data and compute test metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_data_tr, test_data_te = load_tr_te_data(\n",
    "    os.path.join(pro_dir, 'test_tr.csv'),\n",
    "    os.path.join(pro_dir, 'test_te.csv'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_test = test_data_tr.shape[0]\n",
    "idxlist_test = range(N_test)\n",
    "\n",
    "batch_size_test = 2000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "n100_list, r20_list, r50_list = [], [], []\n",
    "\n",
    "\n",
    "for bnum, st_idx in enumerate(range(0, N_test, batch_size_test)):\n",
    "    end_idx = min(st_idx + batch_size_test, N_test)\n",
    "    X = test_data_tr[idxlist_test[st_idx:end_idx]]\n",
    "\n",
    "    #if sparse.isspmatrix(X):\n",
    "    #    X = X.toarray()\n",
    "    #X = X.astype('float32')\n",
    "\n",
    "    pred_val = np.array(X.dot(BBsparse).todense())\n",
    "    # exclude examples from training and validation (if any)\n",
    "    pred_val[X.nonzero()] = -np.inf\n",
    "    n100_list.append(NDCG_binary_at_k_batch(pred_val, test_data_te[idxlist_test[st_idx:end_idx]], k=100))\n",
    "    r20_list.append(Recall_at_k_batch(pred_val, test_data_te[idxlist_test[st_idx:end_idx]], k=20))\n",
    "    r50_list.append(Recall_at_k_batch(pred_val, test_data_te[idxlist_test[st_idx:end_idx]], k=50))\n",
    "\n",
    "n100_list = np.concatenate(n100_list)\n",
    "r20_list = np.concatenate(r20_list)\n",
    "r50_list = np.concatenate(r50_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test NDCG@100=0.38659 (0.00110)\n",
      "Test Recall@20=0.33094 (0.00113)\n",
      "Test Recall@50=0.42356 (0.00118)\n"
     ]
    }
   ],
   "source": [
    "print(\"Test NDCG@100=%.5f (%.5f)\" % (np.mean(n100_list), np.std(n100_list) / np.sqrt(len(n100_list))))\n",
    "print(\"Test Recall@20=%.5f (%.5f)\" % (np.mean(r20_list), np.std(r20_list) / np.sqrt(len(r20_list))))\n",
    "print(\"Test Recall@50=%.5f (%.5f)\" % (np.mean(r50_list), np.std(r50_list) / np.sqrt(len(r50_list))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... accuracy of the sparse approximation (with sparsity 0.1% and parameter r=0.5)"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}