notebooks/04_model_select_and_optimize/tuning_spark_als.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<i>Copyright (c) Microsoft Corporation. All rights reserved.</i>\n",
    "\n",
    "<i>Licensed under the MIT License.</i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hyperparameter tuning (Spark based recommender)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hyperparameter tuning for Spark based recommender algorithm is important to select a model with the optimal performance. This notebook introduces good practices in performing hyperparameter tuning for building recommender models with the utility functions provided in the [Microsoft/Recommenders](https://github.com/Microsoft/Recommenders.git) repository.\n",
    "\n",
    "Three different approaches are introduced and comparatively studied.\n",
    "* Spark native/custom constructs (`ParamGridBuilder`, `TrainValidationSplit`).\n",
    "* `hyperopt` package with Tree of Parzen Estimator algorithm. \n",
    "* Brute-force random search of parameter values sampled with pre-defined space. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0 Global settings and import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "System version: 3.6.0 | packaged by conda-forge | (default, Feb  9 2017, 14:36:55) \n",
      "[GCC 4.8.2 20140120 (Red Hat 4.8.2-15)]\n",
      "Pandas version: 0.23.4\n",
      "PySpark version: 2.3.1\n"
     ]
    }
   ],
   "source": [
    "# set the environment path to find Recommenders\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "sys.path.append(\"../../\")\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "import pyspark\n",
    "from pyspark.sql import SparkSession\n",
    "import pyspark.sql.functions as F\n",
    "from pyspark import SparkContext, SparkConf\n",
    "from pyspark.ml.recommendation import ALS\n",
    "from pyspark.ml.tuning import ParamGridBuilder, TrainValidationSplit\n",
    "from pyspark.ml.evaluation import Evaluator, RegressionEvaluator\n",
    "from pyspark.ml.pipeline import Estimator, Model\n",
    "from pyspark import keyword_only  \n",
    "from pyspark.ml import Transformer\n",
    "from pyspark.ml.param.shared import *\n",
    "from pyspark.ml.util import *\n",
    "from pyspark.mllib.evaluation import RegressionMetrics, RankingMetrics\n",
    "from pyspark.sql.types import ArrayType, IntegerType, StringType\n",
    "\n",
    "from hyperopt import fmin, tpe, hp, STATUS_OK, Trials\n",
    "from hyperopt.pyll.base import scope\n",
    "from hyperopt.pyll.stochastic import sample\n",
    "\n",
    "from reco_utils.common.spark_utils import start_or_get_spark\n",
    "from reco_utils.evaluation.spark_evaluation import SparkRankingEvaluation, SparkRatingEvaluation\n",
    "from reco_utils.dataset.movielens import load_spark_df\n",
    "from reco_utils.dataset.spark_splitters import spark_random_split\n",
    "\n",
    "print(\"System version: {}\".format(sys.version))\n",
    "print(\"Pandas version: {}\".format(pd.__version__))\n",
    "print(\"PySpark version: {}\".format(pyspark.__version__))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "NUMBER_CORES = 1\n",
    "NUMBER_ITERATIONS = 25\n",
    "\n",
    "COL_USER = \"userID\"\n",
    "COL_ITEM = \"itemID\"\n",
    "COL_TIMESTAMP = \"timestamp\"\n",
    "COL_RATING = \"rating\"\n",
    "COL_PREDICTION = \"prediction\"\n",
    "\n",
    "HEADER = {\n",
    "    \"col_user\": COL_USER,\n",
    "    \"col_item\": COL_ITEM,\n",
    "    \"col_rating\": COL_RATING,\n",
    "    \"col_prediction\": COL_PREDICTION,\n",
    "}\n",
    "\n",
    "HEADER_ALS = {\n",
    "    \"userCol\": COL_USER,\n",
    "    \"itemCol\": COL_ITEM,\n",
    "    \"ratingCol\": COL_RATING\n",
    "}\n",
    "\n",
    "RANK = [10, 15, 20, 30, 40]\n",
    "REG = [ 0.1, 0.01, 0.001, 0.0001, 0.00001]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Spark session is created. Note in this case, to study the running time for different approaches, the Spark session in local mode uses only one core for running. This eliminates the impact of parallelization of parameter tuning. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "spark = start_or_get_spark(url=\"local[{}]\".format(NUMBER_CORES))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Movielens 100k dataset is used for running the demonstration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = load_spark_df(spark, size='100k', header=(COL_USER, COL_ITEM, COL_RATING))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset is split into 3 subsets randomly with a given split ratio. The hyperparameter tuning is performed on the training and the validating data, and then the optimal recommender selected is evaluated on the testing dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "train, valid, test = spark_random_split(data, ratio=[3, 1, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Hyper parameter tuning with Azure Machine Learning Services"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `hyperdrive` module in the [Azure Machine Learning Services](https://azure.microsoft.com/en-us/services/machine-learning-service/) runs [hyperparameter tuning and optimizing for machine learning model selection](https://docs.microsoft.com/en-us/azure/machine-learning/service/how-to-tune-hyperparameters). At the moment, the service supports running hyperparameter tuning on heterogenous computing targets such as cluster of commodity compute nodes with or without GPU devices (see detailed documentation [here](https://docs.microsoft.com/en-us/azure/machine-learning/service/how-to-set-up-training-targets)). It is feasible to run parameter tuning on a cluster of VM nodes. In this case, the service containerizes individual and independent Spark session on each node of the cluster to run the parameter tuning job in parallel, instead of inside a single Spark session where the training is executed in a distributed manner.  \n",
    "\n",
    "Detailed instructions of tuning hyperparameter of non-Spark workloads by using Azure Machine Learning Services can be found in [this](./hypertune_aml_wide_and_deep_quickstart.ipynb) notebook. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Hyper parameter tuning with Spark ML constructs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Spark native construct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spark ML lib implements modules such as `CrossValidator` and `TrainValidationSplit` for tuning hyperparameters (see [here](https://spark.apache.org/docs/2.2.0/ml-tuning.html)). However, by default, it does not support custom machine learning algorithms, data splitting methods, and evaluation metrics, like what are offered as utility functions in the Recommenders repository. \n",
    "\n",
    "For example, the Spark native constuct can be used for tuning a recommender against the `rmse` metric which is one of the available regression metrics in Spark."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, a Spark ALS object needs to be created. In this case, for illustration purpose, it is an ALS model object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE the parameters of interest, rank and regParam, are left unset, \n",
    "# because their values will be assigned in the parameter grid builder.\n",
    "als = ALS(\n",
    "    maxIter=15,\n",
    "    implicitPrefs=False,\n",
    "    alpha=0.1,\n",
    "    coldStartStrategy='drop',\n",
    "    nonnegative=False,\n",
    "    **HEADER_ALS\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, a parameter grid can be defined as follows. Without loss of generity, only `rank` and `regParam` are considered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "paramGrid = ParamGridBuilder() \\\n",
    "    .addGrid(als.rank, RANK) \\\n",
    "    .addGrid(als.regParam, REG) \\\n",
    "    .build()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the settings above, a `TrainValidationSplit` constructor can be created for fitting the best model in the given parameter range. In this case, the `RegressionEvaluator` is using `RMSE`, by default, as an evaluation metric. \n",
    "\n",
    "Since the data splitter is embedded in the `TrainValidationSplit` object, to make sure the splitting ratio is consistent across different approaches, the split ratio is set to be 0.75 and in the model training the training dataset and validating dataset are combined. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "tvs = TrainValidationSplit(\n",
    "    estimator=als,\n",
    "    estimatorParamMaps=paramGrid,\n",
    "    # A regression evaluation method is used. \n",
    "    evaluator=RegressionEvaluator(labelCol='rating'),\n",
    "    # 75% of the data will be used for training, 25% for validation.\n",
    "    # NOTE here the splitting is random. The Spark splitting utilities (e.g. chrono splitter)\n",
    "    # are therefore not available here. \n",
    "    trainRatio=0.75\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_start = time.time()\n",
    "\n",
    "# Run TrainValidationSplit, and choose the best set of parameters.\n",
    "# NOTE train and valid is union because in Spark TrainValidationSplit does splitting by itself.\n",
    "model = tvs.fit(train.union(valid))\n",
    "\n",
    "time_spark = time.time() - time_start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model parameters in the grid and the best metrics can be then returned. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Run 0:\n",
      "\tValidation Metric: 0.9482728210259341\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 10\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.1\n",
      "Run 1:\n",
      "\tValidation Metric: 1.2067929009990404\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 10\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.01\n",
      "Run 2:\n",
      "\tValidation Metric: 1.5497272849660675\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 10\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.001\n",
      "Run 3:\n",
      "\tValidation Metric: 1.9726082206094249\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 10\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.0001\n",
      "Run 4:\n",
      "\tValidation Metric: 2.7411419481935124\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 10\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 1e-05\n",
      "Run 5:\n",
      "\tValidation Metric: 0.9490758788255582\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 15\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.1\n",
      "Run 6:\n",
      "\tValidation Metric: 1.2512754564370978\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 15\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.01\n",
      "Run 7:\n",
      "\tValidation Metric: 1.6586340127334998\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 15\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.001\n",
      "Run 8:\n",
      "\tValidation Metric: 2.30626695924491\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 15\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.0001\n",
      "Run 9:\n",
      "\tValidation Metric: 3.3726258779220606\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 15\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 1e-05\n",
      "Run 10:\n",
      "\tValidation Metric: 0.9469224333030561\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 20\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.1\n",
      "Run 11:\n",
      "\tValidation Metric: 1.2807538324207055\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 20\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.01\n",
      "Run 12:\n",
      "\tValidation Metric: 1.69870985430111\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 20\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.001\n",
      "Run 13:\n",
      "\tValidation Metric: 2.3299551882131095\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 20\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.0001\n",
      "Run 14:\n",
      "\tValidation Metric: 2.9949845326363054\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 20\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 1e-05\n",
      "Run 15:\n",
      "\tValidation Metric: 0.9455162431620866\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 30\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.1\n",
      "Run 16:\n",
      "\tValidation Metric: 1.3070897396257575\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 30\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.01\n",
      "Run 17:\n",
      "\tValidation Metric: 1.6655149914634848\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 30\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.001\n",
      "Run 18:\n",
      "\tValidation Metric: 2.153150330244491\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 30\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.0001\n",
      "Run 19:\n",
      "\tValidation Metric: 2.9370601059120793\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 30\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 1e-05\n",
      "Run 20:\n",
      "\tValidation Metric: 0.9433648191344912\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 40\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.1\n",
      "Run 21:\n",
      "\tValidation Metric: 1.2768443759833505\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 40\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.01\n",
      "Run 22:\n",
      "\tValidation Metric: 1.644869524017942\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 40\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.001\n",
      "Run 23:\n",
      "\tValidation Metric: 2.0666533742745883\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 40\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 0.0001\n",
      "Run 24:\n",
      "\tValidation Metric: 2.61814929846137\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='rank', doc='rank of the factorization'): 40\n",
      "\tParam(parent='ALS_4b268dbbfa8855fab501', name='regParam', doc='regularization parameter (>= 0).'): 1e-05\n"
     ]
    }
   ],
   "source": [
    "for idx, item in enumerate(model.getEstimatorParamMaps()):\n",
    "    print('Run {}:'.format(idx))\n",
    "    print('\\tValidation Metric: {}'.format(model.validationMetrics[idx]))\n",
    "    for key, value in item.items():\n",
    "        print('\\t{0}: {1}'.format(repr(key), value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.9482728210259341,\n",
       " 1.2067929009990404,\n",
       " 1.5497272849660675,\n",
       " 1.9726082206094249,\n",
       " 2.7411419481935124,\n",
       " 0.9490758788255582,\n",
       " 1.2512754564370978,\n",
       " 1.6586340127334998,\n",
       " 2.30626695924491,\n",
       " 3.3726258779220606,\n",
       " 0.9469224333030561,\n",
       " 1.2807538324207055,\n",
       " 1.69870985430111,\n",
       " 2.3299551882131095,\n",
       " 2.9949845326363054,\n",
       " 0.9455162431620866,\n",
       " 1.3070897396257575,\n",
       " 1.6655149914634848,\n",
       " 2.153150330244491,\n",
       " 2.9370601059120793,\n",
       " 0.9433648191344912,\n",
       " 1.2768443759833505,\n",
       " 1.644869524017942,\n",
       " 2.0666533742745883,\n",
       " 2.61814929846137]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.validationMetrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the best model, just do"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_best_spark = model.bestModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Custom `Estimator`, `Transformer`, and `Evaluator` for Spark ALS\n",
    "\n",
    "One can also customize Spark modules to allow tuning hyperparameters for a desired model and evaluation metric, given that the native Spark ALS does not allow tuning hyperparameters for ranking metrics such as precision@k, recall@k, etc. This can be done by creating custom `Estimator`, `Transformer` and `Evaluator`. The benefit is that, after the customization, the tuning process can make use of `trainValidSplit` directly, which distributes the tuning in a Spark session."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Customized `Estimator` and `Transformer` for top k recommender based on Spark ALS\n",
    "\n",
    "The following shows how to implement a PySpark `Estimator` and `Transfomer` for recommending top k items from ALS model. The latter generates top k recommendations from the model object. Both of the two are designed by following the protocol of Spark APIs, to make sure that they can be run with the hyperparameter tuning constructs in Spark."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ALSTopK(\n",
    "    ALS,\n",
    "    Estimator,\n",
    "    HasInputCol,\n",
    "    HasPredictionCol\n",
    "):    \n",
    "    rank = Param(Params._dummy(), \"rank\", \"rank of the factorization\",\n",
    "                 typeConverter=TypeConverters.toInt)\n",
    "    numUserBlocks = Param(Params._dummy(), \"numUserBlocks\", \"number of user blocks\",\n",
    "                          typeConverter=TypeConverters.toInt)\n",
    "    numItemBlocks = Param(Params._dummy(), \"numItemBlocks\", \"number of item blocks\",\n",
    "                          typeConverter=TypeConverters.toInt)\n",
    "    implicitPrefs = Param(Params._dummy(), \"implicitPrefs\", \"whether to use implicit preference\",\n",
    "                          typeConverter=TypeConverters.toBoolean)\n",
    "    alpha = Param(Params._dummy(), \"alpha\", \"alpha for implicit preference\",\n",
    "                  typeConverter=TypeConverters.toFloat)\n",
    "    userCol = Param(Params._dummy(), \"userCol\", \"column name for user ids. Ids must be within \" +\n",
    "                    \"the integer value range.\", typeConverter=TypeConverters.toString)\n",
    "    itemCol = Param(Params._dummy(), \"itemCol\", \"column name for item ids. Ids must be within \" +\n",
    "                    \"the integer value range.\", typeConverter=TypeConverters.toString)\n",
    "    ratingCol = Param(Params._dummy(), \"ratingCol\", \"column name for ratings\",\n",
    "                      typeConverter=TypeConverters.toString)\n",
    "    nonnegative = Param(Params._dummy(), \"nonnegative\",\n",
    "                        \"whether to use nonnegative constraint for least squares\",\n",
    "                        typeConverter=TypeConverters.toBoolean)\n",
    "    intermediateStorageLevel = Param(Params._dummy(), \"intermediateStorageLevel\",\n",
    "                                     \"StorageLevel for intermediate datasets. Cannot be 'NONE'.\",\n",
    "                                     typeConverter=TypeConverters.toString)\n",
    "    finalStorageLevel = Param(Params._dummy(), \"finalStorageLevel\",\n",
    "                              \"StorageLevel for ALS model factors.\",\n",
    "                              typeConverter=TypeConverters.toString)\n",
    "    coldStartStrategy = Param(Params._dummy(), \"coldStartStrategy\", \"strategy for dealing with \" +\n",
    "                              \"unknown or new users/items at prediction time. This may be useful \" +\n",
    "                              \"in cross-validation or production scenarios, for handling \" +\n",
    "                              \"user/item ids the model has not seen in the training data. \" +\n",
    "                              \"Supported values: 'nan', 'drop'.\",\n",
    "                              typeConverter=TypeConverters.toString)\n",
    "\n",
    "    @keyword_only\n",
    "    def __init__(\n",
    "        self,\n",
    "        rank=10, maxIter=10, regParam=0.1, numUserBlocks=10, numItemBlocks=10,\n",
    "        implicitPrefs=False, alpha=1.0, userCol=\"user\", itemCol=\"item\", seed=None, k=10,\n",
    "        ratingCol=\"rating\", nonnegative=False, checkpointInterval=10,\n",
    "        intermediateStorageLevel=\"MEMORY_AND_DISK\",\n",
    "        finalStorageLevel=\"MEMORY_AND_DISK\", coldStartStrategy=\"nan\"\n",
    "    ):\n",
    "        super(ALS, self).__init__()\n",
    "        self._java_obj = self._new_java_obj(\"org.apache.spark.ml.recommendation.ALS\", self.uid)\n",
    "        self._setDefault(rank=10, maxIter=10, regParam=0.1, numUserBlocks=10, numItemBlocks=10,\n",
    "                         implicitPrefs=False, alpha=1.0, userCol=\"user\", itemCol=\"item\",\n",
    "                         ratingCol=\"rating\", nonnegative=False, checkpointInterval=10,\n",
    "                         intermediateStorageLevel=\"MEMORY_AND_DISK\",\n",
    "                         finalStorageLevel=\"MEMORY_AND_DISK\", coldStartStrategy=\"nan\")\n",
    "\n",
    "        kwargs = self._input_kwargs \n",
    "        kwargs = {x: kwargs[x] for x in kwargs if x not in {'k'}}\n",
    "        self.setParams(**kwargs)\n",
    "        \n",
    "        # The manually added parameter is not present in ALS Java implementation. \n",
    "        self.k = k\n",
    "        \n",
    "    def setRank(self, value):\n",
    "        \"\"\"\n",
    "        Sets the value of :py:attr:`rank`.\n",
    "        \"\"\"\n",
    "        return self._set(rank=value)\n",
    "\n",
    "    def getRank(self):\n",
    "        \"\"\"\n",
    "        Gets the value of rank or its default value.\n",
    "        \"\"\"\n",
    "        return self.getOrDefault(self.rank)\n",
    "    \n",
    "    def setParams(self, rank=10, maxIter=10, regParam=0.1, numUserBlocks=10, numItemBlocks=10,\n",
    "                  implicitPrefs=False, alpha=1.0, userCol=\"user\", itemCol=\"item\", seed=None,\n",
    "                  ratingCol=\"rating\", nonnegative=False, checkpointInterval=10,\n",
    "                  intermediateStorageLevel=\"MEMORY_AND_DISK\",\n",
    "                  finalStorageLevel=\"MEMORY_AND_DISK\", coldStartStrategy=\"nan\"):\n",
    "        \"\"\"\n",
    "        setParams(self, rank=10, maxIter=10, regParam=0.1, numUserBlocks=10, numItemBlocks=10, \\\n",
    "                 implicitPrefs=False, alpha=1.0, userCol=\"user\", itemCol=\"item\", seed=None, \\\n",
    "                 ratingCol=\"rating\", nonnegative=False, checkpointInterval=10, \\\n",
    "                 intermediateStorageLevel=\"MEMORY_AND_DISK\", \\\n",
    "                 finalStorageLevel=\"MEMORY_AND_DISK\", coldStartStrategy=\"nan\")\n",
    "        Sets params for ALS.\n",
    "        \"\"\"\n",
    "        kwargs = self._input_kwargs\n",
    "        kwargs = {x: kwargs[x] for x in kwargs if x not in {'k'}}\n",
    "        return self._set(**kwargs)\n",
    "        \n",
    "    def _fit(self, dataset):\n",
    "        kwargs = self._input_kwargs    \n",
    "        # Exclude k as it is not a parameter for ALS.\n",
    "        kwargs = {x: kwargs[x] for x in kwargs if x not in {'k'}}\n",
    "        kwargs['rank'] = self.getRank()\n",
    "        kwargs['regParam'] = self.getOrDefault(self.regParam)\n",
    "        als = ALS(\n",
    "            **kwargs\n",
    "        )\n",
    "        als_model = als.fit(dataset)\n",
    "        \n",
    "        user_col = kwargs['userCol']\n",
    "        item_col = kwargs['itemCol']\n",
    "        \n",
    "        k = self.k\n",
    "                     \n",
    "        topk_model = ALSTopKModel()\n",
    "        topk_model.setParams(\n",
    "            als_model,\n",
    "            user_col, \n",
    "            item_col, \n",
    "            k\n",
    "        )\n",
    "        \n",
    "        return topk_model\n",
    "    \n",
    "    \n",
    "class ALSTopKModel(\n",
    "    Model,\n",
    "    HasInputCol,\n",
    "    HasPredictionCol,\n",
    "    HasLabelCol\n",
    "):    \n",
    "    def setParams(self, model, userCol, itemCol, k):\n",
    "        self.model = model\n",
    "        self.userCol = userCol\n",
    "        self.itemCol = itemCol\n",
    "        self.k = k\n",
    "    \n",
    "    def _transform(self, dataset):\n",
    "        predictionCol = self.getPredictionCol()\n",
    "        labelCol = self.getLabelCol()\n",
    "        \n",
    "        users = dataset.select(self.userCol).distinct()\n",
    "        topk_recommendation = self.model.recommendForUserSubset(users, self.k)  \n",
    "        \n",
    "        extract_value = F.udf((lambda x: [y[0] for y in x]), ArrayType(IntegerType()))\n",
    "        topk_recommendation = topk_recommendation.withColumn(predictionCol, extract_value(F.col(\"recommendations\")))        \n",
    "        \n",
    "        dataset = (\n",
    "            dataset\n",
    "            .groupBy(self.userCol)\n",
    "            .agg(F.collect_list(F.col(self.itemCol)).alias(labelCol))\n",
    "        )\n",
    "            \n",
    "        topk_recommendation_all = dataset.join(\n",
    "            topk_recommendation, \n",
    "            on=self.userCol,\n",
    "            how=\"outer\"\n",
    "        )\n",
    "        \n",
    "        return topk_recommendation_all.select(self.userCol, labelCol, predictionCol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Customized precision@k evaluation metric\n",
    "\n",
    "In addition to the custom `Estimator` and `Transformer`, it may also be desired to customize an `Evaluator` to allow \"beyond-rating\" metrics. The codes as following illustrates a precision@k evaluator. Other types of evaluators can be developed in a similar way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a custom Evaulator. Here precision@k is used.\n",
    "class PrecisionAtKEvaluator(Evaluator):\n",
    "\n",
    "    def __init__(self, predictionCol=\"prediction\", labelCol=\"label\", k=10):\n",
    "        self.predictionCol = predictionCol\n",
    "        self.labelCol = labelCol\n",
    "        self.k = k\n",
    "\n",
    "    def _evaluate(self, dataset):\n",
    "        \"\"\"\n",
    "        Returns a random number. \n",
    "        Implement here the true metric\n",
    "        \"\"\"      \n",
    "        # Drop Nulls.\n",
    "        dataset = dataset.na.drop()\n",
    "        metrics = RankingMetrics(dataset.select(self.predictionCol, self.labelCol).rdd)\n",
    "        return metrics.precisionAt(self.k)\n",
    "\n",
    "    def isLargerBetter(self):\n",
    "        return True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then a new ALS top-k recommender can be created, and the Spark native construct, `TrainValidationSplit` module, can be used to find the optimal model w.r.t the precision@k metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "alstopk = ALSTopK(\n",
    "    userCol=COL_USER,\n",
    "    itemCol=COL_ITEM,\n",
    "    ratingCol=COL_RATING,\n",
    "    k=10\n",
    ")\n",
    "\n",
    "# Here for illustration purpose, a small grid is used.\n",
    "paramGrid = ParamGridBuilder() \\\n",
    "    .addGrid(alstopk.rank, [10, 20]) \\\n",
    "    .addGrid(alstopk.regParam, [0.1, 0.01]) \\\n",
    "    .build()\n",
    "\n",
    "tvs = TrainValidationSplit(\n",
    "    estimator=alstopk,\n",
    "    estimatorParamMaps=paramGrid,\n",
    "    # A regression evaluation method is used. \n",
    "    evaluator=PrecisionAtKEvaluator(),\n",
    "    # 75% of the data will be used for training, 25% for validation.\n",
    "    # NOTE here the splitting is random. The Spark splitting utilities (e.g. chrono splitter)\n",
    "    # are therefore not available here. \n",
    "    trainRatio=0.75\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[{Param(parent='ALSTopK_4e2995821d29445956ec', name='rank', doc='rank of the factorization'): 10,\n",
       "  Param(parent='ALSTopK_4e2995821d29445956ec', name='regParam', doc='regularization parameter (>= 0).'): 0.1},\n",
       " {Param(parent='ALSTopK_4e2995821d29445956ec', name='rank', doc='rank of the factorization'): 10,\n",
       "  Param(parent='ALSTopK_4e2995821d29445956ec', name='regParam', doc='regularization parameter (>= 0).'): 0.01},\n",
       " {Param(parent='ALSTopK_4e2995821d29445956ec', name='rank', doc='rank of the factorization'): 20,\n",
       "  Param(parent='ALSTopK_4e2995821d29445956ec', name='regParam', doc='regularization parameter (>= 0).'): 0.1},\n",
       " {Param(parent='ALSTopK_4e2995821d29445956ec', name='rank', doc='rank of the factorization'): 20,\n",
       "  Param(parent='ALSTopK_4e2995821d29445956ec', name='regParam', doc='regularization parameter (>= 0).'): 0.01}]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Run TrainValidationSplit, and choose the best set of parameters.\n",
    "# NOTE train and valid is union because in Spark TrainValidationSplit does splitting by itself.\n",
    "model_precision = tvs.fit(train.union(valid))\n",
    "\n",
    "model_precision.getEstimatorParamMaps()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def best_param(model, is_larger_better=True):\n",
    "    if is_larger_better:\n",
    "        best_metric = max(model.validationMetrics)\n",
    "    else:\n",
    "        best_metric = min(model.validationMetrics)\n",
    "        \n",
    "    parameters = model.getEstimatorParamMaps()[model.validationMetrics.index(best_metric)]\n",
    "     \n",
    "    return list(parameters.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = best_param(model_precision)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------+--------------------+--------------------+\n",
      "|userID|               label|          prediction|\n",
      "+------+--------------------+--------------------+\n",
      "|   148|[50, 98, 114, 164...|[253, 50, 57, 144...|\n",
      "|   463|[20, 93, 112, 117...|[887, 408, 900, 1...|\n",
      "|   471|[8, 99, 588, 627,...|[1159, 982, 681, ...|\n",
      "|   496|[28, 38, 56, 64, ...|[1240, 1368, 89, ...|\n",
      "|   833|[5, 13, 23, 28, 3...|[320, 646, 305, 1...|\n",
      "+------+--------------------+--------------------+\n",
      "\n",
      "Run 0:\n",
      "\tValidation Metric: 0.03347502656748139\n",
      "\tParam(parent='ALSTopK_4e2995821d29445956ec', name='rank', doc='rank of the factorization'): 10\n",
      "\tParam(parent='ALSTopK_4e2995821d29445956ec', name='regParam', doc='regularization parameter (>= 0).'): 0.1\n",
      "Run 1:\n",
      "\tValidation Metric: 0.010095642933049955\n",
      "\tParam(parent='ALSTopK_4e2995821d29445956ec', name='rank', doc='rank of the factorization'): 10\n",
      "\tParam(parent='ALSTopK_4e2995821d29445956ec', name='regParam', doc='regularization parameter (>= 0).'): 0.01\n",
      "Run 2:\n",
      "\tValidation Metric: 0.04261424017003189\n",
      "\tParam(parent='ALSTopK_4e2995821d29445956ec', name='rank', doc='rank of the factorization'): 20\n",
      "\tParam(parent='ALSTopK_4e2995821d29445956ec', name='regParam', doc='regularization parameter (>= 0).'): 0.1\n",
      "Run 3:\n",
      "\tValidation Metric: 0.026461211477151955\n",
      "\tParam(parent='ALSTopK_4e2995821d29445956ec', name='rank', doc='rank of the factorization'): 20\n",
      "\tParam(parent='ALSTopK_4e2995821d29445956ec', name='regParam', doc='regularization parameter (>= 0).'): 0.01\n"
     ]
    }
   ],
   "source": [
    "model_precision.bestModel.transform(valid).limit(5).show()\n",
    "\n",
    "for idx, item in enumerate(model_precision.getEstimatorParamMaps()):\n",
    "    print('Run {}:'.format(idx))\n",
    "    print('\\tValidation Metric: {}'.format(model_precision.validationMetrics[idx]))\n",
    "    for key, value in item.items():\n",
    "        print('\\t{0}: {1}'.format(repr(key), value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Hyperparameter tuning with `hyperopt`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`hyperopt` is an open source Python package that is designed for tuning parameters for generic function with any pre-defined loss. More information about `hyperopt` can be found [here](https://github.com/hyperopt/hyperopt). `hyperopt` supports parallelization on MongoDB but not Spark. In our case, the tuning is performed in a sequential mode on a local computer.\n",
    "\n",
    "In `hyperopt`, an *objective* function is defined for optimizing the hyper parameters. In this case, the objective is similar to that in the Spark native construct situation, which is *to the RMSE metric for an ALS recommender*. Parameters of `rank` and `regParam` are used as hyperparameters. \n",
    "\n",
    "The objective function shown below demonstrates a RMSE loss for an ALS recommender. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Customize an objective function\n",
    "def objective(params):\n",
    "    time_run_start = time.time()\n",
    "    \n",
    "    rank = params['rank']\n",
    "    reg = params['reg']\n",
    "    train = params['train'] \n",
    "    valid = params['valid'] \n",
    "    col_user = params['col_user'] \n",
    "    col_item = params['col_item']\n",
    "    col_rating = params['col_rating'] \n",
    "    col_prediction = params['col_prediction'] \n",
    "    k = params['k']\n",
    "    relevancy_method = params['relevancy_method']\n",
    "    \n",
    "    als = ALS(\n",
    "        rank=rank,\n",
    "        maxIter=15,\n",
    "        implicitPrefs=False,\n",
    "        alpha=0.1,\n",
    "        regParam=reg,\n",
    "        coldStartStrategy='drop',\n",
    "        nonnegative=False,\n",
    "        seed=42,\n",
    "        **HEADER_ALS\n",
    "    )\n",
    "    \n",
    "    model = als.fit(train) \n",
    "    prediction = model.transform(valid)\n",
    "\n",
    "    rating_eval = SparkRatingEvaluation(\n",
    "        valid, \n",
    "        prediction, \n",
    "        **HEADER\n",
    "    )\n",
    "    \n",
    "    rmse = rating_eval.rmse()\n",
    "    \n",
    "    # Return the objective function result.\n",
    "    return {\n",
    "        'loss': rmse,\n",
    "        'status': STATUS_OK,\n",
    "        'eval_time': time.time() - time_run_start\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A search space is usually defined for hyperparameter exploration. Design of search space is empirical, and depends on the understanding of how distribution of parameter of interest affects the model performance measured by the loss function. \n",
    "\n",
    "In the ALS algorithm, the two hyper parameters, rank and reg, affect model performance in a way that\n",
    "* The higher the rank, the better the model performance but also the higher risk of overfitting.\n",
    "* The reg parameter prevents overfitting in certain way. \n",
    "\n",
    "Therefore, in this case, a uniform distribution and a lognormal distribution sampling spaces are used for rank and reg, respectively. A narrow search space is used for illustration purpose, that is, the range of rank is from 10 to 20, while that of reg is from $e^{-5}$ to $e^{-1}$. Together with the randomly sampled hyper parameters, other parameters use for building / evaluating the recommender, like `k`, column names, data, etc., are kept as constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a search space\n",
    "space = {\n",
    "    'rank': hp.quniform('rank', 10, 40, 5),\n",
    "    'reg': hp.loguniform('reg', -5, -1),\n",
    "    'train': train, \n",
    "    'valid': valid, \n",
    "    'col_user': COL_USER, \n",
    "    'col_item': COL_ITEM, \n",
    "    'col_rating': COL_RATING, \n",
    "    'col_prediction': \"prediction\", \n",
    "    'k': 10,\n",
    "    'relevancy_method': \"top_k\"\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Hyperparameter tuning with TPE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`fmin` of `hyperopt` is used for running the trials for searching optimal hyper parameters. In `hyperopt`, there are different strategies for intelligently optimize hyper parameters. For example, `hyperopt` avails [Tree of Parzen Estimators (TPE) method](https://papers.nips.cc/paper/4443-algorithms-for-hyper-parameter-optimization.pdf) for searching optimal parameters. \n",
    "\n",
    "The TPE method models a surface response of $p(x|y)$ by transforming a generative process, replacing the distributions of the configuration prior with non-parametric densities, where $p$ is the probability of configuration space $x$ given the loss $y$. For different configuration space, the TPE method does different replacements. That is, uniform $\\to$ truncated Gaussian mixture, log-uniform $\\to$ exponentiated truncated Gaussian mixture, categorical $\\to$ re-weighted categorical, etc. Using different observations ${x(1), ..., x(k)}$ in the non-parametric densities, these substitutions represent a learning algorithm that can produce a variety of densities over the configuration space $X$. By maintaining sorted lists of observed variables in $H$, the runtime of each iteration of the TPE algorithm can scale linearly in $|H|$ and linearly in the number of variables (dimensions) being optimized. In a nutshell, the algorithm recognizes the irrelevant variables in the configuration space, and thus reduces iterations in searching for the optimal ones. Details of the TPE algorithm can be found in the reference paper.\n",
    "\n",
    "The following runs the trials with the pre-defined objective function and search space. TPE is used as the optimization method. Totally there will be 10 evaluations run for searching the best parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_start = time.time()\n",
    "\n",
    "# Trials for recording each iteration of the hyperparameter searching.\n",
    "trials = Trials()\n",
    "\n",
    "best = fmin(\n",
    "    fn=objective,\n",
    "    space=space,\n",
    "    algo=tpe.suggest,\n",
    "    trials=trials,\n",
    "    max_evals=NUMBER_ITERATIONS\n",
    ")\n",
    "                  \n",
    "time_hyperopt = time.time() - time_start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'state': 2,\n",
       " 'tid': 22,\n",
       " 'spec': None,\n",
       " 'result': {'loss': 0.9371888478720254,\n",
       "  'status': 'ok',\n",
       "  'eval_time': 1551164708.9035273},\n",
       " 'misc': {'tid': 22,\n",
       "  'cmd': ('domain_attachment', 'FMinIter_Domain'),\n",
       "  'workdir': None,\n",
       "  'idxs': {'rank': [22], 'reg': [22]},\n",
       "  'vals': {'rank': [10.0], 'reg': [0.1529975102396072]}},\n",
       " 'exp_key': None,\n",
       " 'owner': None,\n",
       " 'version': 0,\n",
       " 'book_time': datetime.datetime(2019, 2, 26, 7, 4, 41, 559000),\n",
       " 'refresh_time': datetime.datetime(2019, 2, 26, 7, 5, 8, 903000)}"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trials.best_trial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAAE/CAYAAAAOkIE9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+UZGdd5/H3tzJJYAlIYLohZjKZxM1u6HV+BNvALrgGXUP4sQT8sSa2QwQ8OevCAY7LKqiHKJxdWcVd5YCGEdsBpIOKspvjiYRgglkXGzMzTLqT5teQ4E4PIT1JQA3kkFD93T/qzqSmma7unq7qW0/V+3VOna7n3qrq7607Xc98bj3PvZGZSJIkSZLK06i7AEmSJEnSqTHQSZIkSVKhDHSSJEmSVCgDnSRJkiQVykAnSZIkSYUy0EmSJElSoQx0Up+LiE9GxM/WXYckSZL6j4FOkiRJkgploJN6LCI21V2DJEn9zv5SOjUGOqkHIuLLEfGLETEDfCMifiUivhQR/xQRcxHxirbH/kxE/E1EvDMivhYR90bEi5Z53XMiYiYi/suGbYwkST1ykv5ya0T8WUQcrfrD17c99okR8f6qr/xsRPxCRMzXWL7UFwx0Uu9cDbwEeCrweeAHgO8Cfg34o4g4p+2xz6kesxn4DeAPIiLaXywiLgD+Gnh3Zv5m78uXJGlDHOsvnwZ8FLgTOBf4YeCNEfHC6nHXAduAC4EfAX56wyuV+pCBTuqdd2Xm4cx8JDP/NDO/kpmLmfnHwBeBS9se+/eZ+fuZ2QTeD5wDPKNt/RhwG3BdZu7ZsC2QJKn33pWZh4HvBUYy822Z+Whm3gP8PnBV9bj/APy3zPxaZs4D76qpXqmvOFZZ6p3Dx+5ExCuBn6d1ZBHgLFrfxh3z1WN3MvOb1ZdzZ7WtnwAOAR/pUa2SJNXlWH95PvDdEfH1tnWnAf+nuv/dbY9lyX1paPkNndQ7CRAR59M6wvg64OmZ+VTgLiA6PHepXwUeAKYi4rQu1ylJUp2y+nkYuDczn9p2e3Jmvrhafx+wpe15521olVKfMtBJvfckWp3VUYCIeBWtYSVr8RjwE9VrfSAi/NuVJA2avwP+qTpJyhMj4rSI+N6I+P5q/Z8Ab4mIsyPiXFoHSqWh538KpR7LzDngt4C/Be4HtgP/9xRe51HgR2nNrZs01EmSBkk1j/ylwC7gXlojU95H64RiAG8D5qt1n6A1DeFbG1+p1F8iM1d+lCRJktRHIuLngKsy8wfrrkWqk0f4JUmS1Peqa7E+LyIaEfEvgf9M6zIH0lDzLJeSJEkqwRnAe4ELgK8DHwZ+t9aKpD7gkEtJkiRJKpRDLiVJkiSpUAY6SZIkSSpUX86h27x5c27btq3uMiRJPbZ///4HMnOk7jpKYf8oScNjtX1kXwa6bdu2sW/fvrrLkCT1WET8fd01lMT+UZKGx2r7SIdcSpIkSVKhDHSSJEmSVCgDnSRJkiQVykAnSZIkSYUy0EmSJElSoQx0kiRJklQoA50kSZIkFcpAJ0mSJEmFMtBJkiRJUqEMdJIkSZJUKAOdJEmSJBXKQCdJkiRJhTLQSZIkSVKhDHSSJEmSVKhNdRcgafA1m4tMTc0yN3eUsbERJia202h4PEnqF02azDLDAguMMsoOdtLwmK8kFcFAJ6nnpqZmmZw8CMD09BEAdu/eWWdJktrMMsMB9gMwz2EAdnFJnSVJklZpVYffImIyIhYi4q5l1l8ZETMRcTAi9kXE89vWXRMRX6xu13SrcEnlmJs72rEtlWpQ+scFFjq2JUn9a7XjKfYCV3RY/1fAzszcBbwaeB9ARDwNuA54DnApcF1EnH3K1Uoq0tjYSMe2VLC9DED/OMpox7YkqX+tashlZt4eEds6rH+4rfkkIKv7LwRuycyHACLiFlod3w2nUqykMk1MbAc4YQ6dNAgGpX/cQWsIdPscOklSGbo2hy4iXgH8OjAKvKRafC5Ug/Fb5qtlkoZIo9FwzpyGVgn9Y4OGc+YkqVBdO4VVZn40My8GXg68fa3Pj4hrq/kF+44edX6NJGkw2D9Kknqp6+ckzszbgQsjYjNwBDivbfWWatnJnrcnM8czc3xkxPk1kqTBYv8oSeqFrgS6iPjnERHV/WcDZwIPAjcDl0fE2dVk78urZZIkDTz7R0lSr61qDl1E3ABcBmyOiHlaZ+Y6HSAzrwd+DHhlRDwGPAL8ZGYm8FBEvB24o3qptx2bAC5JUunsHyVJdYtWv9JfxsfHc9++fXWXIUnqsYjYn5njdddRCvtHSRoeq+0juz6HTpIkSZK0Mbp22QJtrGZzkamp2ROu69VomM8lSWvTpMksMydcg67h8V5JKoaBrlBTU7NMTh4EYHq6dWI0r/MlSVqrWWY4wH4A5qtL43lNOkkqh4fgCjU3d7RjW5Kk1VhgoWNbktTfDHSFGhsb6diWJGk1Rhnt2JYk9TeHXBZqYmI7wAlz6CRJWqsdtIbrt8+hkySVw0BXqEaj4Zw5SdK6NWg4Z06SCuaQS0mSJEkqlIFOkiRJkgploJMkSZKkQhnoJEmSJKlQBjpJkiRJKpSBTpIkSZIKZaCTJEmSpEIZ6CRJkiSpUAY6SZIkSSqUgU6SJEmSCmWgkyRJkqRCGegkSZIkqVAGOkmSJEkq1Ka6C5A0+JrNRaamZpmbO8rY2AgTE9tpNDyeJPWjJk1mmWGBBUYZZQc7aXj8V5L6loFOUs9NTc0yOXkQgOnpIwDs3r2zzpIkLWOWGQ6wH4B5DgOwi0vqLEmS1IGH3CT13Nzc0Y5tSf1jgYWObUlSf1kx0EXEZEQsRMRdy6yfiIiZiJiNiE9FxM62dV+ulh+MiH3dLFxSOcbGRjq2pVINYh85ymjHtiSpv6xmyOVe4N3AB5ZZfy/wg5n5tYh4EbAHeE7b+hdk5gPrqlJS0SYmtgOcMIdOGhB7GbA+cgetzNk+h06S1L9WDHSZeXtEbOuw/lNtzWlgy/rLkjRIGo2Gc+Y0kAaxj2zQcM6cJBWk23PoXgP8ZVs7gY9HxP6IuLbLv0uSpJLYR0qSuq5rZ7mMiBfQ6qye37b4+Zl5JCJGgVsi4nOZefsyz78WuBZg69at3SpLkqTaraePtH+UJHXSlW/oImIH8D7gysx88NjyzDxS/VwAPgpcutxrZOaezBzPzPGREU+YIEkaDOvtI+0fJUmdrDvQRcRW4M+B3Zn5hbblT4qIJx+7D1wOnPQsYJIkDSL7SElSr6045DIibgAuAzZHxDxwHXA6QGZeD7wVeDrwuxEB8O3MHAeeAXy0WrYJmMrMj/VgGyRJqoV9pCSpbqs5y+XVK6z/WeBnT7L8HvBcx5KkwWUfKUmqW7fPcilJkiRJ2iAGOkmSJEkqlIFOkiRJkgploJMkSZKkQhnoJEmSJKlQBjpJkiRJKpSBTpIkSZIKZaCTJEmSpEIZ6CRJkiSpUJvqLkBaTrO5yNTULHNzRxkbG2FiYjuNhscgJGmjNGkyywwLLDDKKDvYScNjwZLUVwx0hRqGsDM1Ncvk5EEApqePALB79846S5Kkgdce4h7hER7kAQDmOQzALi6pszxJ0hIGukINQ9iZmzvasS1J6r5ZZjjAfgAe4kHO4EzO4iwAFlioszRJ0kkM1lc6Q2QYws7Y2EjHtiSp+9pD2+mcwWM8erw9ymgdJUmSOvAbukKNjY0c/2buWHvQTExsBzhhWKkkqbdGGT0+vPJJPIkRRnkiTzw+h06S1F8MdIUahrDTaDQGbhipJPW7Y6HNE6FIUhkMdIUy7EiSeqFBwxOfSFJBPOQmSZIkSYUy0EmSJElSoQx0kiRJklQoA50kSZIkFcpAJ0mSJEmF8iyXhWo2F5mamj3hsgWNhvlcktQdTZrMMuPlCySpzxnoCjU1Ncvk5EGA4xcY9zIGkqRumWWGA+wHOH6hcS9nIEn9Z8VDbRExGRELEXHXMusnImImImYj4lMRsbNt3RUR8fmIOBQRb+5m4cNubu5ox7YkqfcGuY9cYKFjW5LUH1YzdmIvcEWH9fcCP5iZ24G3A3sAIuI04D3Ai4Ax4OqIGFtXtTpubGykY1uStCH2MqB95CijHduSpP6w4pDLzLw9IrZ1WP+ptuY0sKW6fylwKDPvAYiIDwNXAnOnWqweNzGxHeCEOXSSpI01yH3kDlpfJrbPoZMk9Z9uz6F7DfCX1f1zoRp03zIPPKfLv29oNRoN58xJUlmK6iMbNJwzJ0kF6Fqgi4gX0Oqsnn+Kz78WuBZg69at3SpLkqTaraePtH+UJHXSlfMPR8QO4H3AlZn5YLX4CHBe28O2VMtOKjP3ZOZ4Zo6PjDgfTJI0GNbbR9o/SpI6WXegi4itwJ8DuzPzC22r7gAuiogLIuIM4CrgxvX+PkmSSmEfKUnqtRWHXEbEDcBlwOaImAeuA04HyMzrgbcCTwd+NyIAvl0dSfx2RLwOuBk4DZjMzLt7shWSJNXAPlKSVLfIzLpr+A7j4+O5b9++usuQJPVYROzPzPG66yiF/aMkDY/V9pFdmUMnSZIkSdp4BjpJkiRJKpSBTpIkSZIKZaCTJEmSpEIZ6CRJkiSpUAY6SZIkSSqUgU6SJEmSCmWgkyRJkqRCGegkSZIkqVAGOkmSJEkqlIFOkiRJkgploJMkSZKkQhnoJEmSJKlQBjpJkiRJKpSBTpIkSZIKZaCTJEmSpEIZ6CRJkiSpUAY6SZIkSSqUgU6SJEmSCmWgkyRJkqRCGegkSZIkqVAGOkmSJEkq1Ka6C5CW02wuMjU1y9zcUcbGRpiY2E6j4TEISeoHTZrMMsMCC4wyyg520vA4sSRtuBU/eSNiMiIWIuKuZdZfHBF/GxHfiog3LVn35YiYjYiDEbGvW0VrOExNzTI5eZDp6SNMTh7kQx+arbskSTrBMPeRs8xwgP3Mc5gD7GeGO+suSZKG0moOpe0Fruiw/iHg9cA7l1n/gszclZnja6xNQ25u7mjHtiT1gb0MaR+5wELHtiRpY6wY6DLzdlod0nLrFzLzDuCxbhYmjY2NdGxLUt2GuY8cZbRjW5K0MXo9hy6Bj0dEAu/NzD09/n0aIBMT2wFOmEMnSQOk6D5yBzsBTphDJ0naeL0OdM/PzCMRMQrcEhGfq45mfoeIuBa4FmDr1q09LkslaDQa7N7tfxAkDaxV9ZH92j82aLCLS+ouQ5KGXk9PR5WZR6qfC8BHgUs7PHZPZo5n5vjIiEPrJEmDbbV9pP2jJKmTngW6iHhSRDz52H3gcuCkZwGTJGmY2EdKkrplxSGXEXEDcBmwOSLmgeuA0wEy8/qIeCawD3gKsBgRbwTGgM3ARyPi2O+ZysyP9WIjJEmqg32kJKluKwa6zLx6hfVfBbacZNU/gjOkJUmDyz5SklS3ns6hkyRJkiT1joFOkiRJkgploJMkSZKkQhnoJEmSJKlQBjpJkiRJKpSBTpIkSZIKZaCTJEmSpEIZ6CRJkiSpUAY6SZIkSSqUgU6SJEmSCmWgkyRJkqRCGegkSZIkqVAGOkmSJEkqlIFOkiRJkgploJMkSZKkQhnoJEmSJKlQBjpJkiRJKpSBTpIkSZIKZaCTJEmSpEIZ6CRJkiSpUAY6SZIkSSqUgU6SJEmSCmWgkyRJkqRCrRjoImIyIhYi4q5l1l8cEX8bEd+KiDctWXdFRHw+Ig5FxJu7VbQkSf3APlKSVLfVfEO3F7iiw/qHgNcD72xfGBGnAe8BXgSMAVdHxNiplSlJUl/ai32kJKlGKwa6zLydVoe03PqFzLwDeGzJqkuBQ5l5T2Y+CnwYuHI9xUqS1E/sIyVJdevlHLpzgcNt7flqmSRJw84+UpLUFX1zUpSIuDYi9kXEvqNHj9ZdjiRJfcH+UZLUSS8D3RHgvLb2lmrZSWXmnswcz8zxkZGRHpYlSVLtVt1H2j9KkjrZ1MPXvgO4KCIuoNVJXQX8VA9/HwDN5iJTU7PMzR1lbGyEiYntNBp980WkJElQUx/ZLU2azDLDAguMMsoOdtLon0E/kjRUVgx0EXEDcBmwOSLmgeuA0wEy8/qIeCawD3gKsBgRbwTGMvMfI+J1wM3AacBkZt7dm8143NTULJOTBwGYnm4d7Ny9e2evf60kaQiV1kd2yywzHGA/APPVVMBdXFJnSZI0tFYMdJl59Qrrv0prqMjJ1t0E3HRqpZ2aubmjHduSJHVLaX1ktyyw0LEtSdo4Azc+YmxspGNb6jfN5iIf/OCdvOUtn+CDH7yTxcXFukuSpI5GGe3YbtLkIJ/h49zMQT7DIn6uSVKv9HIOXS0mJrYDnDCHTupnDhOWVJodtD6j2ufQtXNIpiRtnIELdI1Gw/8MqygOE5ZUmgaNjgHNIZmStHEGbsilVBqHCUsaNCsNyZQkdc/AfUMnlcZhwpIGzUpDMiVJ3WOgk2rmMGFJg2alIZmSpO5xyKUkSZIkFcpAJ0mSJEmFMtBJkiRJUqEMdJIkSZJUKAOdJEmSJBXKQCdJkiRJhTLQSZIkSVKhDHSSJEmSVCgDnSRJkiQVykAnSZIkSYUy0EmSJElSoQx0kiRJklQoA50kSZIkFcpAJ0mSJEmFMtBJkiRJUqEMdJIkSZJUKAOdJEmSJBVqxUAXEZMRsRARdy2zPiLiXRFxKCJmIuLZbeuaEXGwut3YzcIlSaqbfaQkqW6r+YZuL3BFh/UvAi6qbtcCv9e27pHM3FXdXnbKVUqS1J/2Yh8pSarRioEuM28HHurwkCuBD2TLNPDUiDinWwVKktSv7CMlSXXrxhy6c4HDbe35ahnAEyJiX0RMR8TLu/C7JEkqiX2kJKmnNvX49c/PzCMRcSFwa0TMZuaXTvbAiLiW1nAUtm7d2uOyJEmq3ar6SPtHSVIn3fiG7ghwXlt7S7WMzDz28x7gk8Aly71IZu7JzPHMHB8ZGelCWZIk1W7dfaT9oySpk24EuhuBV1Zn8nou8A+ZeV9EnB0RZwJExGbgecBcF36fJEmlsI+UJPXUikMuI+IG4DJgc0TMA9cBpwNk5vXATcCLgUPAN4FXVU99FvDeiFikFRzfkZl2VpKkgWEfKUmq24qBLjOvXmF9Aq89yfJPAdtPvTRJkvqbfaQkqW7dGHIpSZIkSaqBgU6SJEmSCmWgkyRJkqRCGegkSZIkqVC9vrC4JNFsLjI1Ncvc3FHGxkaYmNhOo+HxJGmYNWkyywwLLDDKKDvYScPjzJK0ZgY6ST03NTXL5ORBAKanjwCwe/fOOkuSVLNZZjjAfgDmOQzArpNfW12S1IGHwiT13Nzc0Y5tScNngYWObUnS6hjoJPXc2NhIx7ak4TPKaMe2JGl1HHIpqecmJlrXT26fQydpuO2gNey6fQ6dJGntDHSSeq7RaDhnTtIJGjScMydJXeCQS0mSJEkqlIFOkiRJkgploJMkSZKkQhnoJEmSJKlQBjpJkiRJKpSBTpIkSZIKZaCTJEmSpEIZ6CRJkiSpUAY6SZIkSSrUproLkKRB0WwuMjU1y9zcUcbGRpiY2E6j4XEzSZIGTZMms8ywwAKjjLKDnTRq+q7MQCdJXTI1Ncvk5EEApqePALB79846S5IkST0wywwH2A/APIcB2MUltdTioWNJ6pK5uaMd25IkaTAssNCxvZEMdJLUJWNjIx3bkiRpMIwy2rG9kVYV6CJiMiIWIuKuZdZHRLwrIg5FxExEPLtt3TUR8cXqdk23CpekfjMxsZ1Xv3oXz33uubz61buYmNhed0nqMftHSRpOO9jJs/k+tnAez+b72EF9UyxWO4duL/Bu4APLrH8RcFF1ew7we8BzIuJpwHXAOJDA/oi4MTO/tp6iJakfNRoN58wNn73YP0rS0GnQqG3O3FKr+oYuM28HHurwkCuBD2TLNPDUiDgHeCFwS2Y+VHVStwBXrLdoSZL6gf2jJKlu3ZpDdy5Up3dpma+WLbdckqRhYP8oSeqpvjkpSkRcGxH7ImLf0aOeGU6SJLB/lCR11q1AdwQ4r629pVq23PLvkJl7MnM8M8dHRjwznCRpINg/SpJ6qluB7kbgldXZvJ4L/ENm3gfcDFweEWdHxNnA5dUySZKGgf2jJKmnVnWWy4i4AbgM2BwR87TOzHU6QGZeD9wEvBg4BHwTeFW17qGIeDtwR/VSb8vMTpPHJUkqhv2jJKluqwp0mXn1CusTeO0y6yaBybWXpk6azUWmpmaZmzvK2NgIExPbaTT6Zkqk1sB9KZXL/rF/NGkyywwLLDDKKDvYSaN/ThUg6RT5t72y1V6HTn1mamqWycmDAExPt6ZdeP2rMrkvJWn9ZpnhAPsBmK9OINov14iSdOr8216Z8bZQc3NHO7ZVDvelJK3fAgsd25LK5N/2ygx0hRobG+nYVjkuvngz99//MPfc8xD33/8wF1+8ue6SJKk4o4wev7/IIo/wCB/nZg7yGRZZrLEySevR/rd9srYcclmsiYntACfMu1KpkswkE1rTbbLugiSpODtoDVVfYIFHeIQHeQBwiJZUuva/7WNz6HQiA12hGo2G86wGxOc+9yDPfOaTT2hLktamQeN4aPv4kitAOERLKlf737ZOziGXUs0cPitJ3eUQLUnDxG/opJo5fFaSusshWpKGiYFOqpnDZyWpuxyiJWmYOORSkiRJkgploJMkSZKkQjnkUpK0Js3mIlNTsyfM+2w0PD4oSSVq0mSWmRPmnDb8zqcoBjpJ0ppMTc0yOXkQgOnpIwDOA5WkQs0ywwH2A163sVTGb0nSmszNHe3YliSVY+l1Gr1uY3kMdJKkNfHaiZI0OLxuY/kccilJWhOvnShJg8PrNpbPQCdJWhOvnShJg8PrNpbPIZeSJEmSVCgDnSRJkiQVykAnSZIkSYUy0EmSJElSoQx0kiRJklQoA50kSZIkFcrLFkg1azYXmZqaPeGaXo2Gx1okqZ81aTLLzAnX7mp4nFxSDVYV6CLiCuB3gNOA92XmO5asPx+YBEaAh4Cfzsz5al0TmK0e+v8y82Vdql0aCFNTs0xOHgRgevoIgNf4kgph/zi8ZpnhAPsBmOcwgNfyklSLFQNdRJwGvAf4EWAeuCMibszMubaHvRP4QGa+PyJ+CPh1YHe17pHM3NXluqWBMTd3tGNbUn+yfxxuCyx0bEvSRlnN2IBLgUOZeU9mPgp8GLhyyWPGgFur+7edZL2kZYyNjXRsS+pb9o9DbJTRjm1J2iirGXJ5LlRjCVrmgecsecydwI/SGnbyCuDJEfH0zHwQeEJE7AO+DbwjM//X+svWMBiWuWUTE9sBTthOSUWwfxxiO2gNjW+fQ6cyOR9SpevWSVHeBLw7In4GuB04AjSrdedn5pGIuBC4NSJmM/NLS18gIq4FrgXYunVrl8oaXMMQdoZlblmj0RjI7ZIE2D8OrAaNU5ozZ3joP86HVOlWE+iOAOe1tbdUy47LzK/QOgJJRJwF/Fhmfr1ad6T6eU9EfBK4BPiODisz9wB7AMbHx3OtG3LMMAQdGI6w49wySX2uqP5R/eFYeFhkkc8yxywzbGeHwa5GzodU6VbzyXEHcFFEXBARZwBXATe2PyAiNkfEsdd6C60zehERZ0fEmcceAzwPaJ8s3nXHgs709BEmJw/yoQ/NrvykAg1D2HFumaQ+V1T/qP5wLCx8g2/wDR5mgfs5wH5muLPmyoaX8yFVuhW/ocvMb0fE64CbaZ2WeTIz746ItwH7MvNG4DLg1yMiaQ0peW319GcB742IRVrh8R1Lzv7VdcMQdKAVbo59M3esPWicWyapn5XWP6o/jDLKPId5jEcBOJ0zAL8VqpPzIVW6Vc2hy8ybgJuWLHtr2/2PAB85yfM+BWzo/8KHIejAcISdYZlbNizDhKVBVFL/qP5wLCzMMsPX+BpncRbQ+Vsh59311qnOh5T6RbdOitI3fuInnsUf/uFn+NznHuDiizfzkz/5rLpL6olhCDvDEnSGYT6kJKnlWHjYwU5muHNV3wpt9Ek7ehkgH+NRbuImvsp9PJNzeAkvZdPg/XdU2lAD9xf0+td/jE9/+is0m4t8+tNf4XWv+xh79rys7rK6bhjCzrAEnWEZJixJetxavhVa60k71hvIlgbIRRZp0OhKwLuJm7ib1vkNHuQBAK7k5af0WpJaBi7Q3XrrvTz6aOuM0M1mk1tvvbfminpjGMLOsASdYRkmLEk6Ncfm3bW3O1nvN3pLA+Pd3MVjPHbKr9fuq9zXsS1p7QYu0J111hkd24NiGMLOsASdYZgPKUk6dWs9acd6T8O/NECu9Ppr8UzOOf7N3LG2pPUZuED3hjdcyi/90m08/PCjnHXWGbzhDZfWXVJPDEPYGZagMwzzISVJp26tJ+1Y6zd6Sy0NkEnyGQ6c8uu1ewkvBThhDp2k9Rm4QHfNNZewadOmgQ8BwxB2DDqSJK3dek/DvzRALrJIEF05rf8mNjlnTuqygQt0wxIChmU7JUnS2nT7NPye1l/qb4N1WkRJkiRJGiIGOkmSJEkqlIFOkiRJkgploJMkSZKkQhnoJEmSJKlQBjpJkiRJKpSBTpIkSZIKZaCTJEmSpEIZ6CRJkiSpUAY6SZIkSSqUgU6SJEmSCmWgkyRJkqRCGegkSZIkqVAGOkmSJEkqVGRm3TV8h4g4Cvz9Ol9mM/BAF8rpd8OwncOwjTAc2zkM2whu51qcn5kj3ShmGKyzfyz532XJtUPZ9ZdcO5Rdf8m1Q9n190vtq+oj+zLQdUNE7MvM8brr6LVh2M5h2EYYju0chm0Et1P9qeT9VXLtUHb9JdcOZddfcu1Qdv2l1e6QS0mSJEkqlIFOkiRJkgo1yIFuT90FbJBh2M5h2EYYju0chm0Et1P9qeT9VXLtUHb9JdcOZddfcu1Qdv1F1T6wc+gkSZIkadAN8jd0kiRJkjTQBiLQRcRkRCxExF1ty54WEbdExBern2fXWeN6LbONvxoRRyLiYHV7cZ01dkNEnBcRt0XEXETcHRFvqJYPzP7ssI0DtT8j4gkR8XcRcWe1nb9WLb8gIj4dEYci4o8j4oy6az1VHbZxb0Tc27Yvd9VdazdExGkR8ZmI+IuqPTD7snQRcUVEfL7aF28+yfozq310qNpn29rWvaVa/vmIeOFG1l39/lOqPSK2RcQjbX9n1/dh7f82Ig5ExLcj4seXrLsc846TAAAGUUlEQVSm6tO+GBHXbFzVJ9Swnvqbbe/9jRtX9fHfv1LtP1/1szMR8VcRcX7buhLe+0719/t7/x8jYraq728iYqxtXa2fN1UNp1R/P3zmLCszi78B/xZ4NnBX27LfAN5c3X8z8N/rrrMH2/irwJvqrq3L23kO8Ozq/pOBLwBjg7Q/O2zjQO1PIICzqvunA58Gngv8CXBVtfx64OfqrrUH27gX+PG66+vB9v48MAX8RdUemH1Z8g04DfgScCFwBnAnMLbkMf8JuL66fxXwx9X9serxZwIXVK9zWiG1b2vvE/v0fd8G7AA+0P6ZADwNuKf6eXZ1/+xS6q/WPdzn7/0LgH9W3f+5tn83pbz3J62/kPf+KW33XwZ8rLpf6+dNF+qv9TOn020gvqHLzNuBh5YsvhJ4f3X//cDLN7SoLltmGwdOZt6XmQeq+/8EfBY4lwHanx22caBky8NV8/TqlsAPAR+plpe+L5fbxoETEVuAlwDvq9rBAO3Lwl0KHMrMezLzUeDDtD4z27V/hn4E+OFqH14JfDgzv5WZ9wKHqtfbKOupvW4r1p6ZX87MGWBxyXNfCNySmQ9l5teAW4ArNqLoNuupv26rqf22zPxm1ZwGtlT3S3nvl6u/bqup/R/bmk/i8X6x7s8bWF/9fWsgAt0ynpGZ91X3vwo8o85ieuh11dfxk1HwMMSTqYbVXELrW4+B3J9LthEGbH9WQ/QOAgu0Os0vAV/PzG9XD5mn8DC7dBsz89i+/K/VvvyfEXFmjSV2y28Dv8Dj/7F7OgO2Lwt2LnC4rX2yfXH8MdU++wda+3A1z+2l9dQOcEG0hgH/dUT8QK+LXa6uylreu7rf927U8ISI2BcR0xGx0Qdz1lr7a4C/PMXn9sJ66ocC3vuIeG1EfInWCKvXr+W5Pbae+qHez5xlDXKgOy5b35P2fbo+Bb8HfA+wC7gP+K16y+meiDgL+DPgjUuOlAzM/jzJNg7c/szMZmbuonVk8VLg4ppL6rql2xgR3wu8hda2fj+tYT2/WGOJ6xYRLwUWMnN/3bVIbe4DtmbmJVTDgSPiKTXXNEzOz8xx4KeA346I76m7oJOJiJ8GxoHfrLuWU7FM/X3/3mfmezLze2j1f79Sdz1rtUz9ffuZM8iB7v6IOAeg+rlQcz1dl5n3V/+ZXAR+n43/2ronIuJ0WkHnQ5n559XigdqfJ9vGQd2fAJn5deA24F8DT42ITdWqLcCR2grrorZtvKIaVpuZ+S3gDyl/Xz4PeFlEfJnW8JQfAn6HAd2XBToCnNfWPtm+OP6Yap99F/DgKp/bS6dcezVs60GA6mDDl4B/0fOKT1JXZS3vXd3v+7pryMwj1c97gE/SGm2yUVZVe0T8O+CXgZdVn8erfm6Praf+It77Nh/m8eH4xbz3bY7X3wefOcsa5EB3I3DszEXXAP+7xlp64ljAqbwCuGu5x5aimhfxB8BnM/N/tK0amP253DYO2v6MiJGIeGp1/4nAj9CaL3gbcOxsaaXvy5Nt4+faDj4ErY6g6H2ZmW/JzC2ZuY3WSSluzcwJBmhfFu4O4KJonXX0DFr7aOmZ79o/Q3+c1j7MavlV0TqT5AXARcDfbVDdsI7aq7+/0wAi4kJatd+zQXXD6mpfzs3A5RFxdjW8/vJq2UY65fqrus+s7m+mddBnrmeVfqcVa4+IS4D30gpD7QeBi3jvl6u/kPf+orbmS4AvVvfr/ryBddTfB585y1vubCkl3YAbaH0N+hitsbCvoTW+/q9o7YRPAE+ru84ebOMHgVlghtY/xnPqrrML2/l8WsMpZ4CD1e3Fg7Q/O2zjQO1PWmdG+0y1PXcBb62WX0jrA/wQ8KfAmXXX2oNtvLXal3cBf0R1JsxBuAGX8fhZLgdmX5Z+qz5DvkDriPEvV8veRus/gwBPqPbRoWqfXdj23F+unvd54EWl1A78GHB39Rl6APj3fVj791d99jdofSN6d9tzX11t0yHgVX367+ak9QP/pvqMu7P6+Zo+rP0TwP083s/eWNh7f9L6C3nvf6ftb/M24F+1PbfWz5v11N8PnznL3aIqUJIkSZJUmEEecilJkiRJA81AJ0mSJEmFMtBJkiRJUqEMdJIkSZJUKAOdJEmSJBXKQCdJkiRJhTLQSZIkSVKhDHSSJEmSVKj/D0lpTCR870LjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "parameters = ['rank', 'reg']\n",
    "cols = len(parameters)\n",
    "f, axes = plt.subplots(nrows=1, ncols=cols, figsize=(15,5))\n",
    "cmap = plt.cm.jet\n",
    "for i, val in enumerate(parameters):\n",
    "    xs = np.array([t['misc']['vals'][val] for t in trials.trials]).ravel()\n",
    "    ys = [t['result']['loss'] for t in trials.trials]\n",
    "    xs, ys = zip(*sorted(zip(xs, ys)))\n",
    "    ys = np.array(ys)\n",
    "    axes[i].scatter(xs, ys, s=20, linewidth=0.01, alpha=0.75, c=cmap(float(i)/len(parameters)))\n",
    "    axes[i].set_title(val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be seen from the above plot that\n",
    "* The actual impact of rank is in line with the intuition - the smaller the value the better the result.\n",
    "* It is interesting to see that the optimal value of reg is around 0.1 to 0.15. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the best model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "als = ALS(\n",
    "    rank=best[\"rank\"],\n",
    "    regParam=best[\"reg\"],\n",
    "    maxIter=15,\n",
    "    implicitPrefs=False,\n",
    "    alpha=0.1,\n",
    "    coldStartStrategy='drop',\n",
    "    nonnegative=False,\n",
    "    seed=42,\n",
    "    **HEADER_ALS\n",
    ")\n",
    "    \n",
    "model_best_hyperopt = als.fit(train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tuning prameters against other metrics can be simply done by modifying the `objective` function. The following shows an objective function of how to tune \"precision@k\". Since `fmin` in `hyperopt` only supports minimization while the actual objective of the loss is to maximize \"precision@k\", `-precision` instead of `precision` is used in the returned value of the `objective` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Customize an objective function\n",
    "def objective_precision(params):\n",
    "    time_run_start = time.time()\n",
    "    \n",
    "    rank = params['rank']\n",
    "    reg = params['reg']\n",
    "    train = params['train'] \n",
    "    valid = params['valid'] \n",
    "    col_user = params['col_user'] \n",
    "    col_item = params['col_item']\n",
    "    col_rating = params['col_rating'] \n",
    "    col_prediction = params['col_prediction'] \n",
    "    k = params['k']\n",
    "    relevancy_method = params['relevancy_method']\n",
    "    \n",
    "    header = {\n",
    "        \"userCol\": col_user,\n",
    "        \"itemCol\": col_item,\n",
    "        \"ratingCol\": col_rating,\n",
    "    }\n",
    "    \n",
    "    als = ALS(\n",
    "        rank=rank,\n",
    "        maxIter=15,\n",
    "        implicitPrefs=False,\n",
    "        alpha=0.1,\n",
    "        regParam=reg,\n",
    "        coldStartStrategy='drop',\n",
    "        nonnegative=False,\n",
    "        seed=42,\n",
    "        **header\n",
    "    )\n",
    "    \n",
    "    model = als.fit(train)\n",
    "    \n",
    "    users = train.select(col_user).distinct()\n",
    "    items = train.select(col_item).distinct()\n",
    "    user_item = users.crossJoin(items)\n",
    "    dfs_pred = model.transform(user_item)\n",
    "\n",
    "    # Remove seen items.\n",
    "    dfs_pred_exclude_train = dfs_pred.alias(\"pred\").join(\n",
    "        train.alias(\"train\"),\n",
    "        (dfs_pred[col_user] == train[col_user]) & (dfs_pred[col_item] == train[col_item]),\n",
    "        how='outer'\n",
    "    )\n",
    "\n",
    "    top_all = dfs_pred_exclude_train.filter(dfs_pred_exclude_train[\"train.Rating\"].isNull()) \\\n",
    "        .select('pred.' + col_user, 'pred.' + col_item, 'pred.' + \"prediction\")\n",
    "    \n",
    "    top_all.cache().count()\n",
    "\n",
    "    rank_eval = SparkRankingEvaluation(\n",
    "        valid, \n",
    "        top_all, \n",
    "        k=k, \n",
    "        col_user=col_user, \n",
    "        col_item=col_item, \n",
    "        col_rating=\"rating\", \n",
    "        col_prediction=\"prediction\", \n",
    "        relevancy_method=relevancy_method\n",
    "    )\n",
    "    \n",
    "    precision = rank_eval.precision_at_k()\n",
    "    \n",
    "    # Return the objective function result.\n",
    "    return {\n",
    "        'loss': -precision,\n",
    "        'status': STATUS_OK,\n",
    "        'eval_time': time.time() - time_run_start\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Hyperparameter tuning with `hyperopt` sampling methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though `hyperopt` works well in a single node machine, its features (e.g., `Trials` module) do not support Spark environment, which makes it hard to perform the tuning tasks in a distributed/parallel manner. It is useful to use `hyperopt` for sampling parameter values from the pre-defined sampling space, and then parallelize the model training onto Spark cluster with the sampled parameter combinations.\n",
    "\n",
    "The downside of this method is that the intelligent searching algorithm (i.e., TPE) of `hyperopt` cannot be used. The approach introduced here is therefore equivalent to random search."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample the parameters used for model building from the pre-defined space. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_start = time.time()\n",
    "\n",
    "sample_params = [sample(space) for x in range(NUMBER_ITERATIONS)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following runs model building on the sampled parameter values with the pre-defined objective function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_map = list(map(lambda x: objective(x), sample_params))\n",
    "\n",
    "time_sample = time.time() - time_start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[{'loss': 0.9570030109029894, 'status': 'ok', 'eval_time': 26.525237321853638},\n",
       " {'loss': 1.3105647924197812, 'status': 'ok', 'eval_time': 27.681031227111816},\n",
       " {'loss': 0.9924061855935775, 'status': 'ok', 'eval_time': 27.931140184402466},\n",
       " {'loss': 0.9726441114619616, 'status': 'ok', 'eval_time': 27.279237031936646},\n",
       " {'loss': 1.2021897405557644, 'status': 'ok', 'eval_time': 27.261536836624146},\n",
       " {'loss': 1.3204032268351846, 'status': 'ok', 'eval_time': 26.627870321273804},\n",
       " {'loss': 1.197655798181204, 'status': 'ok', 'eval_time': 26.73567533493042},\n",
       " {'loss': 1.0223163707199912, 'status': 'ok', 'eval_time': 27.54908800125122},\n",
       " {'loss': 1.3050029526943239, 'status': 'ok', 'eval_time': 28.92825984954834},\n",
       " {'loss': 1.1397089932872841, 'status': 'ok', 'eval_time': 27.19080638885498},\n",
       " {'loss': 0.9652317477756671, 'status': 'ok', 'eval_time': 26.762747287750244},\n",
       " {'loss': 1.0554259497351064, 'status': 'ok', 'eval_time': 28.918659210205078},\n",
       " {'loss': 1.2077389711046154, 'status': 'ok', 'eval_time': 28.164798498153687},\n",
       " {'loss': 0.9388730625639623, 'status': 'ok', 'eval_time': 28.542875289916992},\n",
       " {'loss': 0.9544511060426731, 'status': 'ok', 'eval_time': 26.06899905204773},\n",
       " {'loss': 0.9487655863729143, 'status': 'ok', 'eval_time': 27.023435831069946},\n",
       " {'loss': 1.2799219029733113, 'status': 'ok', 'eval_time': 27.958611965179443},\n",
       " {'loss': 1.303827696884951, 'status': 'ok', 'eval_time': 28.117038011550903},\n",
       " {'loss': 1.3544475034466399, 'status': 'ok', 'eval_time': 28.15032386779785},\n",
       " {'loss': 0.9856455623489918, 'status': 'ok', 'eval_time': 26.025404930114746},\n",
       " {'loss': 1.2317934060803921, 'status': 'ok', 'eval_time': 26.83657145500183},\n",
       " {'loss': 1.3232453080512732, 'status': 'ok', 'eval_time': 27.046321868896484},\n",
       " {'loss': 0.9381256837211831, 'status': 'ok', 'eval_time': 28.209535837173462},\n",
       " {'loss': 1.01034858445376, 'status': 'ok', 'eval_time': 27.161163568496704},\n",
       " {'loss': 0.948607396466967, 'status': 'ok', 'eval_time': 26.905335187911987}]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the best model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_metrics = np.array([x['loss'] for x in results_map])\n",
    "best_loss = np.where(loss_metrics == min(loss_metrics))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_param = sample_params[best_loss[0].item()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "als = ALS(\n",
    "    rank=best_param[\"rank\"],\n",
    "    regParam=best_param[\"reg\"],\n",
    "    maxIter=15,\n",
    "    implicitPrefs=False,\n",
    "    alpha=0.1,\n",
    "    coldStartStrategy='drop',\n",
    "    nonnegative=False,\n",
    "    seed=42,\n",
    "    **HEADER_ALS\n",
    ")\n",
    "    \n",
    "model_best_sample = als.fit(train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 Evaluation on testing data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimal parameters can then be used for building a recommender, which is then evaluated on the testing data.\n",
    "\n",
    "The following codes generate the evaluation results by using the testing dataset with the optimal model selected against the pre-defined loss. Without loss of generity, in this case, the optimal model that performs the best w.r.t regression loss (i.e., the RMSE metric) is used. One can simply use other metrics like precision@k, as illustrated in the above sections, to evaluate the optimal model on the testing dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get prediction results with the optimal modesl from different approaches.\n",
    "prediction_spark = model_best_spark.transform(test)\n",
    "prediction_hyperopt = model_best_hyperopt.transform(test)\n",
    "prediction_sample = model_best_sample.transform(test)\n",
    "\n",
    "predictions = [prediction_spark, prediction_hyperopt, prediction_sample]\n",
    "elapsed = [time_spark, time_hyperopt, time_sample]\n",
    "\n",
    "approaches = ['spark', 'hyperopt', 'sample']\n",
    "test_evaluations = pd.DataFrame()\n",
    "for ind, approach in enumerate(approaches):    \n",
    "    rating_eval = SparkRatingEvaluation(\n",
    "        test, \n",
    "        predictions[ind],\n",
    "        **HEADER\n",
    "    )\n",
    "    \n",
    "    result = pd.DataFrame({\n",
    "        'Approach': approach,\n",
    "        'RMSE': rating_eval.rmse(),\n",
    "        'MAE': rating_eval.mae(),\n",
    "        'Explained variance': rating_eval.exp_var(),\n",
    "        'R squared': rating_eval.rsquared(),\n",
    "        'Elapsed': elapsed[ind]\n",
    "    }, index=[0])\n",
    "    \n",
    "    test_evaluations = test_evaluations.append(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Approach</th>\n",
       "      <th>RMSE</th>\n",
       "      <th>MAE</th>\n",
       "      <th>Explained variance</th>\n",
       "      <th>R squared</th>\n",
       "      <th>Elapsed</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>spark</td>\n",
       "      <td>0.917751</td>\n",
       "      <td>0.729715</td>\n",
       "      <td>0.346218</td>\n",
       "      <td>0.325242</td>\n",
       "      <td>286.972354</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>hyperopt</td>\n",
       "      <td>0.937614</td>\n",
       "      <td>0.748037</td>\n",
       "      <td>0.326001</td>\n",
       "      <td>0.295247</td>\n",
       "      <td>698.869952</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>sample</td>\n",
       "      <td>0.938145</td>\n",
       "      <td>0.749925</td>\n",
       "      <td>0.328112</td>\n",
       "      <td>0.294449</td>\n",
       "      <td>714.806368</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Approach      RMSE       MAE  Explained variance  R squared     Elapsed\n",
       "0     spark  0.917751  0.729715            0.346218   0.325242  286.972354\n",
       "0  hyperopt  0.937614  0.748037            0.326001   0.295247  698.869952\n",
       "0    sample  0.938145  0.749925            0.328112   0.294449  714.806368"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_evaluations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the results, it can be seen that, *with the same number of iterations*, Spark native construct based approach takes the least amount of time, even if there is no parallel computing. This is simply because Spark native constructs leverage the underlying Java codes for running the actual analytics with high performance efficiency. Interestingly, the run time for `hyperopt` with TPE algorithm and random search methods are almost the same. Possible reasons for this are that, the TPE algorithm searches optimal parameters intelligently but runs the tuning iterations sequentially. Also, the advantage of TPE may become obvious when there is a higher dimensionality of hyperparameters. \n",
    "\n",
    "The three approaches use the same RMSE loss. In this measure, the native Spark construct performs the best. The `hyperopt` based approach performs the second best, but the advantage is very subtle. It should be noted that these differences may be owing to many factors like characteristics of datasets, dimensionality of hyperparameter space, sampling size in the searching, etc. Note the differences in the RMSE metrics may also come from the randomness of the intermediate steps in parameter tuning process. In practice, multiple runs are required for generating statistically robust comparison results. We have tried 5 times for running the same comparison codes above. The results aligned well with each other in terms of objective metric values and elapsed time. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In summary, there are mainly three different approaches for running hyperparameter tuning for Spark based recommendation algorithm. The three different approaches are compared as follows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|Approach|Distributed (on Spark)|Param sampling|Advanced hyperparam searching algo|Custom evaluation metrics|Custom data split|\n",
    "|---------|-------------|--------------|--------------------------|--------------|------------|\n",
    "|AzureML Services|Parallelizing Spark sessions on multi-node cluster or single Spark session on one VM node.)|Random, Grid, Bayesian sampling for discrete and continuous variables.|Bandit policy, Median stopping policy, and truncation selection policy.|Yes|Yes|\n",
    "|Spark native construct|Distributed in single-node standalone Spark environment or multi-node Spark cluster.|No|No|Need to re-engineer Spark modules|Need to re-engineer Spark modules.|\n",
    "|`hyperopt`|No (only support parallelization on MongoDB)|Random sampling for discrete and continuous variables.|Tree Parzen Estimator|Yes|Yes|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "* Azure Machine Learning Services, url: https://azure.microsoft.com/en-us/services/machine-learning-service/\n",
    "* Lisa Li, *et al*, Hyperband: A Novel Bandit-Based Approach to Hyperparameter Optimization, The Journal of Machine Learning Research, Volume 18 Issue 1, pp 6765-6816, January 2017.\n",
    "* James Bergstrat *et al*, Algorithms for Hyper-Parameter Optimization, Procs 25th NIPS 2011. \n",
    "* `hyperopt`, url: http://hyperopt.github.io/hyperopt/.\n",
    "* Bergstra, J., Yamins, D., Cox, D. D. (2013) Making a Science of Model Search: Hyperparameter Optimization in Hundreds of Dimensions for Vision Architectures. Proc. of the 30th International Conference on Machine Learning (ICML 2013).\n",
    "* Kris Wright, \"Hyper parameter tuning with hyperopt\", url:https://districtdatalabs.silvrback.com/parameter-tuning-with-hyperopt"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3.6 (recommender)",
   "language": "python",
   "name": "recommender"
  },
  "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}