diff --git a/pricing_competing_products/price_optimization_with_competing_products.ipynb b/pricing_competing_products/price_optimization_with_competing_products.ipynb index d09d7ab..3b593aa 100644 --- a/pricing_competing_products/price_optimization_with_competing_products.ipynb +++ b/pricing_competing_products/price_optimization_with_competing_products.ipynb @@ -1,1725 +1,1293 @@ { - "cells": [ - { - "cell_type": "markdown", - "id": "7f7d69f9", - "metadata": { - "id": "7f7d69f9" - }, - "source": [ - "Copyright © 2024 Gurobi Optimization, LLC\n", - "\n", - "# Price Optimization with Competing Products\n", - "Finding the delicate balance between price and demand is a difficult problem to solve and one that is prevalent in a number of huge industries like retail, e-commerce, ticketing, and hospitality.\n", - "\n", - "It just so happens that this is a problem that data scientists have recently been getting better and better at addressing. Still, a key piece is missing: How can I make the right pricing decision given all the other constraints and business rules that exist for this problem? That's where optimization comes in. \n", - "\n", - "In this scenario, we have a product that has several categories, but we have a limited amount of \"space\" for our products. This could be shelf or warehouse space in retail, seats for events, or for airline ticketing, or rooms at a hotel. \n", - "\n", - "This problem considers two similar products that are offered where it's our job to use the data available with some optimization know-how to determine the optimal mix of products to offer that maximize revenue while also adhering to a few other business rules. First, we'll create a predictive model to forecast sales based on the prices of each product. Then we'll build an optimization model to find this optimal mix. Finally, we’ll also use the Gurobi-sponsored open-source package Gurobi Machine Learning to seamlessly combine the features of a machine learning model with a decision of an optimization model. \n", - "\n", - "Let’s get started!" - ] - }, - { - "cell_type": "markdown", - "id": "a5846a14", - "metadata": { - "id": "vGL8Rv-uRZ7j" - }, - "source": [ - "## Load required packages\n", - "If you have a Gurobi license you can skip the installation of `gurobipy`, but always make sure you have the [latest version](https://www.gurobi.com/downloads/gurobi-software) available. " - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "id": "1c3276f0", - "metadata": {}, - "outputs": [ + "cells": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Requirement already satisfied: gurobipy in /Users/yurchisin/opt/anaconda3/envs/gurobi_ml/lib/python3.11/site-packages (11.0.0)\n", - "Note: you may need to restart the kernel to use updated packages.\n" - ] - } - ], - "source": [ - "%pip install gurobipy" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "3024295a", - "metadata": { - "id": "3024295a" - }, - "outputs": [], - "source": [ - "import gurobipy as gp\n", - "from gurobipy import GRB\n", - "\n", - "import numpy as np\n", - "import pandas as pd\n", - "import seaborn as sns\n", - "import matplotlib.pyplot as plt\n", - "import warnings\n", - "from sklearn.model_selection import train_test_split\n", - "from sklearn import tree\n", - "\n", - "warnings.filterwarnings(\"ignore\")" - ] - }, - { - "cell_type": "markdown", - "id": "d8c28cfd", - "metadata": { - "id": "d8c28cfd" - }, - "source": [ - "## Start with some data analysis\n", - "\n", - "This data contains prices and sales for two of our competing products and was generated using another script, which can be found [here](). Let's load the data and take a quick look. " - ] - }, - { - "cell_type": "code", - "execution_count": 169, - "id": "eea57a4c", - "metadata": { - "colab": { - "base_uri": "https://localhost:8080/", - "height": 707 + "cell_type": "markdown", + "id": "7f7d69f9", + "metadata": { + "id": "7f7d69f9" + }, + "source": [ + "Copyright © 2024 Gurobi Optimization, LLC\n", + "\n", + "# Price Optimization with Competing Products\n", + "Finding the delicate balance between price and demand is a difficult problem to solve and one that is prevalent in a number of huge industries like retail, e-commerce, ticketing, and hospitality.\n", + "\n", + "It just so happens that this is a problem that data scientists have recently been getting better and better at addressing. Still, a key piece is missing: How can I make the right pricing decision given all the other constraints and business rules that exist for this problem? That's where optimization comes in.\n", + "\n", + "In this scenario, we have a product that has several categories, but we have a limited amount of \"space\" for our products. This could be shelf or warehouse space in retail, seats for events, or for airline ticketing, or rooms at a hotel.\n", + "\n", + "This problem considers two similar products that are offered where it's our job to use the data available with some optimization know-how to determine the optimal mix of products to offer that maximize revenue while also adhering to a few other business rules. First, we'll create a predictive model to forecast sales based on the prices of each product. Then we'll build an optimization model to find this optimal mix. Finally, we’ll also use the Gurobi-sponsored open-source package Gurobi Machine Learning to seamlessly combine the features of a machine learning model with a decision of an optimization model.\n", + "\n", + "Let’s get started!" + ] + }, + { + "cell_type": "markdown", + "id": "a5846a14", + "metadata": { + "id": "a5846a14" + }, + "source": [ + "## Load required packages\n", + "If you have a Gurobi license you can skip the installation of `gurobipy`, but always make sure you have the [latest version](https://www.gurobi.com/downloads/gurobi-software) available." + ] }, - "id": "eea57a4c", - "outputId": "7412df74-2d7e-478a-f7fc-5008c21077a3" - }, - "outputs": [ { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
p[1]p[2]n[1]
0356.12197.67108.0
1358.05189.6866.0
2340.79260.35130.0
3353.76133.5355.0
4341.37229.8091.0
............
995357.63241.5468.0
996352.58212.9587.0
997355.28189.5094.0
998369.75166.3351.0
999349.31222.07114.0
\n", - "

1000 rows × 3 columns

\n", - "
" + "cell_type": "code", + "execution_count": 1, + "id": "1c3276f0", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "1c3276f0", + "outputId": "8b9efd0a-df20-42e8-99b7-28eb2ae2b3e5" + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Requirement already satisfied: gurobipy in /usr/local/lib/python3.10/dist-packages (11.0.0)\n" + ] + } ], - "text/plain": [ - " p[1] p[2] n[1]\n", - "0 356.12 197.67 108.0\n", - "1 358.05 189.68 66.0\n", - "2 340.79 260.35 130.0\n", - "3 353.76 133.53 55.0\n", - "4 341.37 229.80 91.0\n", - ".. ... ... ...\n", - "995 357.63 241.54 68.0\n", - "996 352.58 212.95 87.0\n", - "997 355.28 189.50 94.0\n", - "998 369.75 166.33 51.0\n", - "999 349.31 222.07 114.0\n", - "\n", - "[1000 rows x 3 columns]" + "source": [ + "%pip install gurobipy" ] - }, - "execution_count": 169, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "df = pd.read_csv('https://raw.githubusercontent.com/Gurobi/modeling-examples/master/pricing_competing_products/price_value_data.csv')\n", - "df" - ] - }, - { - "cell_type": "markdown", - "id": "b65707ac", - "metadata": {}, - "source": [ - "### What's in the data?\n", - "The data contains three columns:\n", - "1. `p[1]` is the price (in dollars) of the first category (let's call it Category 1).\n", - "2. `p[2]` is the price (in dollars) of the second category (Category 2).\n", - "3. `n[1]` is the number of the items sold that are of Category 1. \n", - "\n", - "We don't see a column for `n[2]`, which would be the number of items sold that are Category 2. Here is where we make a **pretty big assumption** that we will sell all of the items. This makes our decision to be how to divvy up the limited space we have in order to maximize our revenue. \n", - "The data was created to have a couple of key characteristics.\n", - "1. As the price of Category 1 goes up, the number sold should decrease, so `p[1]` and `n[1]` have a negative correlation.\n", - "2. As the price of Category 2 goes up, the number sold of Category 1 should increase, so `p[2]` and `n[1]` have a positive correlation.\n", - "\n", - "The correlation plot of the columns of the data is below. " - ] - }, - { - "cell_type": "code", - "execution_count": 170, - "id": "4cb22fe7", - "metadata": { - "colab": { - "base_uri": "https://localhost:8080/", - "height": 451 }, - "id": "4cb22fe7", - "outputId": "b6e2ddb2-81ca-485c-e5cc-04b65b292ad9" - }, - "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Warning: environment still referenced so free is deferred (Continue to use WLS)\n", - "Warning: environment still referenced so free is deferred (Continue to use WLS)\n", - "Warning: environment still referenced so free is deferred (Continue to use WLS)\n" - ] + "cell_type": "code", + "execution_count": 2, + "id": "3024295a", + "metadata": { + "id": "3024295a" + }, + "outputs": [], + "source": [ + "import gurobipy as gp\n", + "from gurobipy import GRB\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import seaborn as sns\n", + "import matplotlib.pyplot as plt\n", + "import warnings\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn import tree\n", + "\n", + "warnings.filterwarnings(\"ignore\")" + ] }, { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAABE8AAAHBCAYAAACYIJUYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABRUUlEQVR4nO3deVjUVf//8deAMLgxiqzmbuSGqWACmqWVuC8tLpmoqZSVWqJ33VTeLi3cbWamWS5lpqaVkS0uoanpV9xF00ytNNRA1BS3BMT5/eHPuRuBgWEGkeH5uK7Pdfk5n3MO78/8Meqb9znHYDabzQIAAAAAAECe3Eo6AAAAAAAAgJsZyRMAAAAAAAAbSJ4AAAAAAADYQPIEAAAAAADABpInAAAAAAAANpA8AQAAAAAAsIHkCQAAAAAAgA0kTwAAAAAAAGwgeQIAAAAAAGADyRMAACTt3r1bjz76qOrWrSsvLy9VqlRJoaGhev311/XXX3+VdHgWa9eulcFg0Nq1a+0e+/PPP2vChAk6fPhwrmeDBw9WnTp1HI4PAADAFZE8AQCUebNmzVJYWJi2bt2qf/3rX1qxYoUSEhLUu3dvvf/++xo6dGhJh+gUP//8syZOnJhn8mTcuHFKSEi48UEBAACUAuVKOgAAAEpSUlKSnnjiCXXo0EFfffWVjEaj5VmHDh00ZswYrVixwuGfc/HiRVWoUCFXe05Oji5fvmz1c0tC/fr1S/TnAwAA3MyoPAEAlGmvvvqqDAaDZs6cmWcCw9PTUz169JAkXblyRa+//roaNmwoo9Eof39/DRw4UEePHrUa065dO4WEhOjHH39U69atVaFCBQ0ZMkSHDx+WwWDQ66+/rpdffll169aV0WjUmjVrJEnbtm1Tjx495OPjIy8vL7Vo0UKfffZZge+wbds29evXT3Xq1FH58uVVp04dPfzww/rjjz8sfebOnavevXtLktq3by+DwSCDwaC5c+dKynvZzqVLlxQXF6e6devK09NTt9xyi5566imdOXPGql+dOnXUrVs3rVixQqGhoSpfvrwaNmyoDz/80KrfxYsXNXbsWMvSKB8fH7Vs2VKffvppge8IAABQkqg8AQCUWTk5Ofrhhx8UFhammjVrFtj/iSee0MyZMzVixAh169ZNhw8f1rhx47R27Vrt2LFDvr6+lr6pqakaMGCAnn32Wb366qtyc/vf7yumTp2q2267TW+++aa8vb0VHBysNWvWqFOnTgoPD9f7778vk8mkRYsWqW/fvrp48aIGDx6cb1yHDx9WgwYN1K9fP/n4+Cg1NVUzZszQHXfcoZ9//lm+vr7q2rWrXn31VT3//POaPn26QkNDJeVfcWI2m9WrVy+tXr1acXFxatu2rXbv3q3x48crKSlJSUlJVsmmXbt2acyYMfr3v/+tgIAAzZ49W0OHDtWtt96qu+66S5IUGxurTz75RC+//LJatGihCxcuaM+ePTp16lSBnz0AAECJMgMAUEalpaWZJZn79etXYN99+/aZJZmffPJJq/bNmzebJZmff/55S9vdd99tlmRevXq1Vd9Dhw6ZJZnr169vzsrKsnrWsGFDc4sWLczZ2dlW7d26dTMHBQWZc3JyzGaz2bxmzRqzJPOaNWvyjfXy5cvm8+fPmytWrGh+5513LO2ff/55vmMHDRpkrl27tuV+xYoVZknm119/3arf4sWLzZLMM2fOtLTVrl3b7OXlZf7jjz8sbX///bfZx8fH/Pjjj1vaQkJCzL169co3bgAAgJsVy3YAACiEa0trrq8AadWqlRo1aqTVq1dbtVetWlX33HNPnnP16NFDHh4elvtff/1Vv/zyix555BFJ0uXLly1Xly5dlJqaqv379+cb2/nz5/Xcc8/p1ltvVbly5VSuXDlVqlRJFy5c0L59+4ryuvrhhx8k5X7f3r17q2LFirnet3nz5qpVq5bl3svLS7fddpvV0qFWrVpp+fLl+ve//621a9fq77//LlJsAAAANxrLdgAAZZavr68qVKigQ4cOFdj32tKSoKCgXM+qV69ulSTIr19+z44fPy5JGjt2rMaOHZvnmJMnT+Y7X//+/bV69WqNGzdOd9xxh7y9vWUwGNSlS5ciJyhOnTqlcuXKyc/Pz6rdYDAoMDAw11KbatWq5ZrDaDRa/fypU6eqRo0aWrx4sV577TV5eXmpY8eOeuONNxQcHFykOAEAAG4EkicAgDLL3d1d9957r5YvX66jR4+qRo0a+fa9lhxITU3N1e/PP/+02u9EuppkyM/1z66NjYuL0wMPPJDnmAYNGuTZnpGRoW+//Vbjx4/Xv//9b0t7Zmam/vrrr3xjKEi1atV0+fJlnThxwiqBYjablZaWpjvuuMPuOStWrKiJEydq4sSJOn78uKUKpXv37vrll1+KHCsAAEBxY9kOAKBMi4uLk9lsVkxMjLKysnI9z87O1jfffGNZgjN//nyr51u3btW+fft07733FjmGBg0aKDg4WLt27VLLli3zvCpXrpznWIPBILPZnOukoNmzZysnJ8eq7VqfwlSjXHuf6993yZIlunDhgkPvK0kBAQEaPHiwHn74Ye3fv18XL150aD4AAIDiROUJAKBMi4yM1IwZM/Tkk08qLCxMTzzxhJo0aaLs7Gzt3LlTM2fOVEhIiBISEvTYY4/p3XfflZubmzp37mw5badmzZoaPXq0Q3F88MEH6ty5szp27KjBgwfrlltu0V9//aV9+/Zpx44d+vzzz/Mc5+3trbvuuktvvPGGfH19VadOHa1bt05z5sxRlSpVrPqGhIRIkmbOnKnKlSvLy8tLdevWzXPJTYcOHdSxY0c999xzOnv2rNq0aWM5badFixaKjo62+x3Dw8PVrVs33X777apatar27dunTz75RJGRkapQoYLd8wEAANwoJE8AAGVeTEyMWrVqpbfffluvvfaa0tLS5OHhodtuu039+/fXiBEjJEkzZsxQ/fr1NWfOHE2fPl0mk0mdOnVSfHx8ngkIe7Rv315btmzRK6+8omeeeUanT59WtWrV1LhxY/Xp08fm2IULF+rpp5/Ws88+q8uXL6tNmzZKTExU165drfrVrVtXU6ZM0TvvvKN27dopJydHH330UZ7HIBsMBn311VeaMGGCPvroI73yyivy9fVVdHS0Xn311VyVLoVxzz336Ouvv9bbb7+tixcv6pZbbtHAgQP1wgsv2D0XAADAjWQwm83mkg4CAAAAAADgZsWeJwAAAAAAADaQPAEAAAAAALCB5AkAAAAAAIANJE8AAAAAAIDT/fjjj+revbuqV69u2Yy+IOvWrVNYWJi8vLxUr149vf/++7n6LFmyRI0bN5bRaFTjxo2VkJBQDNFbI3kCAAAAAACc7sKFC2rWrJmmTZtWqP6HDh1Sly5d1LZtW+3cuVPPP/+8Ro0apSVLllj6JCUlqW/fvoqOjtauXbsUHR2tPn36aPPmzcX1GpI4bQcAAAAAABQzg8GghIQE9erVK98+zz33nL7++mvt27fP0jZ8+HDt2rVLSUlJkqS+ffvq7NmzWr58uaVPp06dVLVqVX366afFFj+VJwAAAAAAoFAyMzN19uxZqyszM9MpcyclJSkqKsqqrWPHjtq2bZuys7Nt9tm4caNTYshPuWKd3R57lxTcBwBKkTpdxpZ0CADgVKMWfVXSIQCAU8VGNivpEG4MJ/5/O/7znzRx4kSrtvHjx2vChAkOz52WlqaAgACrtoCAAF2+fFknT55UUFBQvn3S0tIc/vm23DzJEwAAAAAAcFOLi4tTbGysVZvRaHTa/AaDwer+2k4j/2zPq8/1bc5G8gQAAAAAABSK0Wh0arLknwIDA3NVkKSnp6tcuXKqVq2azT7XV6M4G3ueAAAAAADgwsw5OU67ilNkZKQSExOt2r7//nu1bNlSHh4eNvu0bt26WGOj8gQAAAAAADjd+fPn9euvv1ruDx06pOTkZPn4+KhWrVqKi4vTsWPHNG/ePElXT9aZNm2aYmNjFRMTo6SkJM2ZM8fqFJ2nn35ad911l1577TX17NlTS5cu1apVq7Rhw4ZifReSJwAAAAAAuLKcyyXyY7dt26b27dtb7q/tlTJo0CDNnTtXqampSklJsTyvW7euli1bptGjR2v69OmqXr26pk6dqgcffNDSp3Xr1lq0aJFefPFFjRs3TvXr19fixYsVHh5erO9iMF/bfaWkcdoOABfDaTsAXA2n7QBwNWXltJ0r2+c5bS63sIFOm6s0Yc8TAAAAAAAAG1i2AwAAAACAKyvmjV7LApInAAAAAAC4MHMJ7XniSli2AwAAAAAAYAOVJwAAAAAAuDIqTxxG8gQAAAAAABdmvkLyxFEs2wEAAAAAALCByhMAAAAAAFwZp+04jOQJAAAAAAAujNN2HMeyHQAAAAAAABuoPAEAAAAAwJVReeIwkicAAAAAALgw8xX2PHEUy3YAAAAAAABsoPIEAAAAAAAXxoaxjiN5AgAAAACAKyN54jCW7QAAAAAAANhA5QkAAAAAAC6MDWMdR+UJAAAAAACADSRPAAAAAAAAbGDZDgAAAAAArowNYx1G8gQAAAAAABfGUcWOY9kOAAAAAACADVSeAAAAAADgyqg8cRjJEwAAAAAAXBhHFTuOZTsAAAAAAAA2UHkCAAAAAIArY9mOw0ieAAAAAADgwsw5LNtxFMt2AAAAAAAAbKDyBAAAAAAAF2Zm2Y7DSJ4AAAAAAODKrpA8cRTLdgAAAAAAAGyg8gQAAAAAABfGhrGOI3kCAAAAAIArI3niMJbtAAAAAAAA2EDlCQAAAAAALozTdhxH8gQAAAAAAFfGsh2HsWwHAAAAAADABipPAAAAAABwYZy24ziSJwAAAAAAuDDzFZInjmLZDgAAAAAAgA0kTwAAAAAAcGU5Oc67iuC9995T3bp15eXlpbCwMK1fvz7fvoMHD5bBYMh1NWnSxNJn7ty5efa5dOlSkeIrDJInAAAAAACgWCxevFjPPPOMXnjhBe3cuVNt27ZV586dlZKSkmf/d955R6mpqZbryJEj8vHxUe/eva36eXt7W/VLTU2Vl5dXsb0HyRMAAAAAAFAsJk+erKFDh2rYsGFq1KiRpkyZopo1a2rGjBl59jeZTAoMDLRc27Zt0+nTp/Xoo49a9TMYDFb9AgMDi/U9SJ4AAAAAAODCzDk5TrvskZWVpe3btysqKsqqPSoqShs3bizUHHPmzNF9992n2rVrW7WfP39etWvXVo0aNdStWzft3LnTrtjsxWk7AAAAAAC4MHPOFafNlZmZqczMTKs2o9Eoo9GYq+/JkyeVk5OjgIAAq/aAgAClpaUV+LNSU1O1fPlyLVy40Kq9YcOGmjt3rpo2baqzZ8/qnXfeUZs2bbRr1y4FBwcX4a0KRuUJAAAAAAAolPj4eJlMJqsrPj7e5hiDwWB1bzabc7XlZe7cuapSpYp69epl1R4REaEBAwaoWbNmatu2rT777DPddtttevfdd+1+n8Ki8gQAAAAAAFfmxMqTuHFxio2NtWrLq+pEknx9feXu7p6ryiQ9PT1XNcr1zGazPvzwQ0VHR8vT09NmXzc3N91xxx06ePBgId6gaAqdPLn+wymMF198UT4+PnaPAwAAAAAAzmHvXiW25LdEJy+enp4KCwtTYmKi7r//fkt7YmKievbsaXPsunXr9Ouvv2ro0KEF/hyz2azk5GQ1bdq0UHEVRaGTJ1OmTFFkZGSBGZ9rNmzYoBEjRpA8AQAAAACgjIqNjVV0dLRatmypyMhIzZw5UykpKRo+fLgkKS4uTseOHdO8efOsxs2ZM0fh4eEKCQnJNefEiRMVERGh4OBgnT17VlOnTlVycrKmT59ebO9h17KdhIQE+fv7F6pv5cqVixQQAAAAAABwHnOOucR+dt++fXXq1ClNmjRJqampCgkJ0bJlyyyn56SmpiolJcVqTEZGhpYsWaJ33nknzznPnDmjxx57TGlpaTKZTGrRooV+/PFHtWrVqtjew2A2mwv1KX788cfq169foctzFi5cqJ49e6pixYqFi2TvksL1A4BSok6XsSUdAgA41ahFX5V0CADgVLGRzUo6hBvixJhIp83l91aS0+YqTQpdeTJo0CC7Ju7fv7/dwQAAAAAAANxsOG0HZcLWvYc0Z+l67fntmE6cPqfpzw3QfeGNSzosAMjTM888rYf7PyyTyaTknckaN+4/Be4e36lzJ40ZE6tatWopJSVFb77xplau/N7y/Mknn1DHTh1Vv359Xbp0STu279B///uafv/99+J+HQBlnNls1vavPte+dauVeeG8/OsF686BQ+VzS02b43av/E4/r/le50+dlFdlb9VrGa5WD/VXuX/swXjh9F/a9Nl8HdmdrJzsLJkCgnT30CfkV6decb8WUKqYnXjaTlnl5szJdu3aJXd3d2dOCTjFxcwsNagTqP/EdC/pUADApuHDH9fQYUP1n/+MV4/uPXXixAnNX/CJzWWwoaEtNG3au0r4MkFdOndRwpcJmjZ9mpo3b27pEx4erk/mfaL7ez2g6AED5V7OXfM+mafy5cvfgLcCUJbtWrZUu1d+pzYDhuiB8fGqYKqi7954WVl//53vmIMb12vL5wsV1rO3+r76tu4eMly/bUnSli8WWvpkXjivr14eJzf3cuoy5nn1eWWyIh4eKM8KFW7EawGlivmK2WlXWeX0ypNCbqEC3FB3hzbQ3aENSjoMACjQkKFDNH3adK1csVKSNGbMWG3btlU9e/bQwoWf5j1myBBt2LBB7703Q5L03nszFB4eriFDHtWoUU9LkgYNGmw15l9jn9WOndvVtGlTbdmypfheCECZZjab9dP3yxTa/X7VaxkuSWof85TmjYrRr5s2qHH7DnmOO/7bAQUEN1Bw5J2SpMp+/ro1vI3SD/1q6ZP83VJVqlZN7Yc9aWmr7Fe4wy0AwF52JU8eeOABm88zMjJkMBgcCggAgLKqZs2a8vf31/r16y1tWVlZ2rx5s8LCwvJNnrQIbaEP53xo1fbjjz/q0SFD8v1Z107FO3PmjOOBA0A+zp1I18WMM6oR8r9NOd09PBTUsLGO/7o/3+RJYHBDHdy4Xum//yr/erfqbPpxpezeqdva3G3pczh5m2qGNFPitMn6c//PqljVR03uiVKjdvcV+3sBpU1JnrbjKuxKnnzzzTfq0KGDAgIC8nyek5NTqHkyMzOVmZlp1WbMypbR08OecAAAcCl+/n6SpBMnTlq1nzh5UjVuuSX/cX5+OnEy9xg/P998x7w47kVt2bJVBw4ccCBiALDtYsYZSVJ5b5NVe3lvk86fOpnHiKtujWijS+fOaukr4yRJV3Jy1PieKLXo1svS51x6un7+IVFNO3VVi+73K/33X/V/Cz6Su4eHVZIFgGQu3H/VYYNdyZNGjRrpwQcf1NChQ/N8npycrG+//bbAeeLj4zVx4kSrtvFP9NaEp/raEw4AAKVaz1499eqrr1juhzx69e9Xs6x/O2QwGApeFmvOPUb5DJn00iQ1athQDz3U2/6gAcCGgxvX68ePZ1ruO4+Ou/qH66vTC/hO+3PfXu345kvdOXCY/OsF62x6mjYu+EjbTVUU1vOh/z/FFfnVra/wh66e8ulbu65OHzuivT98T/IEgNPZlTwJCwvTjh078k2eGI1G1apVq8B54uLiFBsbaz32t2X2hAIAQKm3KnGVkncmW+49//8JEv5+fjqRfsLS7lutmk6ezP83tCdOnJCfn59Vm2+1armqUSRpwsQJuu++e9WnT1+lpaU5+AYAYK12i5Z6qH6w5T7ncrYk6e+MM6pYpaql/e9zZ1XBZMo1/pqtCYsV3PouNbr7XklStZq1lJ15SevnzlRo9wdkcHNThSpVVbV6DatxVarX0O/bNjvzlQCXwLIdx9mVPHn//fdtLs1p1KiRDh06VOA8RqNRRqPRupElOwCAMubChQu6cOGCVVt6erruvLOt9u79WZLk4eGh8PBw/fe//813np07durOtndqzj/2PWl7V1vt2L7dqt/ESRPVsWOU+vV9WEePHHXimwDAVZ7ly8vzH6d4mc1mVTBV0dG9u+Vbu64kKefyZaX+8rPC+zyS7zyXMzNlcLOuVnFzc5PZfLU2zyApMLiBzqT9adUnI+1PVfa1TiYDkK5wUrHD7Eqe5Ep4AKXEhb8zlZJ2ynJ/NP0v7Tv0p0yVKqi6X5WSCwwArvPhnA/11FNP6vDhQzp06LCeGvGk/r70t5Yu/drS563Jb+l4Wppef/2Nq2M++kiffbZYw4c/rsTERHXo0EFt2rRR74f6WMa89PIk9ezRUzExj+nChfOW/VDOnj2Xax8yAHAWg8GgplFdtPObBJkCgmQKCNTObxNUzmjUrRF3Wvr9MHOaKlb1UXjvq0twajcP0+6V38m3Vl351w/W2eNp2vrlYtVu0VJubm6SpKZRXbX0lXHa8c2Xqt+qtdJ//1X71q7WXYMfK5F3BeDaCp08OXv2rLy9vQs98blz5yw7+QMlbc9vxzTwP7Mt9/EfXV0mdn/7UP135EMlFRYA5PL++x/Iy8tLL738kkzeJiUnJyt6wECrCpVbqleX+R+/QtqxfYdGjhylsWPGKHZMrFJSUjRixEglJydb+kRHR0uSFn+2yOrnjR0zVl98saR4XwpAmdasS09dzsrShnmzlXnhgvzr36quY1+wqlA5f+qk1amdoT0elAwGbf1ykS6c/kvlK3urVvMwtXrwYUsf/3q3KmrkWG35YqF2LF2iyn7+at1/kIJbt72h7weUBmwY6ziDucAd6K5yd3dXamqq/P0Ld3a6t7e3kpOTVa9evcJFspd/uAFwLXW6jC3pEADAqUYt+qqkQwAAp4qNbFZwJxeQEt3SaXPV+mSb0+YqTQpdeWI2mzV79mxVqlSpUP2zs7OLHBQAAAAAAMDNotDJk1q1amnWrFmFnjgwMFAeHmwCCwAAAAAASrdCJ08OHz6cZ/u1VT+G689uBwAAAAAAJY7TdhznVtSBc+bMUUhIiLy8vOTl5aWQkBDNnj274IEAAAAAAACliF1HFV8zbtw4vf322xo5cqQiIyMlSUlJSRo9erQOHz6sl19+2alBAgAAAACAouG0HccVKXkyY8YMzZo1Sw8//L+jwnr06KHbb79dI0eOJHkCAAAAAMBN4soVttlwVJGW7eTk5Khly9xHHYWFheny5csOBwUAAAAAAHCzKFLyZMCAAZoxY0au9pkzZ+qRRx5xOCgAAAAAAOAcV6447yqrirRsR7q6Yez333+viIgISdKmTZt05MgRDRw4ULGxsZZ+kydPdjxKAAAAAABQJOx54rgiJU/27Nmj0NBQSdJvv/0mSfLz85Ofn5/27Nlj6cfxxQAAAAAAoLQrUvJkzZo1zo4DAAAAAAAUAzaMdVyRl+0AAAAAAICb3xWW7TisSBvGAgAAAAAAlBVUngAAAAAA4MJYtuM4kicAAAAAALgwM8kTh7FsBwAAAAAAwAYqTwAAAAAAcGFXrpR0BKUflScAAAAAAAA2UHkCAAAAAIALY8NYx5E8AQAAAADAhZE8cRzLdgAAAAAAAGyg8gQAAAAAABeWQ+WJw0ieAAAAAADgwli24ziW7QAAAAAAANhA5QkAAAAAAC7sipnKE0dReQIAAAAAAGADyRMAAAAAAAAbWLYDAAAAAIALu3KlpCMo/UieAAAAAADgwnLY88RhLNsBAAAAAADF5r333lPdunXl5eWlsLAwrV+/Pt++a9eulcFgyHX98ssvVv2WLFmixo0by2g0qnHjxkpISCjWdyB5AgAAAACAC7tyxeC0y16LFy/WM888oxdeeEE7d+5U27Zt1blzZ6WkpNgct3//fqWmplqu4OBgy7OkpCT17dtX0dHR2rVrl6Kjo9WnTx9t3rzZ7vgKi+QJAAAAAAAuLMdscNplr8mTJ2vo0KEaNmyYGjVqpClTpqhmzZqaMWOGzXH+/v4KDAy0XO7u7pZnU6ZMUYcOHRQXF6eGDRsqLi5O9957r6ZMmWJ3fIVF8gQAAAAAADhdVlaWtm/frqioKKv2qKgobdy40ebYFi1aKCgoSPfee6/WrFlj9SwpKSnXnB07dixwTkewYSwAAAAAAC7sihM3jM3MzFRmZqZVm9FolNFozNX35MmTysnJUUBAgFV7QECA0tLS8pw/KChIM2fOVFhYmDIzM/XJJ5/o3nvv1dq1a3XXXXdJktLS0uya0xlIngAAAAAA4MKcedpOfHy8Jk6caNU2fvx4TZgwId8xBoP1zzebzbnarmnQoIEaNGhguY+MjNSRI0f05ptvWpIn9s7pDCRPAAAAAABAocTFxSk2NtaqLa+qE0ny9fWVu7t7roqQ9PT0XJUjtkRERGj+/PmW+8DAQIfntBd7ngAAAAAA4MJyzM67jEajvL29ra78kieenp4KCwtTYmKiVXtiYqJat25d6Ph37typoKAgy31kZGSuOb///nu75rQXlScAAAAAALgwZ+55Yq/Y2FhFR0erZcuWioyM1MyZM5WSkqLhw4dLulrJcuzYMc2bN0/S1ZN06tSpoyZNmigrK0vz58/XkiVLtGTJEsucTz/9tO666y699tpr6tmzp5YuXapVq1Zpw4YNxfYeJE8AAAAAAECx6Nu3r06dOqVJkyYpNTVVISEhWrZsmWrXri1JSk1NVUpKiqV/VlaWxo4dq2PHjql8+fJq0qSJvvvuO3Xp0sXSp3Xr1lq0aJFefPFFjRs3TvXr19fixYsVHh5ebO9hMJvN5mKb3R57lxTcBwBKkTpdxpZ0CADgVKMWfVXSIQCAU8VGNivpEG6IhND7nDbX/TtWOW2u0oTKEwAAAAAAXFjOzVEyUaqxYSwAAAAAAIANVJ4AAAAAAODCclRyG8a6CpInAAAAAAC4MJbtOI5lOwAAAAAAADZQeQIAAAAAgAvLKekAXADJEwAAAAAAXBjJE8exbAcAAAAAAMAGkicAAAAAAAA2sGwHAAAAAAAXxlHFjqPyBAAAAAAAwAYqTwAAAAAAcGE5ZnNJh1DqkTwBAAAAAMCFcdqO41i2AwAAAAAAYAOVJwAAAAAAuDAqTxxH8gQAAAAAABdG8sRxLNsBAAAAAACwgcoTAAAAAABcWI44bcdRJE8AAAAAAHBhLNtx3E2TPKnTZWxJhwAATnV42ZslHQIAOFXnXgNKOgQAcKrYgz+VdAgoJW6a5AkAAAAAAHC+HDPLdhxF8gQAAAAAABfGsh3HcdoOAAAAAACADVSeAAAAAADgwjhtx3EkTwAAAAAAcGEkTxzHsh0AAAAAAAAbqDwBAAAAAMCFsWGs40ieAAAAAADgwjiq2HEs2wEAAAAAALCByhMAAAAAAFwYG8Y6juQJAAAAAAAujOSJ41i2AwAAAAAAYAPJEwAAAAAAABtYtgMAAAAAgAu7wmk7DqPyBAAAAAAAwAYqTwAAAAAAcGFsGOs4kicAAAAAALgwkieOY9kOAAAAAACADVSeAAAAAADgwnLYMNZhJE8AAAAAAHBhLNtxHMt2AAAAAAAAbCB5AgAAAACAC7tiNjvtKor33ntPdevWlZeXl8LCwrR+/fp8+3755Zfq0KGD/Pz85O3trcjISK1cudKqz9y5c2UwGHJdly5dKlJ8hUHyBAAAAAAAF5Yjs9Muey1evFjPPPOMXnjhBe3cuVNt27ZV586dlZKSkmf/H3/8UR06dNCyZcu0fft2tW/fXt27d9fOnTut+nl7eys1NdXq8vLyKtLnUxjseQIAAAAAAIrF5MmTNXToUA0bNkySNGXKFK1cuVIzZsxQfHx8rv5Tpkyxun/11Ve1dOlSffPNN2rRooWl3WAwKDAwsFhj/ycqTwAAAAAAcGHOrDzJzMzU2bNnra7MzMw8f25WVpa2b9+uqKgoq/aoqCht3LixULFfuXJF586dk4+Pj1X7+fPnVbt2bdWoUUPdunXLVZnibCRPAAAAAABwYc7c8yQ+Pl4mk8nqyquCRJJOnjypnJwcBQQEWLUHBAQoLS2tULG/9dZbunDhgvr06WNpa9iwoebOnauvv/5an376qby8vNSmTRsdPHiw6B9SAVi2AwAAAAAACiUuLk6xsbFWbUaj0eYYg8FgdW82m3O15eXTTz/VhAkTtHTpUvn7+1vaIyIiFBERYblv06aNQkND9e6772rq1KmFeQ27kTwBAAAAAMCFFWWj1/wYjcYCkyXX+Pr6yt3dPVeVSXp6eq5qlOstXrxYQ4cO1eeff6777rvPZl83NzfdcccdxVp5wrIdAAAAAABcWI7Z7LTLHp6engoLC1NiYqJVe2Jiolq3bp3vuE8//VSDBw/WwoUL1bVr1wJ/jtlsVnJysoKCguyKzx5UngAAAAAAgGIRGxur6OhotWzZUpGRkZo5c6ZSUlI0fPhwSVeXAR07dkzz5s2TdDVxMnDgQL3zzjuKiIiwVK2UL19eJpNJkjRx4kRFREQoODhYZ8+e1dSpU5WcnKzp06cX23uQPAEAAAAAwIVdceKyHXv17dtXp06d0qRJk5SamqqQkBAtW7ZMtWvXliSlpqYqJSXF0v+DDz7Q5cuX9dRTT+mpp56ytA8aNEhz586VJJ05c0aPPfaY0tLSZDKZ1KJFC/34449q1apVsb2HwWy2s+6mmNSpXbekQwAApzq87M2SDgEAnKpzrwklHQIAONXygz+VdAg3RK8GLZw211f7i/dI4JsVe54AAAAAAADYwLIdAAAAAABc2JWbY8FJqUblCQAAAAAAgA0kTwAAAAAAAGxg2Q4AAAAAAC4spwRP23EVJE8AAAAAAHBhV8xXSjqEUo9lOwAAAAAAADZQeQIAAAAAgAu7wrIdh5E8AQAAAADAheVwVLHDWLYDAAAAAABgA5UnAAAAAAC4MJbtOI7kCQAAAAAALuwKy3YcxrIdAAAAAAAAG6g8AQAAAADAhV0p6QBcAMkTAAAAAABcGMt2HMeyHQAAAAAAABuoPAEAAAAAwIVx2o7jSJ4AAAAAAODCWLbjOJbtAAAAAAAA2EDlCQAAAAAALoxlO44rdPIkNjbW7slffPFF+fj42D0OAAAAAAA4B8kTxxU6eTJlyhRFRkbK09OzUP03bNigESNGkDwBAAAAAAClml3LdhISEuTv71+ovpUrVy5SQAAAAAAAwHmuUHjisEInTz766COZTKZCT/zBBx8oICCgSEEBAAAAAADnYNmO4wqdPBk0aJBdE/fv39/uYAAAAAAAAG42HFWMUu+ZZ57W5i2b9Mv+fVq06FMFBwcXOKZT505KXPW99h/4RYmrvlfHjlFWz5988gkt/for7dn7k7Zt36qZMz9QvXr1iusVAMBuW/ce0vBX5+nOofFq8MDzWrX555IOCQDy1LV/X330w3It3bNNUxMWq0nL0Hz7VvXz1bOTX9OslV/ru/279PgLz+bq416unPqPGK4PVy/T0j3bNP3rLxTWtk1xvgJQ6l2R2WlXWWV38uS7777TsGHD9Oyzz+qXX36xenb69Gndc889TgsOKMjw4Y9r6LCh+s9/xqtH9546ceKE5i/4RBUrVsx3TGhoC02b9q4SvkxQl85dlPBlgqZNn6bmzZtb+oSHh+uTeZ/o/l4PKHrAQLmXc9e8T+apfPnyN+CtAKBgFzOz1KBOoP4T072kQwGAfN3VpaMef+E5LZoxSyN69tbebdv10uwZ8gsKzLO/h6enMv76S4tmzNKhX/bn2WfQ6JHq3PchzZgUr8c799KyRZ9p3HtTVL9xw+J8FQBlnF3Jk4ULF6pnz55KS0tTUlKSWrRooQULFlieZ2Vlad26dU4PEsjPkKFDNH3adK1csVIHDhzQmDFjVd6rvHr27JH/mCFDtGHDBr333gz99tvveu+9Gdr4fxs1ZMijlj6DBg3WF18s0cGDB7Vv3z79a+yzqlHjFjVt2vRGvBYAFOju0AYa3T9KUREhJR0KAOTr/iED9f0XX2rl51/qyG+H9MErr+tEWpq69u+bZ//0Y3/qg5df0+qvvtGFc+fz7HNPz25a/P5sbV23XmlHjuq7hZ9p+/qNemCIfdsMAIA97EqevPnmm3r77bf17bffav369frkk080fPhwzZkzp7jiA/JVs2ZN+fv7a/369Za2rKwsbd68WWFhYfmOaxHaQut/XG/V9uOPPyrUxphrp0edOXPGsaABAADKiHIe5RTcpLF2bNho1b5jw0Y1Dm1e5Hk9PD2VlZlp1ZaVeUlNwloUeU7A1ZnNzrvKKruOKj5w4IC6detmuX/ooYfk6+urHj16KDs7W/fff7/TAwTy4+fvJ0k6ceKkVfuJkydV45Zb8h/n56cTJ3OP8fPzzXfMi+Ne1JYtW3XgwAEHIgYAACg7vKtWlXu5cjp98pRV+5mTp1TVt1qR592+YaMeGDJQe7ZuV2rKETVvHaGIe9vL3d3d0ZABl1WW9ypxFruSJ97e3jp+/Ljq1q1raWvXrp2++eYbdevWTUePHi3UPJmZmcq8LltsNptlMBjsCQdlTM9ePfXqq69Y7oc8OlSSZL7ui8BgMMhcUErUnHtMft8nk16apEYNG+qhh3rbHzQAAEAZd/0/y67+W63o833w8n816uUJmrnya8lsVmrKESUuWaoOD/Z0LFAAsMGu5EmrVq20fPlyRUREWLXffffdlgRKYcTHx2vixIlWbSZvk6pUqWpPOChjViWuUvLOZMu9p6enJMnfz08n0k9Y2n2rVdPJ6ypL/unEiRPy8/OzavOtVi1XNYokTZg4Qffdd6/69OmrtLQ0B98AAACg7Dh7+rRyLl+Wj591lYmpmo/OnDqVz6iCZfx1Wi89+bQ8PD3lXbWKTh1P15B/jdbxo8ccDRlwWdSdOM6uPU9Gjx4tLy+vPJ+1a9dO3377rQYOHFjgPHFxccrIyLC6TKYq9oSCMujChQv6448/LNfBgweVnp6uO+9sa+nj4eGh8PBwbd++Pd95du7YqTvb3mnV1vauttpx3ZiJkyaqU6eO6v/wIzp6pHBVVQAAALjqcvZlHdz7s1q0ibRqD20TqZ93JDs8f3ZWlk4dT5d7uXJq0/E+Ja1a4/CcgKviqGLH2VV5cvfdd+vuu+/O93m7du3Url27AucxGo0yGo1WbSzZQVF8OOdDPfXUkzp8+JAOHTqsp0Y8qb8v/a2lS7+29Hlr8ls6npam119/4+qYjz7SZ58t1vDhjysxMVEdOnRQmzZt1PuhPpYxL708ST179FRMzGO6cOG8ZT+Us2fP5VpyBgAl4cLfmUpJ+99vbo+m/6V9h/6UqVIFVferUnKBAcA/JHw4T2PfiNfBPXu1b+cude7bW35BQVr26WeSpMFjnla1AH+99ewLljH1GjWQJHlVqCCTj4/qNWqgy9nZSvn1d0lSg2ZNVS3AX7/v269qAf4aMPIJGdzc9MWsj278CwIoMwqdPDl79qy8vb0LPfG5c+csJ5QAxeX99z+Ql5eXXnr5JZm8TUpOTlb0gIG6cOGCpc8t1avLfOWK5X7H9h0aOXKUxo4Zo9gxsUpJSdGIESOVnJxs6RMdHS1JWvzZIqufN3bMWH3xxZLifSkAKIQ9vx3TwP/MttzHf7RMknR/+1D9d+RDJRUWAFj5cdlKVa5SRf2fGi4ffz8dPvCr/hPzpNL/TJUk+fj7yb96kNWY6V9/YfnzbU2bqH2Prjp+9JgGt+8kSfI0GjVo9EgF1qyhvy9c1NZ16/XGv57XhXPnbtyLAaVM2a0XcR6DucCdNa9yd3dXamqq/P39CzWxt7e3kpOTVa9evUL1r1O7bsGdAKAUObzszZIOAQCcqnOvCSUdAgA41fKDP5V0CDfEbbXrOG2uA38cdtpcpUmhK0/MZrNmz56tSpUqFap/dnZ2kYMCAAAAAAC4WRQ6eVKrVi3NmjWr0BMHBgbKw8OjSEEBAAAAAADnKMsbvTpLoZMnhw8fzrP92qofNnwFAAAAAODmQ+rEcXYdVfxPc+bMUUhIiLy8vOTl5aWQkBDNnj274IEAAAAAAAClSJGSJ+PGjdPTTz+t7t276/PPP9fnn3+u7t27a/To0XrxxRedHSMAAAAAACgisxOvonjvvfdUt25deXl5KSwsTOvXr7fZf926dQoLC5OXl5fq1aun999/P1efJUuWqHHjxjIajWrcuLESEhKKGF3hFCl5MmPGDM2aNUvx8fHq0aOHevToofj4eM2cOTPPlwIAAAAAACWjJJMnixcv1jPPPKMXXnhBO3fuVNu2bdW5c2elpKTk2f/QoUPq0qWL2rZtq507d+r555/XqFGjtGTJEkufpKQk9e3bV9HR0dq1a5eio6PVp08fbd68uQgRFk6hjyr+p6pVq2rLli0KDg62aj9w4IBatWqlM2fO2B0IRxUDcDUcVQzA1XBUMQBXU1aOKq5bu7bT5jr0xx929Q8PD1doaKhmzJhhaWvUqJF69eql+Pj4XP2fe+45ff3119q3b5+lbfjw4dq1a5eSkpIkSX379tXZs2e1fPlyS59OnTqpatWq+vTTT+19pUIpUuXJgAEDrF78mpkzZ+qRRx5xOCgAAAAAAOAcJVV5kpWVpe3btysqKsqqPSoqShs3bsxzTFJSUq7+HTt21LZt25SdnW2zT35zOkOhT9u53pw5c/T9998rIiJCkrRp0yYdOXJEAwcOVGxsrKXf5MmTHY8SAAAAAACUuMzMTGVmZlq1GY1GGY3GXH1PnjypnJwcBQQEWLUHBAQoLS0tz/nT0tLy7H/58mWdPHlSQUFB+fbJb05nKFLyZM+ePQoNDZUk/fbbb5IkPz8/+fn5ac+ePZZ+HF8MAAAAAIDriI+P18SJE63axo8frwkTJuQ75vrcgNlstpkvyKv/9e32zumoIiVP1qxZ4+w4AAAAAABAsXBeUiEuLs5qtYmkPKtOJMnX11fu7u65KkLS09NzVY5cExgYmGf/cuXKqVq1ajb75DenMxRpzxMAAAAAAFBaGJx2GY1GeXt7W135JU88PT0VFhamxMREq/bExES1bt06zzGRkZG5+n///fdq2bKlPDw8bPbJb05nKPKeJwAAAAAAALbExsYqOjpaLVu2VGRkpGbOnKmUlBQNHz5c0tVKlmPHjmnevHmSrp6sM23aNMXGxiomJkZJSUmaM2eO1Sk6Tz/9tO666y699tpr6tmzp5YuXapVq1Zpw4YNxfYeJE8AAAAAAECx6Nu3r06dOqVJkyYpNTVVISEhWrZsmWr//+OTU1NTlZKSYulft25dLVu2TKNHj9b06dNVvXp1TZ06VQ8++KClT+vWrbVo0SK9+OKLGjdunOrXr6/FixcrPDy82N7DYL6280oJq1O7bkmHAABOdXjZmyUdAgA4VedeE0o6BABwquUHfyrpEG6IOrXrOW2uw3/87rS5ShP2PAEAAAAAALCBZTsAAAAAALiy4jvBt8wgeQIAAAAAgEtj0Ymj+AQBAAAAAABsoPIEAAAAAAAXZmDdjsNIngAAAAAA4MoMJE8cxbIdAAAAAAAAG6g8AQAAAADAhbFsx3EkTwAAAAAAcGksOnEUnyAAAAAAAIANVJ4AAAAAAODCDGwY6zCSJwAAAAAAuDIDi04cxScIAAAAAABgA5UnAAAAAAC4MAN1Ew4jeQIAAAAAgAtjzxPHkX4CAAAAAACwgcoTAAAAAABcGRvGOozkCQAAAAAALsxA8sRhfIIAAAAAAAA2UHkCAAAAAIAL47Qdx5E8AQAAAADAhbFsx3F8ggAAAAAAADaQPAEAAAAAALCBZTsAAAAAALgwg8G9pEMo9ag8AQAAAAAAsIHKEwAAAAAAXBgbxjqO5AkAAAAAAC6M5Inj+AQBAAAAAABsoPIEAAAAAAAXxoaxjiN5AgAAAACAC2PZjuP4BAEAAAAAAGyg8gQAAAAAABfGsh3HkTwBAAAAAMCFkTxxHMt2AAAAAAAAbKDyBAAAAAAAF+bGhrEOI3kCAAAAAIALY9mO40g/AQAAAAAA2EDlCQAAAAAALozKE8eRPAEAAAAAwIWRPHEcy3YAAAAAAABsoPIEAAAAAAAXZnCj8sRRVJ4AAAAAAODC3AzuTruKy+nTpxUdHS2TySSTyaTo6GidOXMm3/7Z2dl67rnn1LRpU1WsWFHVq1fXwIED9eeff1r1a9eunQwGg9XVr18/u+MjeQIAAAAAAEpU//79lZycrBUrVmjFihVKTk5WdHR0vv0vXryoHTt2aNy4cdqxY4e+/PJLHThwQD169MjVNyYmRqmpqZbrgw8+sDs+lu0AAAAAAODCbvYNY/ft26cVK1Zo06ZNCg8PlyTNmjVLkZGR2r9/vxo0aJBrjMlkUmJiolXbu+++q1atWiklJUW1atWytFeoUEGBgYEOxXjTJE9GLfqqpEMAAKfq3GtASYcAAE61/KsJJR0CAKCEZWZmKjMz06rNaDTKaDQWec6kpCSZTCZL4kSSIiIiZDKZtHHjxjyTJ3nJyMiQwWBQlSpVrNoXLFig+fPnKyAgQJ07d9b48eNVuXJlu2Jk2Q4AAAAAACiU+Ph4y74k1674+HiH5kxLS5O/v3+udn9/f6WlpRVqjkuXLunf//63+vfvL29vb0v7I488ok8//VRr167VuHHjtGTJEj3wwAN2x3jTVJ4AAAAAAADnc+aynbi4OMXGxlq15Vd1MmHCBE2cONHmfFu3bpUkGQyGXM/MZnOe7dfLzs5Wv379dOXKFb333ntWz2JiYix/DgkJUXBwsFq2bKkdO3YoNDS0wLmvIXkCAAAAAIALMxic919/e5bojBgxosCTberUqaPdu3fr+PHjuZ6dOHFCAQEBNsdnZ2erT58+OnTokH744QerqpO8hIaGysPDQwcPHiR5AgAAAAAASpavr698fX0L7BcZGamMjAxt2bJFrVq1kiRt3rxZGRkZat26db7jriVODh48qDVr1qhatWoF/qy9e/cqOztbQUFBhX8RsecJAAAAAAAuzc3g7rSrODRq1EidOnVSTEyMNm3apE2bNikmJkbdunWz2iy2YcOGSkhIkCRdvnxZDz30kLZt26YFCxYoJydHaWlpSktLU1ZWliTpt99+06RJk7Rt2zYdPnxYy5YtU+/evdWiRQu1adPGrhipPAEAAAAAwIUZ3G7uo4qlqyfijBo1SlFRUZKkHj16aNq0aVZ99u/fr4yMDEnS0aNH9fXXX0uSmjdvbtVvzZo1ateunTw9PbV69Wq98847On/+vGrWrKmuXbtq/Pjxcne37zMheQIAAAAAAEqUj4+P5s+fb7OP2Wy2/LlOnTpW93mpWbOm1q1b55T4SJ4AAAAAAODCnLlhbFnFJwgAAAAAgAtz5lHFZRUbxgIAAAAAANhA5QkAAAAAAC6MZTuO4xMEAAAAAMCFFdcRw2UJy3YAAAAAAABsoPIEAAAAAAAXZnDjv/6O4hMEAAAAAMCFseeJ41i2AwAAAAAAYAPpJwAAAAAAXJiBDWMdRvIEAAAAAAAXxrIdx7FsBwAAAAAAwAbSTwAAAAAAuDBO23EclScAAAAAAAA2kH4CAAAAAMCFseeJ4/gEAQAAAABwZSRPHMayHQAAAAAAABtIngAAAAAAANhA7Q4AAAAAAC6M03YcR+UJAAAAAACADaSfAAAAAABwYZy24zg+QQAAAAAAXBnLdhzGsh0AAAAAAAAbSD8BAAAAAODKDO4lHUGpR/IEAAAAAAAXxmk7jmPZDgAAAAAAgA2knwAAAAAAcGWctuMwPkEAAAAAAFyYmWU7DmPZDgAAAAAAgA2knwAAAAAAcGVunLbjKJInAAAAAAC4MpInDmPZDgAAAAAAgA1UngAAAAAA4MLMVJ44jOQJAAAAAAAujOSJ41i2AwAAAAAAYAOVJwAAAAAAuDIqTxxG8gQAAAAAABdmdmPRiaP4BAEAAAAAAGyg8gQAAAAAABfGhrGOo/IEAAAAAADABpInAAAAAAAANpA8AQAAAADAhV1xd3PaVVxOnz6t6OhomUwmmUwmRUdH68yZMzbHDB48WAaDweqKiIiw6pOZmamRI0fK19dXFStWVI8ePXT06FG74yN5AgAAAACACzO7uTntKi79+/dXcnKyVqxYoRUrVig5OVnR0dEFjuvUqZNSU1Mt17Jly6yeP/PMM0pISNCiRYu0YcMGnT9/Xt26dVNOTo5d8bFhLAAAAAAAKDH79u3TihUrtGnTJoWHh0uSZs2apcjISO3fv18NGjTId6zRaFRgYGCezzIyMjRnzhx98sknuu+++yRJ8+fPV82aNbVq1Sp17Nix0DFSeQIAAAAAgAtzZuVJZmamzp49a3VlZmY6FF9SUpJMJpMlcSJJERERMplM2rhxo82xa9eulb+/v2677TbFxMQoPT3d8mz79u3Kzs5WVFSUpa169eoKCQkpcN7rFbryZOrUqXZNLEmPPvqoKleubPc4AAAAAADgHFecuNwmPj5eEydOtGobP368JkyYUOQ509LS5O/vn6vd399faWlp+Y7r3Lmzevfurdq1a+vQoUMaN26c7rnnHm3fvl1Go1FpaWny9PRU1apVrcYFBATYnDcvhU6ePPPMM6pRo4bc3Qt3PvSRI0fUrVs3kicAAAAAALiIuLg4xcbGWrUZjcY8+06YMCFXouV6W7dulSQZDIZcz8xmc57t1/Tt29fy55CQELVs2VK1a9fWd999pwceeCDfcQXNmxe79jzZtm1bntmgvJA0AQAAAACg5JmdeEqO0WjMN1lyvREjRqhfv342+9SpU0e7d+/W8ePHcz07ceKEAgICCh1bUFCQateurYMHD0qSAgMDlZWVpdOnT1tVn6Snp6t169aFnleyI3kyfvx4VapUqdATP//88/Lx8bErGAAAAAAA4FxmN/uqLJzF19dXvr6+BfaLjIxURkaGtmzZolatWkmSNm/erIyMDLuSHKdOndKRI0cUFBQkSQoLC5OHh4cSExPVp08fSVJqaqr27Nmj119/3a53KXT6afz48apQoUKhJ46Li1OVKlXsCgYAAAAAAJQtjRo1UqdOnRQTE6NNmzZp06ZNiomJUbdu3axO2mnYsKESEhIkSefPn9fYsWOVlJSkw4cPa+3aterevbt8fX11//33S5JMJpOGDh2qMWPGaPXq1dq5c6cGDBigpk2bWk7fKSyOKkapZjabtf2rz7Vv3WplXjgv/3rBunPgUPncUtPmuN0rv9PPa77X+VMn5VXZW/VahqvVQ/1VztPT0ufC6b+06bP5OrI7WTnZWTIFBOnuoU/Ir0694n4tAGVY1/599dCwwfLx99MfB3/TB6+8pr3bduTZt6qfr2Li/qXgJo1UvU5tfT1vgT54xfq3KO7lyqnv8GG67/4eqhbgr6O/H9aHb7yt7ev/70a8DgAU2ta9hzRn6Xrt+e2YTpw+p+nPDdB94Y1LOizAJVxxL5nKE3ssWLBAo0aNspyM06NHD02bNs2qz/79+5WRkSFJcnd3108//aR58+bpzJkzCgoKUvv27bV48WKrbUTefvttlStXTn369NHff/+te++9V3Pnzi30fq7XODV5sm/fPnXt2lW///67M6cF8rVr2VLtXvmd2g17UlUCg7Tj6y/13Rsvq2/8FHmWL5/nmIMb12vL5wt199AnFHjrbTpzPFVrZ78nSWrdf7AkKfPCeX318jhVb9REXcY8r/KVvZVx4rg87ai+AgB73dWlox5/4TlNn/Cyft6xU1369dZLs2fo8c49dSI1947wHp6eyvjrLy2aMUv3Pxqd55yDRo9U+x5dNfXFiTry+yGFtW2tce9N0Zi+0frt51+K+5UAoNAuZmapQZ1APXBPqEa+vrCkwwFcSkkt27GHj4+P5s+fb7OP2Wy2/Ll8+fJauXJlgfN6eXnp3Xff1bvvvutQfM7bNUZSVlaW/vjjD2dOCeTLbDbrp++XKbT7/arXMlw+NWqpfcxTupyZqV83bch33PHfDigguIGCI+9UZT9/1QxpplvD2+jE4f8l/ZK/W6pK1aqp/bAn5V/vVlX281eNxk1l8g+8Ea8GoIy6f8hAff/Fl1r5+Zc68tshffDK6zqRlqau/fvm2T/92J/64OXXtPqrb3Th3Pk8+9zTs5sWvz9bW9etV9qRo/pu4Wfavn6jHhgyqDhfBQDsdndoA43uH6WoiJCSDgUAcrGr8uT644iud+LECYeCAexx7kS6LmacUY2QZpY2dw8PBTVsrOO/7lfj9h3yHBcY3FAHN65X+u+/yr/erTqbflwpu3fqtjZ3W/ocTt6mmiHNlDhtsv7c/7MqVvVRk3ui1KidfeviAKCwynmUU3CTxvr8gzlW7Ts2bFTj0OZFntfD01NZmZlWbVmZl9QkrEWR5wQAAKVLaag8udnZlTx555131Lx5c3l7e+f5/Pz5vH/rBRSHixlnJEnlvU1W7eW9TTp/6mS+426NaKNL585q6SvjJElXcnLU+J4otejWy9LnXHq6fv4hUU07dVWL7vcr/fdf9X8LPpK7h4dVkgUAnMW7alW5lyun0ydPWbWfOXlKVX2rFXne7Rs26oEhA7Vn63alphxR89YRiri3vd3rfAEAQOll5q99h9mVPAkODtbo0aM1YMCAPJ8nJycrLCyswHkyMzOVed1vwS5nZVlt1glc7+DG9frx45mW+86j467+wXBdFvUf6+Dy8ue+vdrxzZe6c+Aw+dcL1tn0NG1c8JG2m6oorOdD/3+KK/KrW1/hD/WXJPnWrqvTx45o7w/fkzwBUKyu/wozGAwFfa3Z9MHL/9Wolydo5sqvJbNZqSlHlLhkqTo82NOxQAEAAMoQu5InYWFh2r59e77Jk6v/wCv4X3jx8fGaOHGiVVvUkMfVcdgT9oSDMqZ2i5Z6qH6w5T7ncrYk6e+MM6pYpaql/e9zZ1XBZMo1/pqtCYsV3PouNbr7XklStZq1lJ15SevnzlRo9wdkcHNThSpVVbV6DatxVarX0O/bNjvzlQDA4uzp08q5fFk+ftZVJqZqPjpz6lQ+owqW8ddpvfTk0/Lw9JR31So6dTxdQ/41WsePHnM0ZAAAUEqwbMdxdiVP3nrrrVwVI//UrFkzXblypcB54uLicu2f8v7O/faEgjLIs3x5qxN0zGazKpiq6Oje3fKtXVeSlHP5slJ/+VnhfR7Jd57LmZkyXPfl4ebmJrPZLLMkg6TA4AY6k/anVZ+MtD9V2dfPae8DAP90OfuyDu79WS3aRGpj4g+W9tA2kUpatcbh+bOzsnTqeLrcy5VTm4736cdlBe9ODwAAXIRTj4opm+xKngQGOuekEaPRKKPRaB0IS3ZgJ4PBoKZRXbTzmwSZAoJkCgjUzm8TVM5o1K0Rd1r6/TBzmipW9VF476tLcGo3D9Puld/Jt1Zd+dcP1tnjadr65WLVbtFSbm5Xv1WaRnXV0lfGacc3X6p+q9ZK//1X7Vu7WncNfqxE3hVA2ZDw4TyNfSNeB/fs1b6du9S5b2/5BQVp2aefSZIGj3la1QL89dazL1jG1GvUQJLkVaGCTD4+qteogS5nZyvl16sniDVo1lTVAvz1+779qhbgrwEjn5DBzU1fzProxr8gANhw4e9MpaT9r9LuaPpf2nfoT5kqVVB1vyolFxgAyM7kCXCzadalpy5nZWnDvNnKvHBB/vVvVdexL1hVqJw/dVKGf+yLEtrjQclg0NYvF+nC6b9UvrK3ajUPU6sHH7b08a93q6JGjtWWLxZqx9Ilquznr9b9Bym4ddsb+n4AypYfl61U5SpV1P+p4fLx99PhA7/qPzFPKv3PVEmSj7+f/KsHWY2Z/vUXlj/f1rSJ2vfoquNHj2lw+06SJE+jUYNGj1RgzRr6+8JFbV23Xm/863ldOHfuxr0YABTCnt+OaeB/Zlvu4z9aJkm6v32o/jvyoZIKC3ANbBjrMIO5MJuUSPLx8dGBAwfk6+tbqIlr1aql9evXq3bt2oXqPzlpV6H6AUBpkTgw7/2hAKC0Wv7VhJIOAQCcq8mDJR3BDRExyXn/3970n2ZOm6s0KXTlyZkzZ7R8+XKZbGzE+U+nTp1STk5OkQMDAAAAAAC4Gdi1bGfQoEHFFQcAAAAAACgObBjrsEInTwpzig4AAAAAAICrKfKGsatXr9bq1auVnp5ulVgxGAyaM2eOU4IDAAAAAAAoaUVKnkycOFGTJk1Sy5YtFRQUZHWSCQAAAAAAuHkYWLbjsCIlT95//33NnTtX0dHRzo4HAAAAAAA4kcGtUIfswoYi5Z+ysrLUunVrZ8cCAAAAAABw0ylS8mTYsGFauHChs2MBAAAAAABOZnBz3lVWFWnZzqVLlzRz5kytWrVKt99+uzw8PKyeT5482SnBAQAAAAAAx7i5l3QEpV+Rkie7d+9W8+bNJUl79uyxesbmsQAAAAAAwJUUKXmyZs0aZ8cBAAAAAACKgVsZXm7jLEVKngAAAAAAgNKB03YcR/4JAAAAAADABipPAAAAAABwYSzbcRzJEwAAAAAAXBjJE8fxEQIAAAAAANhA5QkAAAAAAC6MyhPHkTwBAAAAAMCFkTxxHB8hAAAAAACADVSeAAAAAADgwqg8cRzJEwAAAAAAXJi7m7mkQyj1yD8BAAAAAADYQOUJAAAAAAAujGU7jiN5AgAAAACACyN54jg+QgAAAAAAABtIngAAAAAAANjAsh0AAAAAAFyYO2UTDuMjBAAAAAAAsIHKEwAAAAAAXJiboaQjKP1IngAAAAAA4MJYtuM4PkIAAAAAAAAbSJ4AAAAAAODC3NycdxWX06dPKzo6WiaTSSaTSdHR0Tpz5ozNMQaDIc/rjTfesPRp165druf9+vWzOz6W7QAAAAAA4MJKw7Kd/v376+jRo1qxYoUk6bHHHlN0dLS++eabfMekpqZa3S9fvlxDhw7Vgw8+aNUeExOjSZMmWe7Lly9vd3wkTwAAAAAAQInZt2+fVqxYoU2bNik8PFySNGvWLEVGRmr//v1q0KBBnuMCAwOt7pcuXar27durXr16Vu0VKlTI1ddepSD/BAAAAAAAisrdzXlXcUhKSpLJZLIkTiQpIiJCJpNJGzduLNQcx48f13fffaehQ4fmerZgwQL5+vqqSZMmGjt2rM6dO2d3jFSeAAAAAADgwpyZ9MjMzFRmZqZVm9FolNFoLPKcaWlp8vf3z9Xu7++vtLS0Qs3x8ccfq3LlynrggQes2h955BHVrVtXgYGB2rNnj+Li4rRr1y4lJibaFSOVJwAAAAAAoFDi4+Mtm7peu+Lj4/PsO2HChHw3db12bdu2TdLVzV+vZzab82zPy4cffqhHHnlEXl5eVu0xMTG67777FBISon79+umLL77QqlWrtGPHDrvem8oTAAAAAABcmDNPyYmLi1NsbKxVW35VJyNGjCjwZJs6depo9+7dOn78eK5nJ06cUEBAQIExrV+/Xvv379fixYsL7BsaGioPDw8dPHhQoaGhBfa/huQJAAAAAAAuzL1wxRuFYs8SHV9fX/n6+hbYLzIyUhkZGdqyZYtatWolSdq8ebMyMjLUunXrAsfPmTNHYWFhatasWYF99+7dq+zsbAUFBRX8Av/Ash0AAAAAAFBiGjVqpE6dOikmJkabNm3Spk2bFBMTo27dulmdtNOwYUMlJCRYjT179qw+//xzDRs2LNe8v/32myZNmqRt27bp8OHDWrZsmXr37q0WLVqoTZs2dsVI8gQAAAAAABd2s5+2I109Eadp06aKiopSVFSUbr/9dn3yySdWffbv36+MjAyrtkWLFslsNuvhhx/ONaenp6dWr16tjh07qkGDBho1apSioqK0atUqubu72xUfy3YAAAAAAHBhxZn0cBYfHx/Nnz/fZh+z2Zyr7bHHHtNjjz2WZ/+aNWtq3bp1TomvFHyEAAAAAAAAJYfKEwAAAAAAXFg5NyfuGFtGkTwBAAAAAMCFlYZlOzc7PkIAAAAAAAAbqDwBAAAAAMCFubNqx2EkTwAAAAAAcGEs23EcHyEAAAAAAIANJE8AAAAAAABsYNkOAAAAAAAujGU7juMjBAAAAAAAsIHKEwAAAAAAXJi7G8ftOIrkCQAAAAAALoxlO47jIwQAAAAAALCByhMAAAAAAFyYO6t2HEbyBAAAAAAAF8aeJ45j2Q4AAAAAAIANVJ4AAAAAAODC2DDWcQaz2Wwu6SCAGyUzM1Px8fGKi4uT0Wgs6XAAwGF8rwFwNXyvAbgZkTxBmXL27FmZTCZlZGTI29u7pMMBAIfxvQbA1fC9BuBmRPEOAAAAAACADSRPAAAAAAAAbCB5AgAAAAAAYAPJE5QpRqNR48ePZ/MxAC6D7zUArobvNQA3IzaMBQAAAAAAsIHKEwAAAAAAABtIngAAAAAAANhA8gQAAAAAAMAGkidwaXXq1JHBYJDBYNCZM2cKPW7u3LmWcc8880yxxQcA9uJ7DYCrufbdVKVKFbvGTZgwwTJ2ypQpxRIbAFxD8gQub9KkSUpNTZXJZJIkXbp0SYMHD1bTpk1Vrlw59erVK9eYvn37KjU1VZGRkTc4WgAo2PXfa2vXrlXPnj0VFBSkihUrqnnz5lqwYIHVGL7XANzMPvroIx04cMByn5qaqv79+6tBgwZyc3PLM+k7duxYpaamqkaNGjcwUgBlFckTuLzKlSsrMDBQBoNBkpSTk6Py5ctr1KhRuu+++/IcU758eQUGBsrT0/NGhgoAhXL999rGjRt1++23a8mSJdq9e7eGDBmigQMH6ptvvrGM4XsNwM2sSpUq8vf3t9xnZmbKz89PL7zwgpo1a5bnmEqVKikwMFDu7u43KkwAZVi5kg4AcES7du0UEhIiSZo/f77c3d31xBNP6KWXXrL8p+J6FStW1IwZMyRJ//d//2dX2TsAFLeifK89//zzVvejRo3SypUrlZCQoO7duxd7zABgS7t27XT77bfLy8tLs2fPlqenp4YPH64JEybkO6ZOnTp65513JEkffvjhDYoUAPJH5QlKvY8//ljlypXT5s2bNXXqVL399tuaPXt2SYcFAEXmjO+1jIwM+fj4FFOEAGCfjz/+WBUrVtTmzZv1+uuva9KkSUpMTCzpsACg0Kg8QalXs2ZNvf322zIYDGrQoIF++uknvf3224qJiSnp0ACgSBz9Xvviiy+0detWffDBB8UcKQAUzu23367x48dLkoKDgzVt2jStXr1aHTp0KOHIAKBwqDxBqRcREWFVyh4ZGamDBw8qJyenBKMCgKJz5Htt7dq1Gjx4sGbNmqUmTZoUZ5gAUGi333671X1QUJDS09NLKBoAsB/JEwAAXMS6devUvXt3TZ48WQMHDizpcADAwsPDw+reYDDoypUrJRQNANiP5AlKvU2bNuW6Dw4OZud1AKVWUb7X1q5dq65du+q///2vHnvsseIOEQAAoExhzxOUekeOHFFsbKwef/xx7dixQ++++67eeustm2N+/vlnZWVl6a+//tK5c+eUnJwsSWrevHnxBwwABbD3e+1a4uTpp5/Wgw8+qLS0NEmSp6cnm8YCKLWu/fvs/PnzOnHihJKTk+Xp6anGjRuXbGAAyiSSJyj1Bg4cqL///lutWrWSu7u7Ro4cWeBvXbt06aI//vjDct+iRQtJktlsLtZYAaAw7P1emzt3ri5evKj4+HjFx8db2u+++26tXbv2BkQMAM537d9nkrR9+3YtXLhQtWvX1uHDh0suKABlFskTlHoeHh6aMmWKZsyYUegx/KUL4GZm7/fa3LlzNXfu3OINCgCKKK8k7ldffVXgOH6pBeBmwp4ncHnPPfecKlWqpIyMjEKPWbBggSpVqqT169cXY2QAUDR8rwFwNQ8//LBq1Khh15hXX31VlSpVUkpKSjFFBQD/Q+UJXNq6deuUnZ0tSapcuXKhx/Xo0UPh4eGSpCpVqhRHaABQJHyvAXA1Bw8elCS7N/sfPny4+vTpI0ny8/NzelwA8E8GM/VwAAAAAAAA+WLZDgAAAAAAgA0kTwAAAAAAAGwgeQIAAAAAAGADyRMAAAAAAAAbSJ4AAAAAAADYQPIEAAAAAADABpInAAAAAAAANpA8AQAAAAAAsIHkCQAAAAAAgA3/D8I+NT6uHehsAAAAAElFTkSuQmCC", - "text/plain": [ - "
" + "cell_type": "markdown", + "id": "d8c28cfd", + "metadata": { + "id": "d8c28cfd" + }, + "source": [ + "## Start with some data analysis\n", + "\n", + "This data contains prices and sales for two of our competing products and was generated using another script, which can be found [here](). Let's load the data and take a quick look." ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(15, 5))\n", - "sns.heatmap(df[['p[1]','p[2]','n[1]']].corr(),annot=True, center=0,ax=axes)\n", - "\n", - "axes.set_title('Correlations')\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "a4aa96d4", - "metadata": {}, - "source": [ - "In this problem, we've assumed that the amount of space we have available for the products is 200 units. In retail, this could be the amount of warehouse space, or for ticketing this could represent the number of seats available." - ] - }, - { - "cell_type": "markdown", - "id": "6e9c6157", - "metadata": { - "id": "6e9c6157" - }, - "source": [ - "### Building regressors to predict sales\n", - "\n", - "The prices for each category item will be used to predict the number of Category 1 items sold. Here we build a regression model to form this relationship which will later be used as part of the optimization model. " - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "a70d5f2c", - "metadata": { - "id": "a70d5f2c" - }, - "outputs": [], - "source": [ - "from sklearn.compose import make_column_transformer\n", - "from sklearn.linear_model import LinearRegression\n", - "from sklearn.pipeline import make_pipeline\n", - "from sklearn.metrics import r2_score\n", - "from sklearn.model_selection import train_test_split #importing scikit-learn's function for data splitting\n", - "from sklearn.ensemble import GradientBoostingRegressor #importing scikit-learn's gradient booster regressor function\n", - "from sklearn.model_selection import cross_validate #improting scikit-learn's cross validation function" - ] - }, - { - "cell_type": "code", - "execution_count": 171, - "id": "3095da93", - "metadata": { - "colab": { - "base_uri": "https://localhost:8080/" }, - "id": "3095da93", - "outputId": "e6415f1b-34b8-4c25-ebf5-f4d1fa2223e1" - }, - "outputs": [], - "source": [ - "X = df[[\"p[1]\",\"p[2]\"]]\n", - "y = df[\"n[1]\"]\n", - "# Split the data for training and testing\n", - "X_train, X_test, y_train, y_test = train_test_split(\n", - " X, y, train_size=0.75, random_state=1\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "f8d713c2", - "metadata": {}, - "source": [ - "First we'll start with a linear regression model. " - ] - }, - { - "cell_type": "code", - "execution_count": 172, - "id": "33f4eaec", - "metadata": { - "colab": { - "base_uri": "https://localhost:8080/", - "height": 92 + { + "cell_type": "code", + "execution_count": 3, + "id": "eea57a4c", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 424 + }, + "id": "eea57a4c", + "outputId": "28a70699-d016-4f52-96e3-bc9badf6e0b0" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " p[1] p[2] n[1]\n", + "0 356.12 197.67 108.0\n", + "1 358.05 189.68 66.0\n", + "2 340.79 260.35 130.0\n", + "3 353.76 133.53 55.0\n", + "4 341.37 229.80 91.0\n", + ".. ... ... ...\n", + "995 357.63 241.54 68.0\n", + "996 352.58 212.95 87.0\n", + "997 355.28 189.50 94.0\n", + "998 369.75 166.33 51.0\n", + "999 349.31 222.07 114.0\n", + "\n", + "[1000 rows x 3 columns]" + ], + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
p[1]p[2]n[1]
0356.12197.67108.0
1358.05189.6866.0
2340.79260.35130.0
3353.76133.5355.0
4341.37229.8091.0
............
995357.63241.5468.0
996352.58212.9587.0
997355.28189.5094.0
998369.75166.3351.0
999349.31222.07114.0
\n", + "

1000 rows × 3 columns

\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n", + "\n", + "
\n", + " \n", + "\n", + "\n", + "\n", + " \n", + "
\n", + "\n", + "
\n", + " \n", + " \n", + " \n", + "
\n", + "\n", + "
\n", + "
\n" + ], + "application/vnd.google.colaboratory.intrinsic+json": { + "type": "dataframe", + "variable_name": "df", + "summary": "{\n \"name\": \"df\",\n \"rows\": 1000,\n \"fields\": [\n {\n \"column\": \"p[1]\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 14.409836071934338,\n \"min\": 300.0,\n \"max\": 400.0,\n \"num_unique_values\": 911,\n \"samples\": [\n 350.69,\n 359.04,\n 351.26\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"p[2]\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 29.093624863457553,\n \"min\": 100.0,\n \"max\": 300.0,\n \"num_unique_values\": 947,\n \"samples\": [\n 204.43,\n 211.29,\n 175.08\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"n[1]\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 29.65004858986544,\n \"min\": 0.0,\n \"max\": 200.0,\n \"num_unique_values\": 148,\n \"samples\": [\n 70.0,\n 105.0,\n 48.0\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n }\n ]\n}" + } + }, + "metadata": {}, + "execution_count": 3 + } + ], + "source": [ + "df = pd.read_csv('https://raw.githubusercontent.com/Gurobi/modeling-examples/master/pricing_competing_products/price_value_data.csv')\n", + "df" + ] }, - "id": "33f4eaec", - "outputId": "d696a20d-da2a-4f48-8224-9501c4328615" - }, - "outputs": [ { - "data": { - "text/plain": [ - "(array([0.80665368, 0.81004987, 0.7941077 , 0.79780172, 0.79495142]),\n", - " array([0.77650521, 0.74227784, 0.82105606, 0.81102818, 0.82166193]))" + "cell_type": "markdown", + "id": "b65707ac", + "metadata": { + "id": "b65707ac" + }, + "source": [ + "### What's in the data?\n", + "The data contains three columns:\n", + "1. `p[1]` is the price (in dollars) of the first category (let's call it Category 1).\n", + "2. `p[2]` is the price (in dollars) of the second category (Category 2).\n", + "3. `n[1]` is the number of the items sold that are of Category 1.\n", + "\n", + "We don't see a column for `n[2]`, which would be the number of items sold that are Category 2. Here is where we make a **pretty big assumption** that we will sell all of the items. This makes our decision to be how to divvy up the limited space we have in order to maximize our revenue.\n", + "The data was created to have a couple of key characteristics.\n", + "1. As the price of Category 1 goes up, the number sold should decrease, so `p[1]` and `n[1]` have a negative correlation.\n", + "2. As the price of Category 2 goes up, the number sold of Category 1 should increase, so `p[2]` and `n[1]` have a positive correlation.\n", + "\n", + "The correlation plot of the columns of the data is below." ] - }, - "execution_count": 172, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "linear_regressor = make_pipeline(LinearRegression())\n", - "linear_regressor.fit(X_train, y_train)\n", - "linear_regression_validation = cross_validate(linear_regressor, X_train, y_train, cv=5, return_train_score=True, return_estimator=True)\n", - "\n", - "linear_regression_validation['train_score'],linear_regression_validation['test_score']" - ] - }, - { - "cell_type": "markdown", - "id": "cc9dd59b", - "metadata": {}, - "source": [ - "Let's try a gradient boosting model as well. " - ] - }, - { - "cell_type": "code", - "execution_count": 173, - "id": "9b8eb814", - "metadata": { - "colab": { - "base_uri": "https://localhost:8080/" }, - "id": "9b8eb814", - "outputId": "de999b5f-807f-4ff6-8356-bbdb9f88db51" - }, - "outputs": [ { - "data": { - "text/plain": [ - "(array([0.83316977, 0.83568563, 0.82395805, 0.82872996, 0.82289456]),\n", - " array([0.75951446, 0.75090466, 0.79206779, 0.79360144, 0.80517282]))" + "cell_type": "code", + "execution_count": 4, + "id": "4cb22fe7", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 468 + }, + "id": "4cb22fe7", + "outputId": "37183ae6-b780-4225-cc3d-bbaa2f4c7924" + }, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": [ + "
" + ], + "image/png": "iVBORw0KGgoAAAANSUhEUgAABE8AAAHDCAYAAADV6DQTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABSKklEQVR4nO3de1zUVf7H8feAOoDITW7iBbyUaaYopmK1WpJ4S21t09JQ8lpZJlaGmXjZYk0z7WrmvdWfdrVNyzTTrqSmS95NTSVLQFFABQeE+f3RNruTMDLMIDK8no/H97F9z/ecM5/vPPYx5cfzOcdgNpvNAgAAAAAAQIncKjsAAAAAAACAaxnJEwAAAAAAABtIngAAAAAAANhA8gQAAAAAAMAGkicAAAAAAAA2kDwBAAAAAACwgeQJAAAAAACADSRPAAAAAAAAbCB5AgAAAAAAYAPJEwAAqpClS5fKYDDo2LFjTpvz2LFjMhgMWrp0qdPmBAAAcCUkTwAAkHTkyBGNHj1aTZo0kYeHh3x8fHTLLbdo3rx5ys/Pr+zwnGLlypWaO3duZYcBAABQ5dSo7AAAAKhs69at09/+9jcZjUbFxcWpVatWKigo0DfffKMnn3xSe/fu1YIFCyo7TIetXLlSe/bs0eOPP27VHh4ervz8fNWsWbNyAgMAALjGkTwBAFRrR48e1aBBgxQeHq4vvvhC9erVszx75JFHdPjwYa1bt86hzzCbzbp48aI8PT0ve3bx4kXVqlVLbm6VtxjUYDDIw8Oj0j4fAADgWkfZDgCgWnvhhRd0/vx5LVq0yCpx8odmzZpp3LhxkqRLly5pxowZatq0qYxGoyIiIjRp0iSZTCarMREREerTp48+++wztW/fXp6ennrzzTe1ZcsWGQwGrVq1SpMnT1b9+vXl5eWl3NxcSdLWrVvVo0cP+fr6ysvLS126dNG33357xXf46KOP1Lt3b4WFhcloNKpp06aaMWOGioqKLH26du2qdevW6fjx4zIYDDIYDIqIiJBU+p4nX3zxhW677TbVrl1bfn5+6tevn/bv32/VZ+rUqTIYDDp8+LCGDRsmPz8/+fr6Kj4+Xnl5eVZ9N27cqFtvvVV+fn7y9vZW8+bNNWnSpCu+HwAAQGVj5QkAoFr7+OOP1aRJE3Xu3PmKfUeMGKFly5bpnnvu0YQJE7R161YlJydr//79+vDDD636Hjx4UPfdd59Gjx6tkSNHqnnz5pZnM2bMUK1atfTEE0/IZDKpVq1a+uKLL9SzZ09FRUUpKSlJbm5uWrJkie644w59/fXX6tChQ6lxLV26VN7e3kpISJC3t7e++OILTZkyRbm5uZo1a5Yk6ZlnnlFOTo5OnDihl156SZLk7e1d6pyff/65evbsqSZNmmjq1KnKz8/XK6+8oltuuUU7d+60JF7+cO+996px48ZKTk7Wzp07tXDhQgUHB2vmzJmSpL1796pPnz5q3bq1pk+fLqPRqMOHD5cpOQQAAFDpzAAAVFM5OTlmSeZ+/fpdsW9qaqpZknnEiBFW7U888YRZkvmLL76wtIWHh5slmdevX2/Vd/PmzWZJ5iZNmpjz8vIs7cXFxebrrrvOHBsbay4uLra05+XlmRs3bmy+8847LW1LliwxSzIfPXrUqt+fjR492uzl5WW+ePGipa13797m8PDwy/oePXrULMm8ZMkSS1tkZKQ5ODjYnJWVZWn78ccfzW5ubua4uDhLW1JSklmS+cEHH7Sa8+677zbXrVvXcv/SSy+ZJZlPnTp12ecDAABc6yjbAQBUW3+Uy9SpU+eKfT/55BNJUkJCglX7hAkTJOmyfVEaN26s2NjYEucaOnSo1f4nqampOnTokO6//35lZWXp9OnTOn36tC5cuKBu3brpq6++UnFxcamx/e9c586d0+nTp3XbbbcpLy9PBw4cuOK7/dnJkyeVmpqqYcOGKSAgwNLeunVr3XnnnZbv4n+NGTPG6v62225TVlaW5Tv28/OT9HuJka13AQAAuBaRPAEAVFs+Pj6Sfk84XMnx48fl5uamZs2aWbWHhobKz89Px48ft2pv3LhxqXP9+dmhQ4ck/Z5UCQoKsroWLlwok8mknJycUufbu3ev7r77bvn6+srHx0dBQUEaMmSIJNkcV5o/3uV/S43+0KJFC0ti5381atTI6t7f31+SdPbsWUnSwIEDdcstt2jEiBEKCQnRoEGD9M4775BIAQAAVQJ7ngAAqi0fHx+FhYVpz549ZR5jMBjK1K+kk3VKe/ZHAmHWrFmKjIwscUxp+5NkZ2erS5cu8vHx0fTp09W0aVN5eHho586dmjhx4lVLTri7u5fYbjabJf3+zl999ZU2b96sdevWaf369Vq9erXuuOMObdiwodTxAAAA1wKSJwCAaq1Pnz5asGCBUlJSFB0dXWq/8PBwFRcX69ChQ2rRooWlPSMjQ9nZ2QoPDy93DE2bNpX0ezInJibGrrFbtmxRVlaWPvjgA/3lL3+xtB89evSyvmVN/PzxLgcPHrzs2YEDBxQYGKjatWvbFackubm5qVu3burWrZvmzJmj559/Xs8884w2b95s93sDAABcTZTtAACqtaeeekq1a9fWiBEjlJGRcdnzI0eOaN68eerVq5ckae7cuVbP58yZI0nq3bt3uWOIiopS06ZNNXv2bJ0/f/6y56dOnSp17B8rNv5Y4SFJBQUFev311y/rW7t27TKV8dSrV0+RkZFatmyZsrOzLe179uzRhg0bLN+FPc6cOXNZ2x+rbP581DMAAMC1hpUnAIBqrWnTplq5cqUGDhyoFi1aKC4uTq1atVJBQYG+++47vfvuuxo2bJjGjRunoUOHasGCBZZSmW3btmnZsmXq37+/br/99nLH4ObmpoULF6pnz5668cYbFR8fr/r16+vXX3/V5s2b5ePjo48//rjEsZ07d5a/v7+GDh2qxx57TAaDQW+//bZVMuUPUVFRWr16tRISEnTzzTfL29tbd911V4nzzpo1Sz179lR0dLSGDx9uOarY19dXU6dOtfsdp0+frq+++kq9e/dWeHi4MjMz9frrr6tBgwa69dZb7Z4PAADgaiJ5AgCo9vr27atdu3Zp1qxZ+uijj/TGG2/IaDSqdevWevHFFzVy5EhJ0sKFC9WkSRMtXbpUH374oUJDQ5WYmKikpCSHY+jatatSUlI0Y8YMvfrqqzp//rxCQ0PVsWNHjR49utRxdevW1dq1azVhwgRNnjxZ/v7+GjJkiLp163bZaT8PP/ywUlNTtWTJEr300ksKDw8vNXkSExOj9evXKykpSVOmTFHNmjXVpUsXzZw50+ZmuKXp27evjh07psWLF+v06dMKDAxUly5dNG3aNPn6+to9HwAAwNVkMJf0V1MAAAAAAACQxJ4nAAAAAAAANpE8AQAAAAAAsIHkCQAAAAAAgA0kTwAAAAAAgNN99dVXuuuuuxQWFiaDwaA1a9ZcccyWLVvUrl07GY1GNWvWTEuXLr2sz2uvvaaIiAh5eHioY8eO2rZtm/OD/xOSJwAAAAAAwOkuXLigNm3a6LXXXitT/6NHj6p37966/fbblZqaqscff1wjRozQZ599ZumzevVqJSQkKCkpSTt37lSbNm0UGxurzMzMinoNSZy2AwAAAAAAKpjBYNCHH36o/v37l9pn4sSJWrdunfbs2WNpGzRokLKzs7V+/XpJUseOHXXzzTfr1VdflSQVFxerYcOGevTRR/X0009XWPysPAEAAAAAAGViMpmUm5trdZlMJqfMnZKSopiYGKu22NhYpaSkSJIKCgq0Y8cOqz5ubm6KiYmx9KkoNSp0dnvsfb+yIwAAp4ro9URlhwAATvXYqjWVHQIAOFVCdJvKDuHqcOKft5Pf3a1p06ZZtSUlJWnq1KkOz52enq6QkBCrtpCQEOXm5io/P19nz55VUVFRiX0OHDjg8Ofbcu0kTwAAAAAAwDUtMTFRCQkJVm1Go7GSorl6SJ4AAAAAAIAyMRqNFZYsCQ0NVUZGhlVbRkaGfHx85OnpKXd3d7m7u5fYJzQ0tEJi+gN7ngAAAAAA4MLMRUVOuypSdHS0Nm3aZNW2ceNGRUdHS5Jq1aqlqKgoqz7FxcXatGmTpU9FIXkCAAAAAACc7vz580pNTVVqaqqk348iTk1NVVpamqTfS4Di4uIs/ceMGaOff/5ZTz31lA4cOKDXX39d77zzjsaPH2/pk5CQoLfeekvLli3T/v379dBDD+nChQuKj4+v0HehbAcAAAAAAFdWdKlSPvaHH37Q7bffbrn/Y6+UoUOHaunSpTp58qQlkSJJjRs31rp16zR+/HjNmzdPDRo00MKFCxUbG2vpM3DgQJ06dUpTpkxRenq6IiMjtX79+ss2kXU2g9lsNlfoJ5QVp+0AcDGctgPA1XDaDgBXU11O2ynesdxpc7lFxV25kwuibAcAAAAAAMAGynYAAAAAAHBlFbzRa3VA8gQAAAAAABdmrqQ9T1wJZTsAAAAAAAA2sPIEAAAAAABXxsoTh5E8AQAAAADAhZmLSZ44irIdAAAAAAAAG1h5AgAAAACAK+O0HYeRPAEAAAAAwIVx2o7jKNsBAAAAAACwgZUnAAAAAAC4MlaeOIzkCQAAAAAALsxczJ4njqJsBwAAAAAAwAZWngAAAAAA4MLYMNZxJE8AAAAAAHBlJE8cRtkOAAAAAACADaw8AQAAAADAhbFhrONYeQIAAAAAAGADyRMAAAAAAAAbKNsBAAAAAMCVsWGsw0ieAAAAAADgwjiq2HGU7QAAAAAAANjAyhMAAAAAAFwZK08cRvIEAAAAAAAXxlHFjqNsBwAAAAAAwAZWngAAAAAA4Moo23EYyRMAAAAAAFyYuYiyHUdRtgMAAAAAAGADK08AAAAAAHBhZsp2HEbyBAAAAAAAV1ZM8sRRlO0AAAAAAADYwMoTAAAAAABcGBvGOo7kCQAAAAAArozkicMo2wEAAAAAALCBlScAAAAAALgwTttxHMkTAAAAAABcGWU7DqNsBwAAAAAAwAZWngAAAAAA4MI4bcdxJE8AAAAAAHBh5mKSJ46ibAcAAAAAAMAGkicAAAAAALiyoiLnXeXw2muvKSIiQh4eHurYsaO2bdtWat+uXbvKYDBcdvXu3dvSZ9iwYZc979GjR7liKyvKdgAAAAAAQIVYvXq1EhISNH/+fHXs2FFz585VbGysDh48qODg4Mv6f/DBByooKLDcZ2VlqU2bNvrb3/5m1a9Hjx5asmSJ5d5oNFbcS4iVJwAAAAAAoILMmTNHI0eOVHx8vFq2bKn58+fLy8tLixcvLrF/QECAQkNDLdfGjRvl5eV1WfLEaDRa9fP396/Q9yB5AgAAAACACzMXFTntskdBQYF27NihmJgYS5ubm5tiYmKUkpJSpjkWLVqkQYMGqXbt2lbtW7ZsUXBwsJo3b66HHnpIWVlZdsVmL8p2AAAAAABwYeaiYqfNZTKZZDKZrNqMRmOJZTOnT59WUVGRQkJCrNpDQkJ04MCBK37Wtm3btGfPHi1atMiqvUePHvrrX/+qxo0b68iRI5o0aZJ69uyplJQUubu7l+OtroyVJwAAAAAAoEySk5Pl6+trdSUnJ1fIZy1atEg33XSTOnToYNU+aNAg9e3bVzfddJP69++vtWvXavv27dqyZUuFxCGRPAEAAAAAwLUVFTvtSkxMVE5OjtWVmJhY4scGBgbK3d1dGRkZVu0ZGRkKDQ21GfKFCxe0atUqDR8+/Iqv16RJEwUGBurw4cNl/07sVOaynYSEBLsnnzx5sgICAuweBwAAAAAAnMPevUpsKa1EpyS1atVSVFSUNm3apP79+0uSiouLtWnTJo0dO9bm2HfffVcmk0lDhgy54uecOHFCWVlZqlevXpniKo8yJ0/mzp2r6Oho1apVq0z9v/nmG40dO5bkCQAAAAAA1VRCQoKGDh2q9u3bq0OHDpo7d64uXLig+Ph4SVJcXJzq169/WenPokWL1L9/f9WtW9eq/fz585o2bZoGDBig0NBQHTlyRE899ZSaNWum2NjYCnsPuzaM/fDDD0s8h7kkderUKVdAAAAAAADAecxF5kr77IEDB+rUqVOaMmWK0tPTFRkZqfXr11s2kU1LS5Obm/WOIgcPHtQ333yjDRs2XDafu7u7du3apWXLlik7O1thYWHq3r27ZsyYUeYVMeVR5uTJkiVL5OvrW+aJ33zzzct21AUAAAAAAFeXM0/bKY+xY8eWWqZT0iavzZs3l9lccsLH09NTn332mTPDK5MyJ0+GDh1q18T333+/3cEAFWX73qNa9NHX2nPkV506e06vTRyimI4tKzssACjV+ITxuu++QfLx8dEPP/ygyc88q2PHjtkc80DcAxo9apSCgoK0f/9+JSVN1Y8//ihJ8vX11fiE8brttttUv36YsrKytGHDRs15cY7OnTt3Fd4IQHVlNpv1w4fv6MCXm2TKu6DQ627QbXEj5Bta+t4ExcXF2vHhOzqU8rXycrJV2y9A19/aRe36DpDBYLD0O/vbCW19Z4VOHtyn4qJi+ddvoDvHTlCduoFX49UAVCOctoNqIc9UoOYRoUoa2beyQwGAKxozZrTihw3TM5Mmq3+/u5Wfl6/lby+T0Vj6vmN9+vTW5MnPaN68eerdp4/27d+v5W8vs9QJh4SEKCQkWM8/97y63xmrJ554Ul26dNHMF2ZerdcCUE39+MlH2rPxU902dKTunvK8ahiNWvfic7pUUFDqmNR1a7Rv80bdMmS4Bj7/kjreO1g/fvov7fn8U0ufnMx0ffTcFPnVq6+7np6qe/4+S+36DlCNmjWvxmsBVYq5qNhpV3Xl1OTJjz/+KHd3d2dOCThFl3bNNf7+7rqz042VHQoAXNGDwx/UK6++qo0bN+rAgQNKSJigkOAQde/evdQxI0aM0KpVq/Xuu+/p8KHDembSM8rPz9e99/5NkvTTTz/poTEPa9OmTUpLS1PKdymaPWu2unW7g393A6gwZrNZuzd8onZ9/6qIdjerbsNw3T5yrPLOntWxndtLHZdx+CeFt22v8Mh2qhMUrCY3d1KDG1sr8+f/HkO6/b1VatS6rToNHKLA8MbyDQ5VRNv28vQp+1YDQHVhLjY77aqunL7ypLS6JAAAcGUNGzZUcHCwvv3mG0vbuXPnlJqaqnbt2pU4pmbNmmp1UyurMWazWd9+822pYySpjk8dnT9/XkVOPL4QAP7XuVOZysvJVv2WrS1tRi8vBTdtpowjP5U6LqTZ9fp13x5lp/8mScpKO6b0QwfV6Ka2kiRzcbHSdu2Ub2g9rZv9nJY9OkIfTp+kozu2VewLAai27Dpt569//avN5zk5OVY1iAAAwD5BwUGSpFOnT1u1nzp9WkFBQSWO8ff3V40aNXS6hDFNmzYtdcyjjz6q//u/VU6IGgBKlpeTLUny/NPBE54+vpZnJWnbu78K8/O1OnG83NzcVFxcrA4DBum6zrdJkvJzc1V48aJS132kmwcMVMe/DdYvu1O14dUXddfEJIXdwN52wP+qzNN2XIVdyZOPP/5Yd955Z6mn6JT1b65MJpNMJpNVm7GgUMZa1CcCAKqXfv376fnnn7PcPxg/vMI/09vbW0uWLNbhw4c096W5Ff55AKqPQ999ra+WLbDc9xyfWK55jmxL0aHvv1G30Y/Jv35DZaUd03crl8rLz1/Nb+0qs/n3fRci2rVX69g+kqTA8AhlHD6ofZs3kDwB/sTMIlOH2ZU8adGihQYMGKDhw0v+D7vU1FStXbv2ivMkJydr2rRpVm1JD/1NUx8ZaE84AABUeZ9v/Fyp/0613Neq9fumsEGBgTqVecrSHhQYqH379pU4x9mzZ3Xp0iUFBlqfLhEUGKhTp05ZtdWuXVvLli/V+QvnNXrUaF26dMlJbwIAUnjb9rqn6XWW+6JLhZKk/Jwc1fbzt7Tn5+aobqOIUuf5/p1/KrJXPzXrdIskqW7DRjqfdUqpa9eo+a1d5VHHR27u7vIPa2A1zi+svtJ/OujENwKA39m150lUVJR27txZ6nOj0ahGjRpdcZ7ExETl5ORYXYkjbZcEAQDgii5cuKDjx49brkOHDikzM1Odb7nF0sfb21uRkZGl/ju4sLBQe3bvsRpjMBjU+ZbOVmO8vb319j+Xq7CgUCOGj5TJVPpJFwBQHrU8PeUbEmq5/MMayMvXT7/u223pU5Cfp8wjhxXS9PpS57lkMsngZv1HFYObm2V/RfcaNRTUuKmyT/5m1Scn/aTqBHJMMfBn5iKz067qyq6VJ/Pnz7dZmtOiRQsdPXr0ivMYjUYZjUbrRkp2UIEu5JuUlp5luT+ReUb7j/4mX28vhQX5VV5gAFCCxYsW69FHx+rY0WP65ZdfNGFCgjIyM7RhwwZLnxUr/6nPPtug5cuWS5IWLlyoF198Ubt37VLqjz9q+IMPysvLS++++56k/yRO3l4uD09PPT5uvOrU8VadOt6SpKysMyourr5HDwKoOAaDQTd176WdH38g39B6qhMYrB8+WCUvf39FtLvZ0u/jmdPVOKqDWsX0kCSFR0bp3x9/IO+AQAXUb6DTace067O1an7b7ZYxbXr21eevv6R6zVsorEUr/bI7VcdTd+iup6de7dcErnn8a95xdiVPLkt4AFXEniO/Km7KQst98pJPJEl3395O/3j0nsoKCwBKNH/+m/L08lJy8vPy8fHR9h+2a2jcMKuVIuGNwhXg/98l8GvXrlNA3boan5CgoKBA7d+3X0Pjhlk2kW3V6ka1bff7KRVfff2l1efdesutOnHi16vwZgCqoza9+qnQZNJXS95UQV6eQq+/Qb0mTFKN/5QpSlJuZoYunsu13N8y5EFt/2C1vnl7ofJzc1TbL0Atut6pqH7//e+2xlEddNvQkfr3ujX6dsUS+YWGqfvYCap3/Q1X9f0AVA8GcxnPFs7NzZWPj0+ZJz537pzq1KlT9kj2vl/2vgBQBUT0eqKyQwAAp3ps1ZrKDgEAnCohuk1lh3BV/DKsvdPmarj0B6fNVZWUec8Tf39/ZWZmlnni+vXr6+effy5XUAAAAAAAwDnMRc67qqsyl+2YzWYtXLhQ3t7eZepfWFhY7qAAAAAAAACuFWVOnjRq1EhvvfVWmScODQ1VzZpsAgsAAAAAAKq2MidPjh07VmL7H1umGAwGpwQEAAAAAACch9N2HFfmPU/+bNGiRWrVqpU8PDzk4eGhVq1aaeHChVceCAAAAAAAUIXYdVTxH6ZMmaI5c+bo0UcfVXR0tCQpJSVF48ePV1pamqZPn+7UIAEAAAAAQPlU541enaVcyZM33nhDb731lu677z5LW9++fdW6dWs9+uijJE8AAAAAALhGFBezzYajylW2U1hYqPbtLz8nOioqSpcuXXI4KAAAAAAAgGtFuZInDzzwgN54443L2hcsWKDBgwc7HBQAAAAAAHCO4mLnXdVVucp2pN83jN2wYYM6deokSdq6davS0tIUFxenhIQES785c+Y4HiUAAAAAACgX9jxxXLmSJ3v27FG7du0kSUeOHJEkBQYGKjAwUHv27LH04/hiAAAAAABQ1ZUrebJ582ZnxwEAAAAAACoAG8Y6rtxlOwAAAAAA4NpXTNmOw8q1YSwAAAAAAEB1wcoTAAAAAABcGGU7jiN5AgAAAACACzOTPHEYZTsAAAAAAAA2sPIEAAAAAAAXVlxc2RFUfaw8AQAAAAAAsIGVJwAAAAAAuDA2jHUcyRMAAAAAAFwYyRPHUbYDAAAAAABgAytPAAAAAABwYUWsPHEYyRMAAAAAAFwYZTuOo2wHAAAAAADABlaeAAAAAADgworNrDxxFCtPAAAAAAAAbCB5AgAAAAAAYANlOwAAAAAAuLDi4sqOoOojeQIAAAAAgAsrYs8Th1G2AwAAAAAAKsxrr72miIgIeXh4qGPHjtq2bVupfZcuXSqDwWB1eXh4WPUxm82aMmWK6tWrJ09PT8XExOjQoUMV+g4kTwAAAAAAcGHFxQanXfZavXq1EhISlJSUpJ07d6pNmzaKjY1VZmZmqWN8fHx08uRJy3X8+HGr5y+88IJefvllzZ8/X1u3blXt2rUVGxurixcv2h1fWZE8AQAAAADAhRWZDU677DVnzhyNHDlS8fHxatmypebPny8vLy8tXry41DEGg0GhoaGWKyQkxPLMbDZr7ty5mjx5svr166fWrVtr+fLl+u2337RmzZryfD1lQvIEAAAAAAA4XUFBgXbs2KGYmBhLm5ubm2JiYpSSklLquPPnzys8PFwNGzZUv379tHfvXsuzo0ePKj093WpOX19fdezY0eacjmLDWAAAAAAAXFixEzeMNZlMMplMVm1Go1FGo/GyvqdPn1ZRUZHVyhFJCgkJ0YEDB0qcv3nz5lq8eLFat26tnJwczZ49W507d9bevXvVoEEDpaenW+b485x/PKsIrDwBAAAAAMCFObNsJzk5Wb6+vlZXcnKy02KNjo5WXFycIiMj1aVLF33wwQcKCgrSm2++6bTPKA9WngAAAAAAgDJJTExUQkKCVVtJq04kKTAwUO7u7srIyLBqz8jIUGhoaJk+r2bNmmrbtq0OHz4sSZZxGRkZqlevntWckZGRZX0Nu7HyBAAAAAAAF1Zkdt5lNBrl4+NjdZWWPKlVq5aioqK0adMmS1txcbE2bdqk6OjossVeVKTdu3dbEiWNGzdWaGio1Zy5ubnaunVrmecsD1aeAAAAAADgwpy554m9EhISNHToULVv314dOnTQ3LlzdeHCBcXHx0uS4uLiVL9+fUvpz/Tp09WpUyc1a9ZM2dnZmjVrlo4fP64RI0ZI+v0knscff1x///vfdd1116lx48Z69tlnFRYWpv79+1fYe5A8AQAAAAAAFWLgwIE6deqUpkyZovT0dEVGRmr9+vWWDV/T0tLk5vbfopizZ89q5MiRSk9Pl7+/v6KiovTdd9+pZcuWlj5PPfWULly4oFGjRik7O1u33nqr1q9fLw8Pjwp7D4PZbDZX2Oz22Pt+ZUcAAE4V0euJyg4BAJzqsVVrKjsEAHCqhOg2lR3CVfFhu5grdyqju3d+7rS5qhJWngAAAAAA4MKKro0lE1UaG8YCAAAAAADYwMoTAAAAAABcWJEqb8NYV0HyBAAAAAAAF0bZjuMo2wEAAAAAALCBlScAAAAAALiwosoOwAWQPAEAAAAAwIWRPHEcZTsAAAAAAAA2kDwBAAAAAACwgbIdAAAAAABcGEcVO46VJwAAAAAAADaw8gQAAAAAABdWZDZXdghVHskTAAAAAABcGKftOI6yHQAAAAAAABtYeQIAAAAAgAtj5YnjSJ4AAAAAAODCSJ44jrIdAAAAAAAAG1h5AgAAAACACysSp+04iuQJAAAAAAAujLIdx10zyZOIXk9UdggA4FTHPpld2SEAgFP17D+kskMAAKdKOLS7skNAFXHNJE8AAAAAAIDzFZkp23EUyRMAAAAAAFwYZTuO47QdAAAAAAAAG1h5AgAAAACAC+O0HceRPAEAAAAAwIWRPHEcZTsAAAAAAAA2sPIEAAAAAAAXxoaxjiN5AgAAAACAC+OoYsdRtgMAAAAAAGADK08AAAAAAHBhbBjrOJInAAAAAAC4MJInjqNsBwAAAAAAwAaSJwAAAAAAADZQtgMAAAAAgAsr5rQdh7HyBAAAAAAAwAZWngAAAAAA4MLYMNZxJE8AAAAAAHBhJE8cR9kOAAAAAACADaw8AQAAAADAhRWxYazDSJ4AAAAAAODCKNtxHGU7AAAAAAAANpA8AQAAAADAhRWbzU67yuO1115TRESEPDw81LFjR23btq3Uvm+99ZZuu+02+fv7y9/fXzExMZf1HzZsmAwGg9XVo0ePcsVWViRPAAAAAABwYUUyO+2y1+rVq5WQkKCkpCTt3LlTbdq0UWxsrDIzM0vsv2XLFt13333avHmzUlJS1LBhQ3Xv3l2//vqrVb8ePXro5MmTluv//u//yvXdlBXJEwAAAAAAUCHmzJmjkSNHKj4+Xi1bttT8+fPl5eWlxYsXl9h/xYoVevjhhxUZGakbbrhBCxcuVHFxsTZt2mTVz2g0KjQ01HL5+/tX6HuQPAEAAAAAwIU5c+WJyWRSbm6u1WUymUr83IKCAu3YsUMxMTGWNjc3N8XExCglJaVMsefl5amwsFABAQFW7Vu2bFFwcLCaN2+uhx56SFlZWeX/gsqA5AkAAAAAAC7MmXueJCcny9fX1+pKTk4u8XNPnz6toqIihYSEWLWHhIQoPT29TLFPnDhRYWFhVgmYHj16aPny5dq0aZNmzpypL7/8Uj179lRRUVH5v6Qr4KhiAAAAAABQJomJiUpISLBqMxqNFfJZ//jHP7Rq1Spt2bJFHh4elvZBgwZZ/vmmm25S69at1bRpU23ZskXdunWrkFhIngAAAAAA4MLKs9FraYxGY5mTJYGBgXJ3d1dGRoZVe0ZGhkJDQ22OnT17tv7xj3/o888/V+vWrW32bdKkiQIDA3X48OEKS55QtgMAAAAAgAsrMpuddtmjVq1aioqKstrs9Y/NX6Ojo0sd98ILL2jGjBlav3692rdvf8XPOXHihLKyslSvXj274rMHyRMAAAAAAFAhEhIS9NZbb2nZsmXav3+/HnroIV24cEHx8fGSpLi4OCUmJlr6z5w5U88++6wWL16siIgIpaenKz09XefPn5cknT9/Xk8++aS+//57HTt2TJs2bVK/fv3UrFkzxcbGVth7ULYDAAAAAIALK3Zi2Y69Bg4cqFOnTmnKlClKT09XZGSk1q9fb9lENi0tTW5u/13X8cYbb6igoED33HOP1TxJSUmaOnWq3N3dtWvXLi1btkzZ2dkKCwtT9+7dNWPGjArbe0WSDGaznetuKkhEeOPKDgEAnOrYJ7MrOwQAcKqe/adWdggA4FSfHtpd2SFcFf2bt3XaXGsO/ttpc1UllO0AAAAAAADYQNkOAAAAAAAurPjaKDip0lh5AgAAAAAAYAPJEwAAAAAAABso2wEAAAAAwIUVVeJpO66C5AkAAAAAAC6s2Fxc2SFUeZTtAAAAAAAA2MDKEwAAAAAAXFgxZTsOI3kCAAAAAIALK+KoYodRtgMAAAAAAGADK08AAAAAAHBhlO04juQJAAAAAAAurJiyHYdRtgMAAAAAAGADK08AAAAAAHBhxZUdgAsgeQIAAAAAgAujbMdxlO0AAAAAAADYwMoTAAAAAABcGKftOI7kCQAAAAAALoyyHcdRtgMAAAAAAGADK08AAAAAAHBhlO04rszJk5dfftnuyePj41WnTh27xwEAAAAAAOcgeeK4MidPHn/8cTVo0EDu7u5l6v/LL7+oT58+JE8AAAAAAECVZlfZzg8//KDg4OAy9SVpAgAAAABA5Stm4YnDypw8SUpKkre3d5knnjRpkgICAsoVFAAAAAAAcA7KdhxnV/LEHomJiXYHAwAAAAAAcK3hqGK4hPEJ47Vt+1YdOLhf/1zxtiIiIq445oG4B/TNN1/r4MEDWrPmQ7Vp08byzNfXV1OnTdWmLzbpwMH9+va7b5Q0NYlyNADXjO17j2rM88t16/BkNf/rJH2+dV9lhwQAJeozeJCWbl6vj/b8oJfeW6HrW7cqtW+jZk31zKtztHTzen16aLf6DxtyWR/P2l4a/cxTWrrlM63ZvV0vrn5b1990Y0W+AlDlFcvstKu6sjt58sknn2jEiBF66qmndODAAatnZ8+e1R133OG04ICyGDNmtOKHDdMzkyarf7+7lZ+Xr+VvL5PRWKvUMX369Nbkyc9o3rx56t2nj/bt36/lby9T3bp1JUkhISEKCQnW8889r+53xuqJJ55Uly5dNPOFmVfrtQDApjxTgZpHhCppZN/KDgUASvWXXrEaNelJrXh1vh7tf6+O7v9Jf1/8pnxLKe/38PRQ+i8ntGT2XJ3JPFVin3HPTVPbW6I1+8lJeqj3X7Xzm+/0/LK3VDekbHszAkB52JU8Wblypfr27av09HSlpKSobdu2WrFiheV5QUGBvvzyS6cHCdjy4PAH9cqrr2rjxo06cOCAEhImKCQ4RN27dy91zIgRI7Rq1Wq9++57OnzosJ6Z9Izy8/N1771/kyT99NNPemjMw9q0aZPS0tKU8l2KZs+arW7d7ijziVMAUJG6tGuu8fd3152d+NtWANeuux+M06er39fG99co7fDPemXKdJny89X9nrtL7P/T7r1aNHOOvly3XoUFBZc9r2U06tbYGC16YY72bN+hk2m/aMUrb+i347+o9/0DK/p1AFRjdiVPZs2apTlz5mjt2rX6+uuvtWzZMo0ePVqLFi2qqPgAmxo2bKjg4GB9+803lrZz584pNTVV7dq1K3FMzZo11eqmVlZjzGazvv3m21LHSFIdnzo6f/68ioqKnPcCAAAALqpGzRq67saWSv3ue0ub2WxW6nffq0XbNjZGls69hrvca9RQock6sVJw8aJujGrrULyAKzObnXdVV3YdVXzo0CHdddddlvt7771XQUFB6tu3rwoLC3X33SVnkIGKEhQcJEk6dfq0Vfup06cVFBRU4hh/f3/VqFFDp0sY07Rp01LHPProo/q//1vlhKgBAABcn4+/v9xr1NDZ01lW7WezstSgaeNyzZl/IU/7dqbqvkdGK+3Iz8o+naUufXrphrZtdPJ4mjPCBlxSdd6rxFnsWnni4+OjjIwMq7bbb79da9eu1ZNPPqlXXnmlTPOYTCbl5uZaXebqnMJCmfXr30979+2xXDVr1Kzwz/T29taSJYt1+PAhzX1pboV/HgAAAEo3+8lEGQwGrfj2C/1r7w71i7tfX679VMX8eQJABbJr5UmHDh306aefqlOnTlbtXbp00ccff6w+ffqUaZ7k5GRNmzbNqs3Xx1d+fv72hINq6PONnyv136mW+1q1ft8UNigwUKf+Z1OxoMBA7dtX8skTZ8+e1aVLlxQYGGjVHhQYqFOnrDcmq127tpYtX6rzF85r9KjRunTpkpPeBAAAwLXlnj2rokuX5B9Y16rdv25dnT2VVcqoKzuZdkJPDY6X0dNTXt61dfbUaT09d5bSfznhaMiAyyK16Di7Vp6MHz9eHh4eJT7r2rWrPv74Y8XFxV1xnsTEROXk5Fhdvr5+9oSCaurChQs6fvy45Tp06JAyMzPV+ZZbLH28vb0VGRmpnTt3ljhHYWGh9uzeYzXGYDCo8y2drcZ4e3vr7X8uV2FBoUYMHymT6fJNywAAAFCyS4WXdGjvPkVGd7S0GQwGRXbupP3//tHh+U35+Tp76rS8fXwUdVtnff/5ZofnBFwVRxU7zq6VJ126dFGXLl1KfX777bfr9ttvv+I8RqNRRqPRqs1gMNgTCmCxeNFiPfroWB07eky//PKLJkxIUEZmhjZs2GDps2LlP/XZZxu0fNlySdLChQv14osvaveuXUr98UcNf/BBeXl56d1335P0n8TJ28vl4empx8eNV5063qpTx1uSlJV1RsXFxVf/RQHgf1zINykt/b9/c3si84z2H/1Nvt5eCgvyq7zAAOB/fLh4uSa88JwO7dmrg7t2q/+wB2T09NTG99dIkia88JyyMjK19MV5kn7fZLZRs6b/+eeaqhsSrCYtmiv/Qp5Opv0iSWp3a2cZDAadOHpMYeGNNHxigk78fFQb/jMnAFSEMidPcnNz5ePjU+aJz507pzp16pQrKMAe8+e/KU8vLyUnPy8fHx9t/2G7hsYNs1opEt4oXAH+/y0LW7t2nQLq1tX4hAQFBQVq/779Gho3zLKJbKtWN6ptu993bP/qa+vjt2+95VadOPHrVXgzACjdniO/Km7KQst98pJPJEl3395O/3j0nsoKCwCsfPXJZ/INCNCQcY8oIChQR/Yf0LPDxyg76/fkb3BYPau9DwOCg/Xav96z3N8zIl73jIjXrq3bNXHIg5Kk2nXqKP6JcQoMDdG57Bx989nnWjbnZRVRXg2UqvquF3Eeg7mMO7W6u7vr5MmTCg4OLtPEPj4+Sk1NVZMmTcrUPyK8fDtuA8C16tgnsys7BABwqp79p1Z2CADgVJ8e2l3ZIVwV14dHOG2un44fc9pcVUmZV56YzWYtXLhQ3t7eZepfWFhY7qAAAAAAAACuFWVOnjRq1EhvvfVWmScODQ1VzZoVf4wsAAAAAAAoXXXe6NVZypw8OXbsWIntf1T9sOErAAAAAADXHlInjrPrqOL/tWjRIrVq1UoeHh7y8PBQq1attHDhwisPBAAAAAAAqELKlTyZMmWKxo0bp7vuukvvvvuu3n33Xd11110aP368pkyZ4uwYAQAAAABAOZmdeJXHa6+9poiICHl4eKhjx47atm2bzf7vvvuubrjhBnl4eOimm27SJ598Yv0+ZrOmTJmievXqydPTUzExMTp06FA5oyubciVP3njjDb311ltKTk5W37591bdvXyUnJ2vBggV6/fXXnR0jAAAAAAAop8pMnqxevVoJCQlKSkrSzp071aZNG8XGxiozM7PE/t99953uu+8+DR8+XP/+97/Vv39/9e/fX3v27LH0eeGFF/Tyyy9r/vz52rp1q2rXrq3Y2FhdvHixHBGWTbmSJ4WFhWrfvv1l7VFRUbrE+eoAAAAAAEDSnDlzNHLkSMXHx6tly5aaP3++vLy8tHjx4hL7z5s3Tz169NCTTz6pFi1aaMaMGWrXrp1effVVSb+vOpk7d64mT56sfv36qXXr1lq+fLl+++03rVmzpsLeo1zJkwceeEBvvPHGZe0LFizQ4MGDHQ4KAAAAAAA4R2WtPCkoKNCOHTsUExNjaXNzc1NMTIxSUlJKHJOSkmLVX5JiY2Mt/Y8ePar09HSrPr6+vurYsWOpczpDmU/b+bNFixZpw4YN6tSpkyRp69atSktLU1xcnBISEiz95syZ43iUAAAAAACg0plMJplMJqs2o9Eoo9F4Wd/Tp0+rqKhIISEhVu0hISE6cOBAifOnp6eX2D89Pd3y/I+20vpUhHIlT/bs2aN27dpJko4cOSJJCgwMVGBgoFUdEscXAwAAAADgOpKTkzVt2jSrtqSkJE2dOrVyArpKypU82bx5s7PjAAAAAAAAFcJ5CxsSExOtqk0klbjqRPp9kYW7u7syMjKs2jMyMhQaGlrimNDQUJv9//jfjIwM1atXz6pPZGSkXe9ij3LteQIAAAAAAKoKg9Muo9EoHx8fq6u05EmtWrUUFRWlTZs2WdqKi4u1adMmRUdHlzgmOjraqr8kbdy40dK/cePGCg0NteqTm5urrVu3ljqnM5R7zxMAAAAAAABbEhISNHToULVv314dOnTQ3LlzdeHCBcXHx0uS4uLiVL9+fSUnJ0uSxo0bpy5duujFF19U7969tWrVKv3www9asGCBpN+3B3n88cf197//Xdddd50aN26sZ599VmFhYerfv3+FvQfJEwAAAAAAUCEGDhyoU6dOacqUKUpPT1dkZKTWr19v2fA1LS1Nbm7/LYrp3LmzVq5cqcmTJ2vSpEm67rrrtGbNGrVq1crS56mnntKFCxc0atQoZWdn69Zbb9X69evl4eFRYe9hMJvN9p42VCEiwhtXdggA4FTHPpld2SEAgFP17D+1skMAAKf69NDuyg7hqogIb+K0uY4d/9lpc1Ul7HkCAAAAAABgA2U7AAAAAAC4MucdtlNtkTwBAAAAAMClUXTiKL5BAAAAAAAAG1h5AgAAAACACzNQt+MwkicAAAAAALgyA8kTR1G2AwAAAAAAYAMrTwAAAAAAcGGU7TiO5AkAAAAAAC6NohNH8Q0CAAAAAADYwMoTAAAAAABcmIENYx1G8gQAAAAAAFdmoOjEUXyDAAAAAAAANrDyBAAAAAAAF2Zg3YTDSJ4AAAAAAODC2PPEcaSfAAAAAAAAbGDlCQAAAAAArowNYx1G8gQAAAAAABdmIHniML5BAAAAAAAAG1h5AgAAAACAC+O0HceRPAEAAAAAwIVRtuM4vkEAAAAAAAAbSJ4AAAAAAADYQNkOAAAAAAAuzGBwr+wQqjxWngAAAAAAANjAyhMAAAAAAFwYG8Y6juQJAAAAAAAujOSJ4/gGAQAAAAAAbGDlCQAAAAAALowNYx1H8gQAAAAAABdG2Y7j+AYBAAAAAABsYOUJAAAAAAAujLIdx5E8AQAAAADAhZE8cRxlOwAAAAAAADaw8gQAAAAAABfmxoaxDiN5AgAAAACAC6Nsx3GknwAAAAAAAGxg5QkAAAAAAC6MlSeOI3kCAAAAAIALI3niOMp2AAAAAAAAbGDlCQAAAAAALszgxsoTR7HyBAAAAAAAF+ZmcHfaVVHOnDmjwYMHy8fHR35+fho+fLjOnz9vs/+jjz6q5s2by9PTU40aNdJjjz2mnJwcq34Gg+Gya9WqVXbHx8oTAAAAAABQqQYPHqyTJ09q48aNKiwsVHx8vEaNGqWVK1eW2P+3337Tb7/9ptmzZ6tly5Y6fvy4xowZo99++03vvfeeVd8lS5aoR48elns/Pz+74yN5AgAAAACAC7vWN4zdv3+/1q9fr+3bt6t9+/aSpFdeeUW9evXS7NmzFRYWdtmYVq1a6f3337fcN23aVM8995yGDBmiS5cuqUaN/6Y7/Pz8FBoa6lCM10zy5LFVayo7BABwqp79h1R2CADgVJ+umVrZIQAAKpnJZJLJZLJqMxqNMhqN5Z4zJSVFfn5+lsSJJMXExMjNzU1bt27V3XffXaZ5cnJy5OPjY5U4kaRHHnlEI0aMUJMmTTRmzBjFx8fLYDDYFSN7ngAAAAAAgDJJTk6Wr6+v1ZWcnOzQnOnp6QoODrZqq1GjhgICApSenl6mOU6fPq0ZM2Zo1KhRVu3Tp0/XO++8o40bN2rAgAF6+OGH9corr9gd4zWz8gQAAAAAADifM8t2EhMTlZCQYNVW2qqTp59+WjNnzrQ53/79+x2OKTc3V71791bLli01depUq2fPPvus5Z/btm2rCxcuaNasWXrsscfs+gySJwAAAAAAuDCDwXl/9LenRGfChAkaNmyYzT5NmjRRaGioMjMzrdovXbqkM2fOXHGvknPnzqlHjx6qU6eOPvzwQ9WsWdNm/44dO2rGjBkymUx2lRqRPAEAAAAAAE4XFBSkoKCgK/aLjo5Wdna2duzYoaioKEnSF198oeLiYnXs2LHUcbm5uYqNjZXRaNS//vUveXh4XPGzUlNT5e/vb/ceLSRPAAAAAABwYW7X+Gk7LVq0UI8ePTRy5EjNnz9fhYWFGjt2rAYNGmQ5aefXX39Vt27dtHz5cnXo0EG5ubnq3r278vLy9M9//lO5ubnKzc2V9HvSxt3dXR9//LEyMjLUqVMneXh4aOPGjXr++ef1xBNP2B0jyRMAAAAAAFyYwe3aTp5I0ooVKzR27Fh169ZNbm5uGjBggF5++WXL88LCQh08eFB5eXmSpJ07d2rr1q2SpGbNmlnNdfToUUVERKhmzZp67bXXNH78eJnNZjVr1kxz5szRyJEj7Y6P5AkAAAAAAKhUAQEBWrlyZanPIyIiZDabLfddu3a1ui9Jjx491KNHD6fER/IEAAAAAAAX5swNY6srvkEAAAAAAFyYM48qrq7cKjsAAAAAAACAaxkrTwAAAAAAcGGU7TiObxAAAAAAABd2rR9VXBVQtgMAAAAAAGADK08AAAAAAHBhBjf+6O8ovkEAAAAAAFwYe544jrIdAAAAAAAAG0g/AQAAAADgwgxsGOswkicAAAAAALgwynYcR9kOAAAAAACADaSfAAAAAABwYZy24zhWngAAAAAAANhA+gkAAAAAABfGnieO4xsEAAAAAMCVkTxxGGU7AAAAAAAANpA8AQAAAAAAsIG1OwAAAAAAuDBO23EcK08AAAAAAABsIP0EAAAAAIAL47Qdx/ENAgAAAADgyijbcRhlOwAAAAAAADaQfgIAAAAAwJUZ3Cs7giqP5AkAAAAAAC6M03YcR9kOAAAAAACADaSfAAAAAABwZZy24zC+QQAAAAAAXJiZsh2HUbYDAAAAAABgA+knAAAAAABcmRun7TiK5AkAAAAAAK6M5InDKNsBAAAAAACwgZUnAAAAAAC4MDMrTxxG8gQAAAAAABdG8sRxlO0AAAAAAADYwMoTAAAAAABcGStPHEbyBAAAAAAAF2Z2o+jEUXyDAAAAAAAANrDyBAAAAAAAF8aGsY5j5QkAAAAAAIANJE8AAAAAAABsIHkCAAAAAIALK3Z3c9pVUc6cOaPBgwfLx8dHfn5+Gj58uM6fP29zTNeuXWUwGKyuMWPGWPVJS0tT79695eXlpeDgYD355JO6dOmS3fGx5wkAAAAAAC6sKpy2M3jwYJ08eVIbN25UYWGh4uPjNWrUKK1cudLmuJEjR2r69OmWey8vL8s/FxUVqXfv3goNDdV3332nkydPKi4uTjVr1tTzzz9vV3wkTwAAAAAAQKXZv3+/1q9fr+3bt6t9+/aSpFdeeUW9evXS7NmzFRYWVupYLy8vhYaGlvhsw4YN2rdvnz7//HOFhIQoMjJSM2bM0MSJEzV16lTVqlWrzDFe++knAAAAAABQbmY3N6ddJpNJubm5VpfJZHIovpSUFPn5+VkSJ5IUExMjNzc3bd261ebYFStWKDAwUK1atVJiYqLy8vKs5r3pppsUEhJiaYuNjVVubq727t1rV4xlXnny8ssv2zWxJMXHx6tOnTp2jwMAAAAAAM5R7MSyneTkZE2bNs2qLSkpSVOnTi33nOnp6QoODrZqq1GjhgICApSenl7quPvvv1/h4eEKCwvTrl27NHHiRB08eFAffPCBZd7/TZxIstzbmrckZU6ePP7442rQoIHc3ct2PvQvv/yiPn36kDwBAAAAAMBFJCYmKiEhwarNaDSW2Pfpp5/WzJkzbc63f//+cscyatQoyz/fdNNNqlevnrp166YjR46oadOm5Z63JHbtefLDDz9clg0qDUkTAAAAAAAqn9mJp+QYjcZSkyV/NmHCBA0bNsxmnyZNmig0NFSZmZlW7ZcuXdKZM2dK3c+kJB07dpQkHT58WE2bNlVoaKi2bdtm1ScjI0OS7JpXsiN5kpSUJG9v7zJPPGnSJAUEBNgVDAAAAAAAcC6zm6FSPjcoKEhBQUFX7BcdHa3s7Gzt2LFDUVFRkqQvvvhCxcXFloRIWaSmpkqS6tWrZ5n3ueeeU2ZmpmUhyMaNG+Xj46OWLVva9S5lTj8lJSVZHflzJYmJifLz87MrGAAAAAAAUL20aNFCPXr00MiRI7Vt2zZ9++23Gjt2rAYNGmQ5aefXX3/VDTfcYFlJcuTIEc2YMUM7duzQsWPH9K9//UtxcXH6y1/+otatW0uSunfvrpYtW+qBBx7Qjz/+qM8++0yTJ0/WI488UubVM3/gqGJUaWazWT98+I4OfLlJprwLCr3uBt0WN0K+ofVKHVNcXKwdH76jQylfKy8nW7X9AnT9rV3Uru8AGQz/zcie/e2Etr6zQicP7lNxUbH86zfQnWMnqE7dwKvxagCqqT6DB+meEcPkHxSonw8c1BvTk/XTrj0l9m3UrKkeePwRXXdjS4U0qK83n5upNUv/adXHs7aX4h4fq+g7u8mvboCO7DugN//+D/20274d5gGgom3fe1SLPvpae478qlNnz+m1iUMU09G+vxkGULJi98pZeWKPFStWaOzYserWrZvc3Nw0YMAAq4NrCgsLdfDgQctpOrVq1dLnn3+uuXPn6sKFC2rYsKEGDBigyZMnW8a4u7tr7dq1euihhxQdHa3atWtr6NChmj59ut3xOTV5sn//fvXu3Vs///yzM6cFSvXjJx9pz8ZPdfvIR1QnKFjbP1itdS8+p3ufm6MapZzZnbpujfZt3qiuIx5RQP0GOnXsZ21Z9LpqeXnppjt7SZJyMtP10XNTdMNf7lD7u+9VTU9Pnf31hGrUrHk1Xw9ANfOXXrEaNelJvTJlhg7+uEv9hz6gvy9+UyO736WcM2cu6+/h6aH0X07om083aNSkp0qcc9xz0xRxfTPNfnKSsjIydUe/Pnp+2Vsa3bO/sjIySxwDAJUhz1Sg5hGhGnBHlMa+sKKywwFcSmWV7dgjICBAK1euLPV5RESEzGaz5b5hw4b68ssvrzhveHi4PvnkE4fjc96uMZIKCgp0/PhxZ04JlMpsNmv3hk/Uru9fFdHuZtVtGK7bR45V3tmzOrZze6njMg7/pPC27RUe2U51goLV5OZOanBja2X+fNjSZ/t7q9SodVt1GjhEgeGN5Rscqoi27eXp43s1Xg1ANXX3g3H6dPX72vj+GqUd/lmvTJkuU36+ut9zd4n9f9q9V4tmztGX69arsKDgsue1jEbdGhujRS/M0Z7tO3Qy7ReteOUN/Xb8F/W+f2BFvw4A2KVLu+Yaf3933dnpxsoOBQAuY9fKkz8fR/Rnp06dcigYwB7nTmUqLydb9Vu2trQZvbwU3LSZMo78pGadbilxXEiz67V/yyZlp/8mv9AwZaUdU/qhg4oeFCdJMhcXK23XTrXp2VfrZj+n08ePyicoWJG9+6txVIer8m4Aqp8aNWvouhtb6p35iyxtZrNZqd99rxZt25RrTvca7nKvUUOFJuvESsHFi7oxqq1D8QIAgKqjKqw8udbZlTyZN2+eIiMj5ePjU+Lz8+fPOyUooCzycrIlSZ6+1qtBPH18Lc9K0rZ3fxXm52t14ni5ubmpuLhYHQYM0nWdb5Mk5efmqvDiRaWu+0g3Dxiojn8brF92p2rDqy/qrolJCruB2lsAzufj7y/3GjV09nSWVfvZrCw1aNq4XHPmX8jTvp2puu+R0Uo78rOyT2epS59euqFtG508nuaMsAEAQBVgdq/sCKo+u5InzZo10/jx4zVkyJASn6emplqOFbLFZDLJZDJZtV0qKCh1jwpAkg5997W+WrbAct9zfGK55jmyLUWHvv9G3UY/Jv/6DZWVdkzfrVwqLz9/Nb+1q8zmYklSRLv2ah3bR5IUGB6hjMMHtW/zBpInAKqU2U8manzyDK349gsVXbqkw3v368u1n6pZK37LAAAAysqu5En79u21Y8eOUpMnBoPBagOX0iQnJ2vatGlWbd0fHK3YEQ/ZEw6qmfC27XVP0+ss90WXCiVJ+Tk5qu3nb2nPz81R3UYRpc7z/Tv/VGSvfpaynroNG+l81imlrl2j5rd2lUcdH7m5u8s/rIHVOL+w+kr/6aAT3wgA/iv37FkVXbok/8C6Vu3+devq7KmsUkZd2cm0E3pqcLyMnp7y8q6ts6dO6+m5s5T+ywlHQwYAAFUEZTuOs2vD2BdffFGPP/54qc/btGmj4uLiK86TmJionJwcq6tb3HB7QkE1VMvTU74hoZbLP6yBvHz99Ou+3ZY+Bfl5yjxyWCFNry91nksmkwxu1v/XN7i5WRJ/7jVqKKhxU2Wf/M2qT076SdUJ5JhiABXjUuElHdq7T5HRHS1tBoNBkZ07af+/f3R4flN+vs6eOi1vHx9F3dZZ33++2eE5AQBAFeHmxKuasmvlSWhoqFM+1Gg0ymg0WgdCyQ7sZDAYdFP3Xtr58QfyDa2nOoHB+uGDVfLy91dEu5st/T6eOV2NozqoVUwPSVJ4ZJT+/fEH8g4IVED9Bjqddky7Plur5rfdbhnTpmdfff76S6rXvIXCWrTSL7tTdTx1h+56eurVfk0A1ciHi5drwgvP6dCevTq4a7f6D3tARk9PbXx/jSRpwgvPKSsjU0tfnCfp901mGzVr+p9/rqm6IcFq0qK58i/k6WTaL5Kkdrd2lsFg0ImjxxQW3kjDJyboxM9HteE/cwLAteJCvklp6f9daXci84z2H/1Nvt5eCgvyq7zAAEB2Jk+Aa02bXv1UaDLpqyVvqiAvT6HX36BeEyZZJeNyMzN08Vyu5f6WIQ9q+wer9c3bC5Wfm6PafgFq0fVORfW7x9KncVQH3TZ0pP69bo2+XbFEfqFh6j52gupdf8NVfT8A1ctXn3wm34AADRn3iAKCAnVk/wE9O3yMsrN+/8NEcFg9q/LYgOBgvfav9yz394yI1z0j4rVr63ZNHPKgJKl2nTqKf2KcAkNDdC47R9989rmWzXlZRZcuXd2XA4Ar2HPkV8VNWWi5T17yiSTp7tvb6R+P3lPaMABlwYaxDjOYy7JJiaSAgAD99NNPCixj2UKjRo309ddfKzw8vEz956Q4viQZAK4lG+NK3h8KAKqqT9dMrewQAMC5bhxQ2RFcFZ2mO+/P299PaeO0uaqSMq88yc7O1qeffirfPx0LW5qsrCwVFRWVOzAAAAAAAIBrgV1lO0OHDq2oOAAAAAAAQEWoxhu9OkuZkydlOUUHAAAAAADA1ZR7w9hNmzZp06ZNyszMtEqsGAwGLVq0yCnBAQAAAAAAVLZyJU+mTZum6dOnq3379qpXr54MBoOz4wIAAAAAAE5goGzHYeVKnsyfP19Lly7VAw884Ox4AAAAAACAExncynTILmwoV/6poKBAnTt3dnYsAAAAAAAA15xyJU9GjBihlStXOjsWAAAAAADgZAY3513VVbnKdi5evKgFCxbo888/V+vWrVWzZk2r53PmzHFKcAAAAAAAwDFu7pUdQdVXruTJrl27FBkZKUnas2eP1TM2jwUAAAAAAK6kXMmTzZs3OzsOAAAAAABQAdyqcbmNs5QreQIAAAAAAKoGTttxHPknAAAAAAAAG1h5AgAAAACAC6Nsx3EkTwAAAAAAcGEkTxzHVwgAAAAAAGADK08AAAAAAHBhrDxxHMkTAAAAAABcGMkTx/EVAgAAAAAA2MDKEwAAAAAAXBgrTxxH8gQAAAAAABfm7mau7BCqPPJPAAAAAAAANrDyBAAAAAAAF0bZjuNIngAAAAAA4MJInjiOrxAAAAAAAMAGkicAAAAAAAA2ULYDAAAAAIALc2fZhMP4CgEAAAAAAGxg5QkAAAAAAC7MzVDZEVR9JE8AAAAAAHBhlO04jq8QAAAAAADABpInAAAAAAC4MDc3510V5cyZMxo8eLB8fHzk5+en4cOH6/z586X2P3bsmAwGQ4nXu+++a+lX0vNVq1bZHR9lOwAAAAAAuLCqULYzePBgnTx5Uhs3blRhYaHi4+M1atQorVy5ssT+DRs21MmTJ63aFixYoFmzZqlnz55W7UuWLFGPHj0s935+fnbHR/IEAAAAAABUmv3792v9+vXavn272rdvL0l65ZVX1KtXL82ePVthYWGXjXF3d1doaKhV24cffqh7771X3t7eVu1+fn6X9bVXFcg/AQAAAACA8nJ3c95VEVJSUuTn52dJnEhSTEyM3NzctHXr1jLNsWPHDqWmpmr48OGXPXvkkUcUGBioDh06aPHixTKbzXbHyMoTAAAAAABcmDOTHiaTSSaTyarNaDTKaDSWe8709HQFBwdbtdWoUUMBAQFKT08v0xyLFi1SixYt1LlzZ6v26dOn64477pCXl5c2bNighx9+WOfPn9djjz1mV4ysPAEAAAAAAGWSnJwsX19fqys5ObnEvk8//XSpm7r+cR04cMDhmPLz87Vy5coSV508++yzuuWWW9S2bVtNnDhRTz31lGbNmmX3Z7DyBAAAAAAAF+bMU3ISExOVkJBg1VbaqpMJEyZo2LBhNudr0qSJQkNDlZmZadV+6dIlnTlzpkx7lbz33nvKy8tTXFzcFft27NhRM2bMkMlksmu1DMkTAAAAAABcmLvBeXPZU6ITFBSkoKCgK/aLjo5Wdna2duzYoaioKEnSF198oeLiYnXs2PGK4xctWqS+ffuW6bNSU1Pl7+9vd5kRyRMAAAAAAFBpWrRooR49emjkyJGaP3++CgsLNXbsWA0aNMhy0s6vv/6qbt26afny5erQoYNl7OHDh/XVV1/pk08+uWzejz/+WBkZGerUqZM8PDy0ceNGPf/883riiSfsjpHkCQAAAAAALqyiTslxphUrVmjs2LHq1q2b3NzcNGDAAL388suW54WFhTp48KDy8vKsxi1evFgNGjRQ9+7dL5uzZs2aeu211zR+/HiZzWY1a9ZMc+bM0ciRI+2Oz2Auzxk9FWBOyo+VHQIAONXGuCGVHQIAONWna6ZWdggA4Fw3DqjsCK6Kcet2Om2ueb3bOW2uqqQK5J8AAAAAAAAqD2U7AAAAAAC4sBpuTtwxtpoieQIAAAAAgAurCnueXOv4CgEAAAAAAGxg5QkAAAAAAC7Mnaodh5E8AQAAAADAhVG24zi+QgAAAAAAABtIngAAAAAAANhA2Q4AAAAAAC6Msh3H8RUCAAAAAADYwMoTAAAAAABcmLsbx+04iuQJAAAAAAAujLIdx/EVAgAAAAAA2MDKEwAAAAAAXJg7VTsOI3kCAAAAAIALY88Tx1G2AwAAAAAAYAMrTwAAAAAAcGFsGOs4g9lsNld2EMDVYjKZlJycrMTERBmNxsoOBwAcxu8aAFfD7xqAaxHJE1Qrubm58vX1VU5Ojnx8fCo7HABwGL9rAFwNv2sArkUs3gEAAAAAALCB5AkAAAAAAIANJE8AAAAAAABsIHmCasVoNCopKYnNxwC4DH7XALgaftcAXIvYMBYAAAAAAMAGVp4AAAAAAADYQPIEAAAAAADABpInAAAAAAAANpA8gUuLiIiQwWCQwWBQdnZ2mcctXbrUMu7xxx+vsPgAwF7l/V2bOnWqZdzcuXMrLD4AsNcfv01+fn52jeN3DcDVRPIELm/69Ok6efKkfH19JUkXL17UsGHDdNNNN6lGjRrq37//ZWMGDhyokydPKjo6+ipHCwBX9ufftS1btqhfv36qV6+eateurcjISK1YscJqzBNPPKGTJ0+qQYMGlREyANi0ZMkS/fTTT5b7kydP6v7779f1118vNze3Ev8yi981AFcTyRO4vDp16ig0NFQGg0GSVFRUJE9PTz322GOKiYkpcYynp6dCQ0NVq1atqxkqAJTJn3/XvvvuO7Vu3Vrvv/++du3apfj4eMXFxWnt2rWWMd7e3goNDZW7u3tlhQ0ApfLz81NwcLDl3mQyKSgoSJMnT1abNm1KHMPvGoCrqUZlBwA4omvXrmrVqpUk6e2331bNmjX10EMPafr06ZY/VPxZ7dq19cYbb0iSvv32W7uWvQNARSvP79qkSZOs7seNG6cNGzbogw8+UJ8+fSo8ZgCwpWvXrmrdurU8PDy0cOFC1apVS2PGjNHUqVNLHRMREaF58+ZJkhYvXnyVIgWA0rHyBFXesmXLVKNGDW3btk3z5s3TnDlztHDhwsoOCwDKzRm/azk5OQoICKigCAHAPsuWLVPt2rW1detWvfDCC5o+fbo2btxY2WEBQJmx8gRVXsOGDfXSSy/JYDCoefPm2r17t1566SWNHDmyskMDgHJx9HftnXfe0fbt2/Xmm29WcKQAUDatW7dWUlKSJOm6667Tq6++qk2bNunOO++s5MgAoGxYeYIqr1OnTlZL2aOjo3Xo0CEVFRVVYlQAUH6O/K5t3rxZ8fHxeuutt3TjjTdWZJgAUGatW7e2uq9Xr54yMzMrKRoAsB/JEwAAXMSXX36pu+66Sy+99JLi4uIqOxwAsKhZs6bVvcFgUHFxcSVFAwD2I3mCKm/r1q1W999//72uu+46dl4HUGWV53dty5Yt6t27t2bOnKlRo0ZVdIgAAADVCnueoMpLS0tTQkKCRo8erZ07d+qVV17Riy++aHPMvn37VFBQoDNnzujcuXNKTU2VJEVGRlZ8wABwBfb+rm3evFl9+vTRuHHjNGDAAKWnp0uSatWqxaaxAKqsP/777Pz58zp16pRSU1NVq1YttWzZsnIDA1AtkTxBlRcXF6f8/Hx16NBB7u7uGjdu3BX/1rVXr146fvy45b5t27aSJLPZXKGxAkBZ2Pu7tmzZMuXl5Sk5OVnJycmW9i5dumjLli1XIWIAcL4//vtMknbs2KGVK1cqPDxcx44dq7ygAFRbJE9Q5dWsWVNz587VG2+8UeYx/EsXwLXM3t+1pUuXaunSpRUbFACUU0lJ3DVr1lxxHH+pBeBawp4ncHkTJ06Ut7e3cnJyyjxmxYoV8vb21tdff12BkQFA+ZTnd+3555+Xt7e30tLSKjAyACif++67Tw0aNLBrDL9rAK4mVp7ApX355ZcqLCyUJNWpU6fM4/r27auOHTtKkvz8/CoiNAAol/L+ro0ZM0b33nuvJCkoKKhCYgOA8jh06JAk2b3ZP79rAK4mg5n1cAAAAAAAAKWibAcAAAAAAMAGkicAAAAAAAA2kDwBAAAAAACwgeQJAAAAAACADSRPAAAAAAAAbCB5AgAAAAAAYAPJEwAAAAAAABtIngAAAAAAANhA8gQAAAAAAMCG/weAyEmeOcThwgAAAABJRU5ErkJggg==\n" + }, + "metadata": {} + } + ], + "source": [ + "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(15, 5))\n", + "sns.heatmap(df[['p[1]','p[2]','n[1]']].corr(),annot=True, center=0,ax=axes)\n", + "\n", + "axes.set_title('Correlations')\n", + "plt.show()" ] - }, - "execution_count": 173, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "from sklearn.ensemble import GradientBoostingRegressor\n", - "xgb_regressor = make_pipeline(GradientBoostingRegressor(n_estimators=25)) \n", - "xgb_regressor.fit(X_train, y_train) \n", - "xgb_regressor_validation = cross_validate(xgb_regressor, X_train, y_train, cv=5, return_train_score=True, return_estimator=True)\n", - "\n", - "xgb_regressor_validation['train_score'], xgb_regressor_validation['test_score']\n" - ] - }, - { - "cell_type": "markdown", - "id": "mTquaGiJF2pO", - "metadata": { - "id": "mTquaGiJF2pO" - }, - "source": [ - "## Price optimization model with competing products\n", - "\n", - "Our problem is to:\n", - "1.\tDetermine the number of each category of product to make available given the overall restriction of what we can offer. \n", - "2.\tWe are also instructed to make sure there are a minimum number of each category made available as well as a minimum and maximum price for each category.\n", - "3.\tLastly, the product categories should be decreasing in price, meaning Category 1 should be the most expensive, and so on. Specifically, we must make sure there is at least a $50 gap between categories, but no more than $100. \n", - " \n", - "With the predictive part in place, it's time to build the optimization model. The model is formulated (i.e. the mathematical representation) for an unspecified number of categories, but the code will reflect that we have two categories of products in this problem. We start by setting some parameter values (not to be confused with ML hyperparameters) and initialize the optimization model. \n", - "\n" - ] - }, - { - "cell_type": "markdown", - "id": "95b966b2", - "metadata": {}, - "source": [ - "### Initialize model and set input parameters\n", - "- $C$: Number of product categories\n", - "- $N$: Total amount of space available\n", - "- $\\lambda$: Price control parameter\n", - "\n", - "Here is the first mention of a price control parameter. It is fairly common in optimization modeling to add penalty terms to try and prevent undesirable outcomes. This is akin to using penalty terms in machine learning and applied statistics to prevent overfitting, with [Lasso](https://en.wikipedia.org/wiki/Lasso_(statistics)) and [Ridge](https://en.wikipedia.org/wiki/Ridge_regression) regression as a couple of common examples. " - ] - }, - { - "cell_type": "code", - "execution_count": 174, - "id": "67db9804", - "metadata": {}, - "outputs": [ + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "Set parameter Username\n", - "Set parameter WLSAccessID\n", - "Set parameter WLSSecret\n", - "Set parameter LicenseID to value 874302\n", - "WLS license 874302 - registered to Gurobi Optimization LLC\n" - ] - } - ], - "source": [ - "#### Initialize the model\n", - "m = gp.Model(\"price optimization\")\n", - "\n", - "products = [1,2] #### Category 1 and Category 2\n", - "N = 200 #### limit on available space\n", - "l = 0 #### price control, we'll start this at 0" - ] - }, - { - "cell_type": "markdown", - "id": "ac5a875b", - "metadata": {}, - "source": [ - "### Decision variables\n", - "- $p_c$: price per item in category $c = 1,2,\\dots, C$\n", - "- $n_c$: number of items allocated to category, predicted using features $p_c$, $c = 1,2,\\dots, C$" - ] - }, - { - "cell_type": "code", - "execution_count": 175, - "id": "bbe94a3c", - "metadata": {}, - "outputs": [], - "source": [ - "p = m.addVars(products, name=\"p\") #### price decision variables\n", - "n = m.addVars(products, name=\"n\") #### decision variable for number of items in each category" - ] - }, - { - "cell_type": "markdown", - "id": "721614e0", - "metadata": {}, - "source": [ - "### Constraints" - ] - }, - { - "cell_type": "markdown", - "id": "8403f941", - "metadata": {}, - "source": [ - "We need to have a minimum number of each category available. \n", - "\\begin{align*}\n", - "n_c \\ge l_{n_c}\n", - "\\end{align*}\n", - "\n", - "We also set lower and upper bounds on the prices.\n", - "\\begin{align*}\n", - "l_{n_c} \\le p_c \\le u_{p_c}\n", - "\\end{align*}" - ] - }, - { - "cell_type": "code", - "execution_count": 176, - "id": "edec4432", - "metadata": {}, - "outputs": [], - "source": [ - "min_items = {1:50,2:50}\n", - "price_bounds = {1:[300,400], 2:[100,300]}\n", - "m.addConstrs(n[c] >= min_items[c] for c in products) #### we could hardcode 50 instead of min_items, but this is more flexible\n", - "m.addConstr(p[1] == [300,400]) #### this is a shorthand way to code 300 <= p[1] <= 400 \n", - "m.addConstr(p[2] == [100,300]);" - ] - }, - { - "cell_type": "markdown", - "id": "3bef9e0d", - "metadata": {}, - "source": [ - "Another note: each of the above constraints can be addressed when defining the decision variables. Here is an example for the decision variable $n$." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e2065d83", - "metadata": {}, - "outputs": [], - "source": [ - "#price_lb = {1:300, 2:100}\n", - "#price_ub = {1:400, 2:300}\n", - "#p = m.addVars(products, lb = price_lb, ub = price_ub, name=\"p\") #### each price is now bounded\n", - "#n = m.addVars(products, lb = min_items, name=\"n\") " - ] - }, - { - "cell_type": "markdown", - "id": "0578e271", - "metadata": {}, - "source": [ - "In general, the number of items allocated must equal the total available space. \n", - "\\begin{equation*}\n", - "n_1 + n_2 + \\dots + n_C = \\sum_{c}n_c = N \\\\\n", - "\\end{equation*}\n", - "Note that this, along with the constraint on the minimum number available means we don't have to specify an upper bound for each $n_c$." - ] - }, - { - "cell_type": "code", - "execution_count": 177, - "id": "140e739f", - "metadata": {}, - "outputs": [], - "source": [ - "m.addConstr(n.sum() == N); #### remember we set N = 200 earlier" - ] - }, - { - "cell_type": "markdown", - "id": "b17a1e78", - "metadata": {}, - "source": [ - "The last set of constraints are for price ordering. This requires the subsequent category to be cost between $50 and $100 less than the previous. \n", - "\\begin{equation*}\n", - "50 \\le p_c - p_{c+1} \\le 100\n", - "\\end{equation*}" - ] - }, - { - "cell_type": "code", - "execution_count": 178, - "id": "1428be22", - "metadata": {}, - "outputs": [], - "source": [ - "m.addConstr(p[1]-p[2] == [50,100]);" - ] - }, - { - "cell_type": "markdown", - "id": "51ee2204", - "metadata": {}, - "source": [ - "### Objective function" - ] - }, - { - "cell_type": "markdown", - "id": "6eb69e54", - "metadata": {}, - "source": [ - "We want to maximize total revenue with the portion of total revenue coming from category $c$ being $p_cn_c$. This makes the total revenue $\\sum_{c} p_c n_c$. That is the first part of the objective. Earlier a price control parameter was introduced which is the second part of the objective. The lambda parameter captures the trade-off between the revenue and price-control pieces. This term penalizes the model from setting too high of prices since doing so could lose sales. Our model assumed we'll sell all of the items so having this penalty term can make this assumption more realistic. For reference, [here is a good source](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4565407).\n", - "\n", - "This term will be defined as $λ (\\sum_{c} p_c^2)$ for this problem. So, the complete objective is:\n", - "\\begin{equation*}\n", - "\\textrm{maximize} \\sum_{c} p_c n_c - λ (\\sum_{c} p_c^2)\n", - "\\end{equation*}" - ] - }, - { - "cell_type": "code", - "execution_count": 179, - "id": "a2cfcabc", - "metadata": {}, - "outputs": [], - "source": [ - "revenue = gp.quicksum(p[c]*n[c] for c in products) #### you could also use the more simple p.prod(n)\n", - "penalty = l*(p[1]**2+p[1]**2) #### we used l as the lambda parameter earlier\n", - "m.setObjective(revenue - penalty, sense = GRB.MAXIMIZE)" - ] - }, - { - "cell_type": "markdown", - "id": "c707af60", - "metadata": {}, - "source": [ - "### Integrate the ML model\n" - ] - }, - { - "cell_type": "markdown", - "id": "38d41039", - "metadata": {}, - "source": [ - "Right now, if we were to run the optimization, the solution would be to set the price for Category 1 to $400, Category 2 to $300, and sell 150 and 50 of each item, respectively. That's because we have yet to add in the relationship between price and demand that was derived from the ML model. To integrate the machine learning model into the optimization model, we'll use the Gurobi Machine Learning package. The magic happens using `add_predictor_constr` function. " - ] - }, - { - "cell_type": "code", - "execution_count": 180, - "id": "b113eda2", - "metadata": {}, - "outputs": [ + "cell_type": "markdown", + "id": "a4aa96d4", + "metadata": { + "id": "a4aa96d4" + }, + "source": [ + "In this problem, we've assumed that the amount of space we have available for the products is 200 units. In retail, this could be the amount of warehouse space, or for ticketing this could represent the number of seats available." + ] + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "Requirement already satisfied: gurobi-machinelearning in /Users/yurchisin/opt/anaconda3/envs/gurobi_ml/lib/python3.11/site-packages (1.3.3)\n", - "Requirement already satisfied: numpy>=1.22.0 in /Users/yurchisin/opt/anaconda3/envs/gurobi_ml/lib/python3.11/site-packages (from gurobi-machinelearning) (1.26.2)\n", - "Requirement already satisfied: gurobipy>=10.0.0 in /Users/yurchisin/opt/anaconda3/envs/gurobi_ml/lib/python3.11/site-packages (from gurobi-machinelearning) (11.0.0)\n", - "Requirement already satisfied: scipy>=1.9.3 in /Users/yurchisin/opt/anaconda3/envs/gurobi_ml/lib/python3.11/site-packages (from gurobi-machinelearning) (1.11.4)\n", - "Note: you may need to restart the kernel to use updated packages.\n" - ] - } - ], - "source": [ - "#### install the package and load the required function\n", - "%pip install gurobi-machinelearning\n", - "from gurobi_ml import add_predictor_constr" - ] - }, - { - "cell_type": "markdown", - "id": "119a8f7e", - "metadata": {}, - "source": [ - "This additional package is useful when we have **decision variables** that are also **features** of a machine learning model. First, we need a data frame that contains these decision variables. It is important to make sure the indices of the data frame have the **same name** as the training data for the machine learning model. " - ] - }, - { - "cell_type": "code", - "execution_count": 182, - "id": "0b1ac8cc", - "metadata": {}, - "outputs": [], - "source": [ - "m_feats = pd.DataFrame({\"p[1]\":[p[1]],\"p[2]\":[p[2]]})" - ] - }, - { - "cell_type": "markdown", - "id": "3263a396", - "metadata": {}, - "source": [ - "Adding the predictive model to the optimization model requires specifying the model we want to use `(m)`, regression object `(xgb_regressor)`, feature data frame `(m_feats)`, and the output decision variable `(n[1])`. Remember `n[2]` is **NOT** the output of the regression. We can then print the number of variables and constraints added to the model using `print_stats`.\n" - ] - }, - { - "cell_type": "code", - "execution_count": 186, - "id": "1fe6c59f", - "metadata": {}, - "outputs": [ + "cell_type": "markdown", + "id": "6e9c6157", + "metadata": { + "id": "6e9c6157" + }, + "source": [ + "### Building regressors to predict sales\n", + "\n", + "The prices for each category item will be used to predict the number of Category 1 items sold. Here we build a regression model to form this relationship which will later be used as part of the optimization model." + ] + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model for pipe0:\n", - "0 variables\n", - "1 constraints\n", - "Input has shape (1, 2)\n", - "Output has shape (1, 1)\n", - "\n", - "Pipeline has 1 steps:\n", - "\n", - "--------------------------------------------------------------------------------\n", - "Step Output Shape Variables Constraints \n", - " Linear Quadratic General\n", - "================================================================================\n", - "lin_reg (1, 1) 0 1 0 0\n", - "\n", - "--------------------------------------------------------------------------------\n" - ] - } - ], - "source": [ - "pred_constr = add_predictor_constr(m, linear_regressor, m_feats, n[1])\n", - "pred_constr.print_stats()" - ] - }, - { - "cell_type": "markdown", - "id": "53ce50f6", - "metadata": {}, - "source": [ - "### Solve the optimization and get the solution" - ] - }, - { - "cell_type": "markdown", - "id": "ec80341d", - "metadata": {}, - "source": [ - "Since this is a quadratic, non-convex problem we set the `NonConvex` parameter to 2. See the [documentation](https://www.gurobi.com/documentation/current/refman/nonconvex.html) for more information. We'll also print out the optimal solution. " - ] - }, - { - "cell_type": "code", - "execution_count": 187, - "id": "d04d5e5b", - "metadata": {}, - "outputs": [ + "cell_type": "code", + "execution_count": 5, + "id": "a70d5f2c", + "metadata": { + "id": "a70d5f2c" + }, + "outputs": [], + "source": [ + "from sklearn.compose import make_column_transformer\n", + "from sklearn.linear_model import LinearRegression\n", + "from sklearn.pipeline import make_pipeline\n", + "from sklearn.metrics import r2_score\n", + "from sklearn.model_selection import train_test_split #importing scikit-learn's function for data splitting\n", + "from sklearn.ensemble import GradientBoostingRegressor #importing scikit-learn's gradient booster regressor function\n", + "from sklearn.model_selection import cross_validate #improting scikit-learn's cross validation function" + ] + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (mac64[rosetta2] - Darwin 22.6.0 22G91)\n", - "\n", - "CPU model: Apple M1\n", - "Thread count: 8 physical cores, 8 logical processors, using up to 8 threads\n", - "\n", - "WLS license 874302 - registered to Gurobi Optimization LLC\n", - "Optimize a model with 33 rows, 232 columns and 240 nonzeros\n", - "Model fingerprint: 0x27842501\n", - "Model has 2 quadratic objective terms\n", - "Model has 648 general constraints\n", - "Variable types: 32 continuous, 200 integer (200 binary)\n", - "Coefficient statistics:\n", - " Matrix range [1e-01, 2e+00]\n", - " Objective range [0e+00, 0e+00]\n", - " QObjective range [2e+00, 2e+00]\n", - " Bounds range [1e+00, 2e+02]\n", - " RHS range [1e+00, 7e+02]\n", - " GenCon rhs range [3e-03, 4e+02]\n", - " GenCon coe range [1e+00, 1e+00]\n", - "\n", - "MIP start from previous solve produced solution with objective 62753.9 (0.04s)\n", - "Loaded MIP start from previous solve with objective 62753.9\n", - "\n", - "Presolve added 111 rows and 0 columns\n", - "Presolve removed 0 rows and 68 columns\n", - "Presolve time: 0.03s\n", - "Presolved: 149 rows, 167 columns, 749 nonzeros\n", - "Presolved model has 2 bilinear constraint(s)\n", - "\n", - "Solving non-convex MIQCP\n", - "\n", - "Variable types: 33 continuous, 134 integer (134 binary)\n", - "\n", - "Root relaxation: objective 6.762746e+04, 191 iterations, 0.00 seconds (0.00 work units)\n", - "\n", - " Nodes | Current Node | Objective Bounds | Work\n", - " Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n", - "\n", - " 0 0 67627.4598 0 31 62753.8592 67627.4598 7.77% - 0s\n", - " 0 0 67627.4598 0 33 62753.8592 67627.4598 7.77% - 0s\n", - " 0 0 67627.4598 0 35 62753.8592 67627.4598 7.77% - 0s\n", - "H 0 0 63247.729112 67627.4598 6.92% - 0s\n", - " 0 0 66606.1147 0 17 63247.7291 66606.1147 5.31% - 0s\n", - " 0 0 66561.7393 0 21 63247.7291 66561.7393 5.24% - 0s\n", - " 0 0 66219.9473 0 45 63247.7291 66219.9473 4.70% - 0s\n", - " 0 0 65902.8627 0 38 63247.7291 65902.8627 4.20% - 0s\n", - " 0 0 65755.2174 0 23 63247.7291 65755.2174 3.96% - 0s\n", - " 0 0 65563.5628 0 11 63247.7291 65563.5628 3.66% - 0s\n", - " 0 0 65456.9640 0 7 63247.7291 65456.9640 3.49% - 0s\n", - "H 0 0 65292.830679 65456.9640 0.25% - 0s\n", - " 0 0 65456.9640 0 7 65292.8307 65456.9640 0.25% - 0s\n", - "\n", - "Cutting planes:\n", - " Learned: 1\n", - " Cover: 4\n", - " Implied bound: 1\n", - " Clique: 35\n", - " MIR: 3\n", - " Flow cover: 4\n", - " Relax-and-lift: 14\n", - "\n", - "Explored 1 nodes (728 simplex iterations) in 0.21 seconds (0.06 work units)\n", - "Thread count was 8 (of 8 available processors)\n", - "\n", - "Solution count 3: 65292.8 63247.7 62753.9 \n", - "\n", - "Optimal solution found (tolerance 1.00e-04)\n", - "Best objective 6.529283067918e+04, best bound 6.529283067918e+04, gap 0.0000%\n" - ] - } - ], - "source": [ - "m.Params.NonConvex = 2\n", - "m.optimize()" - ] - }, - { - "cell_type": "code", - "execution_count": 188, - "id": "ec56f1da", - "metadata": {}, - "outputs": [ + "cell_type": "code", + "execution_count": 6, + "id": "3095da93", + "metadata": { + "id": "3095da93" + }, + "outputs": [], + "source": [ + "X = df[[\"p[1]\",\"p[2]\"]]\n", + "y = df[\"n[1]\"]\n", + "# Split the data for training and testing\n", + "X_train, X_test, y_train, y_test = train_test_split(\n", + " X, y, train_size=0.75, random_state=1\n", + ")" + ] + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "Optimal price for the two categories:\n", - " 378.77 300.0\n", - "\n", - "Optimal number of space assigned to the two categories:\n", - " 67 133\n", - "\n", - "Total revenue:\n", - " 65292.83\n" - ] - } - ], - "source": [ - "print(\"\\nOptimal price for the two categories:\\n\",round(p[1].X,2),round(p[2].X,2))\n", - "print(\"\\nOptimal number of space assigned to the two categories:\\n\",round(n[1].X), round(n[2].X))\n", - "print(\"\\nTotal revenue:\\n\",round(revenue.getValue(),2))" - ] - }, - { - "cell_type": "markdown", - "id": "5e9e2bbd", - "metadata": {}, - "source": [ - "## Build an interactive model" - ] - }, - { - "cell_type": "markdown", - "id": "e9e00776", - "metadata": {}, - "source": [ - "The penalty portion of the objective can have a big impact on the solution. We can combine the parts of the model into a function to give us the ability to see what different values of $\\lambda$ do to the solution. \n" - ] - }, - { - "cell_type": "code", - "execution_count": 111, - "id": "c0c84944", - "metadata": {}, - "outputs": [], - "source": [ - "from ipywidgets import interact, interactive, fixed, interact_manual\n", - "import ipywidgets as widgets\n", - "\n", - "def solve(x = 0):\n", - "\n", - " model = gp.Model(\"Price optimization\")\n", - " model.Params.OutputFlag = 0\n", - "\n", - " l = x\n", - " N = 200\n", - " price_lb = {1:300, 2:100}\n", - " price_ub = {1:400, 2:300}\n", - " min_items = {1:50,2:50}\n", - " price_bounds = {1:[300,400], 2:[100,300]}\n", - " p = model.addVars(products, lb = price_lb, ub = price_ub, name=\"p\") \n", - " n = model.addVars(products, lb = min_items, name=\"n\") \n", - " \n", - " model.addConstrs(n[c] >= min_items[c] for c in products) \n", - " model.addConstr(p[1] == [300,400]) \n", - " model.addConstr(p[2] == [100,300])\n", - "\n", - " revenue = gp.quicksum(p[c]*n[c] for c in products) \n", - " penalty = l*(p[1]**2+p[1]**2) \n", - " model.setObjective(revenue - penalty, sense = GRB.MAXIMIZE)\n", - "\n", - " model.addConstr(n.sum() == N)\n", - "\n", - " model.addConstr(p[1]-p[2] == [50,100])\n", - "\n", - " m_feats = pd.DataFrame({\"p[1]\":[p[1]],\"p[2]\":[p[2]]})\n", - "\n", - " pred_constr = add_predictor_constr(m, xgb_regressor, m_feats, n[1])\n", - "\n", - " model.Params.NonConvex = 2\n", - " model.optimize()\n", - " print(\"\\nOptimal price for the two categories:\\n\",round(p[1].X,2),round(p[2].X,2))\n", - " print(\"\\nOptimal number of space assigned to the two categories:\\n\",round(n[1].X), round(n[2].X))\n", - " print(\"\\nTotal revenue:\\n\",round(revenue.getValue(),2))\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e3b012c2", - "metadata": {}, - "outputs": [], - "source": [ - "upper_bound_for_lamda = 0.25\n", - "print(upper_bound_for_lamda)\n", - "print(\"Select a value for regularization parameter lambda:\\n\")\n", - "interact(solve, x=(0,upper_bound_for_lamda,0.01))" - ] - }, - { - "cell_type": "markdown", - "id": "e0888fb6", - "metadata": { - "id": "e0888fb6" - }, - "source": [ - "Copyright © 2024 Gurobi Optimization, LLC" - ] - } - ], - "metadata": { - "colab": { - "provenance": [] - }, - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.7" - }, - "widgets": { - "application/vnd.jupyter.widget-state+json": { - "04c449bbca144b7db6821cf019ebb48b": { - "model_module": "@jupyter-widgets/output", - "model_module_version": "1.0.0", - "model_name": "OutputModel", - "state": { - "_dom_classes": [], - "_model_module": "@jupyter-widgets/output", - "_model_module_version": "1.0.0", - "_model_name": "OutputModel", - "_view_count": null, - "_view_module": "@jupyter-widgets/output", - "_view_module_version": "1.0.0", - "_view_name": "OutputView", - "layout": "IPY_MODEL_c18fa29a67c643e79a1f0dffb1b591b6", - "msg_id": "", + "cell_type": "markdown", + "id": "f8d713c2", + "metadata": { + "id": "f8d713c2" + }, + "source": [ + "First we'll start with a linear regression model." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "33f4eaec", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "33f4eaec", + "outputId": "d02c9203-4bef-4a4e-cb16-bb948a8fca00" + }, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model for pipe:\n", - "180 variables\n", - "21 constraints\n", - "559 general constraints\n", - "Input has shape (1, 2)\n", - "Output has shape (1, 1)\n", - "\n", - "Pipeline has 1 steps:\n", - "\n", - "--------------------------------------------------------------------------------\n", - "Step Output Shape Variables Constraints \n", - " Linear Quadratic General\n", - "================================================================================\n", - "gbtree_reg (1, 1) 180 21 0 559\n", - "\n", - "--------------------------------------------------------------------------------\n", - "Set parameter NonConvex to value 2\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Gurobi Optimizer version 10.0.2 build v10.0.2rc0 (linux64)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Optimize a model with 30 rows, 187 columns and 193 nonzeros\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model fingerprint: 0xb29f40d6\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model has 3 quadratic objective terms\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model has 2 quadratic constraints\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model has 559 general constraints\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Variable types: 27 continuous, 160 integer (160 binary)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Coefficient statistics:\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " Matrix range [1e-01, 1e+00]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " QMatrix range [1e+00, 1e+00]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " QLMatrix range [1e+00, 1e+00]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " Objective range [0e+00, 0e+00]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " QObjective range [3e-01, 2e+00]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " Bounds range [6e-02, 1e+00]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " RHS range [5e-01, 4e+02]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " GenCon rhs range [5e-04, 4e+02]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " GenCon coe range [1e+00, 1e+00]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Presolve added 113 rows and 0 columns\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Presolve time: 0.06s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Presolved: 156 rows, 191 columns, 1024 nonzeros\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Presolved model has 1 quadratic constraint(s)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Presolved model has 4 bilinear constraint(s)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Variable types: 31 continuous, 160 integer (160 binary)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Root relaxation: objective 4.668508e+04, 188 iterations, 0.00 seconds (0.00 work units)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " Nodes | Current Node | Objective Bounds | Work\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 46685.0841 0 43 - 46685.0841 - - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "H 0 0 25589.338929 46685.0841 82.4% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 40215.3725 0 21 25589.3389 40215.3725 57.2% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "H 0 0 37190.120990 40215.3725 8.13% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 40207.4274 0 19 37190.1210 40207.4274 8.11% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 39746.2851 0 29 37190.1210 39746.2851 6.87% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "H 0 0 37399.288369 39746.2851 6.28% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 39724.9188 0 27 37399.2884 39724.9188 6.22% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 39531.4442 0 37 37399.2884 39531.4442 5.70% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 39516.6123 0 50 37399.2884 39516.6123 5.66% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 39436.2486 0 43 37399.2884 39436.2486 5.45% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 39433.6216 0 50 37399.2884 39433.6216 5.44% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 39423.4882 0 39 37399.2884 39423.4882 5.41% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 39405.9128 0 53 37399.2884 39405.9128 5.37% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 39401.2897 0 50 37399.2884 39401.2897 5.35% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 0 39033.2506 0 21 37399.2884 39033.2506 4.37% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0 1 39033.0955 0 21 37399.2884 39033.0955 4.37% - 0s\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Cutting planes:\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " Learned: 13\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " Cover: 12\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " Clique: 33\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " MIR: 6\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " RLT: 3\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " Relax-and-lift: 15\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Explored 22 nodes (590 simplex iterations) in 0.43 seconds (0.10 work units)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Thread count was 2 (of 2 available processors)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Solution count 3: 37399.3 37190.1 25589.3 \n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Optimal solution found (tolerance 1.00e-04)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Best objective 3.739928836871e+04, best bound 3.739928836871e+04, gap 0.0000%\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "Optimal price for the two categories:\n", - " 325.04908752441406 275.04908752441406\n", - "\n", - "Optimal number of seats assigned to the two categories:\n", - " 150.00000000000006 49.99999999999994\n" - ] - } + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "(array([0.80665368, 0.81004987, 0.7941077 , 0.79780172, 0.79495142]),\n", + " array([0.77650521, 0.74227784, 0.82105606, 0.81102818, 0.82166193]))" + ] + }, + "metadata": {}, + "execution_count": 7 + } + ], + "source": [ + "linear_regressor = make_pipeline(LinearRegression())\n", + "linear_regressor.fit(X_train, y_train)\n", + "linear_regression_validation = cross_validate(linear_regressor, X_train, y_train, cv=5, return_train_score=True, return_estimator=True)\n", + "\n", + "linear_regression_validation['train_score'],linear_regression_validation['test_score']" + ] + }, + { + "cell_type": "markdown", + "id": "cc9dd59b", + "metadata": { + "id": "cc9dd59b" + }, + "source": [ + "Let's try a gradient boosting model as well." ] - } }, - "0653e902f546425e9506d8340450dbd2": { - "model_module": "@jupyter-widgets/controls", - "model_module_version": "1.5.0", - "model_name": "VBoxModel", - "state": { - "_dom_classes": [ - "widget-interact" + { + "cell_type": "code", + "execution_count": 8, + "id": "9b8eb814", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "9b8eb814", + "outputId": "ca16f36f-2399-4e29-c3b3-9646cdb22a38" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "(array([0.71068886, 0.71262322, 0.70264476, 0.70564959, 0.7008412 ]),\n", + " array([0.65706085, 0.69348053, 0.66690978, 0.67664668, 0.67992755]))" + ] + }, + "metadata": {}, + "execution_count": 8 + } ], - "_model_module": "@jupyter-widgets/controls", - "_model_module_version": "1.5.0", - "_model_name": "VBoxModel", - "_view_count": null, - "_view_module": "@jupyter-widgets/controls", - "_view_module_version": "1.5.0", - "_view_name": "VBoxView", - "box_style": "", - "children": [ - "IPY_MODEL_e212c9758a0e4460950cb22f911671e1", - "IPY_MODEL_04c449bbca144b7db6821cf019ebb48b" + "source": [ + "from sklearn.ensemble import GradientBoostingRegressor\n", + "xgb_regressor = make_pipeline(GradientBoostingRegressor(n_estimators=10))\n", + "xgb_regressor.fit(X_train, y_train)\n", + "xgb_regressor_validation = cross_validate(xgb_regressor, X_train, y_train, cv=5, return_train_score=True, return_estimator=True)\n", + "\n", + "xgb_regressor_validation['train_score'], xgb_regressor_validation['test_score']\n" + ] + }, + { + "cell_type": "markdown", + "id": "mTquaGiJF2pO", + "metadata": { + "id": "mTquaGiJF2pO" + }, + "source": [ + "## Price optimization model with competing products\n", + "\n", + "Our problem is to:\n", + "1.\tDetermine the number of each category of product to make available given the overall restriction of what we can offer.\n", + "2.\tWe are also instructed to make sure there are a minimum number of each category made available as well as a minimum and maximum price for each category.\n", + "3.\tLastly, the product categories should be decreasing in price, meaning Category 1 should be the most expensive, and so on. Specifically, we must make sure there is at least a $50 gap between categories, but no more than $100.\n", + "\n", + "With the predictive part in place, it's time to build the optimization model. The model is formulated (i.e. the mathematical representation) for an unspecified number of categories, but the code will reflect that we have two categories of products in this problem. We start by setting some parameter values (not to be confused with ML hyperparameters) and initialize the optimization model.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "95b966b2", + "metadata": { + "id": "95b966b2" + }, + "source": [ + "### Initialize model and set input parameters\n", + "- $C$: Number of product categories\n", + "- $N$: Total amount of space available\n", + "- $\\lambda$: Price control parameter\n", + "\n", + "Here is the first mention of a price control parameter. It is fairly common in optimization modeling to add penalty terms to try and prevent undesirable outcomes. This is akin to using penalty terms in machine learning and applied statistics to prevent overfitting, with [Lasso](https://en.wikipedia.org/wiki/Lasso_(statistics)) and [Ridge](https://en.wikipedia.org/wiki/Ridge_regression) regression as a couple of common examples." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "67db9804", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "67db9804", + "outputId": "2a1edc1c-0b23-4649-88da-53d34a5f2b03" + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Restricted license - for non-production use only - expires 2025-11-24\n" + ] + } ], - "layout": "IPY_MODEL_c2a976389eb046728c1cf683b4920a85" - } + "source": [ + "#### Initialize the model\n", + "m = gp.Model(\"price optimization\")\n", + "\n", + "products = [1,2] #### Category 1 and Category 2\n", + "N = 200 #### limit on available space\n", + "l = 0 #### price control, we'll start this at 0" + ] + }, + { + "cell_type": "markdown", + "id": "ac5a875b", + "metadata": { + "id": "ac5a875b" + }, + "source": [ + "### Decision variables\n", + "- $p_c$: price per item in category $c = 1,2,\\dots, C$\n", + "- $n_c$: number of items allocated to category, predicted using features $p_c$, $c = 1,2,\\dots, C$" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "bbe94a3c", + "metadata": { + "id": "bbe94a3c" + }, + "outputs": [], + "source": [ + "p = m.addVars(products, name=\"p\") #### price decision variables\n", + "n = m.addVars(products, name=\"n\") #### decision variable for number of items in each category" + ] + }, + { + "cell_type": "markdown", + "id": "721614e0", + "metadata": { + "id": "721614e0" + }, + "source": [ + "### Constraints" + ] + }, + { + "cell_type": "markdown", + "id": "8403f941", + "metadata": { + "id": "8403f941" + }, + "source": [ + "We need to have a minimum number of each category available.\n", + "\\begin{align*}\n", + "n_c \\ge l_{n_c}\n", + "\\end{align*}\n", + "\n", + "We also set lower and upper bounds on the prices.\n", + "\\begin{align*}\n", + "l_{n_c} \\le p_c \\le u_{p_c}\n", + "\\end{align*}" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "edec4432", + "metadata": { + "id": "edec4432" + }, + "outputs": [], + "source": [ + "min_items = {1:50,2:50}\n", + "price_bounds = {1:[300,400], 2:[100,300]}\n", + "m.addConstrs(n[c] >= min_items[c] for c in products) #### we could hardcode 50 instead of min_items, but this is more flexible\n", + "m.addConstr(p[1] == [300,400]) #### this is a shorthand way to code 300 <= p[1] <= 400\n", + "m.addConstr(p[2] == [100,300]);" + ] + }, + { + "cell_type": "markdown", + "id": "3bef9e0d", + "metadata": { + "id": "3bef9e0d" + }, + "source": [ + "Another note: each of the above constraints can be addressed when defining the decision variables. Here is an example for the decision variable $n$." + ] }, - "39b9c940f508451b9292f29639f5d128": { - "model_module": "@jupyter-widgets/controls", - "model_module_version": "1.5.0", - "model_name": "SliderStyleModel", - "state": { - "_model_module": "@jupyter-widgets/controls", - "_model_module_version": "1.5.0", - "_model_name": "SliderStyleModel", - "_view_count": null, - "_view_module": "@jupyter-widgets/base", - "_view_module_version": "1.2.0", - "_view_name": "StyleView", - "description_width": "", - "handle_color": null - } + { + "cell_type": "code", + "execution_count": 12, + "id": "e2065d83", + "metadata": { + "id": "e2065d83" + }, + "outputs": [], + "source": [ + "#price_lb = {1:300, 2:100}\n", + "#price_ub = {1:400, 2:300}\n", + "#p = m.addVars(products, lb = price_lb, ub = price_ub, name=\"p\") #### each price is now bounded\n", + "#n = m.addVars(products, lb = min_items, name=\"n\")" + ] }, - "c18fa29a67c643e79a1f0dffb1b591b6": { - "model_module": "@jupyter-widgets/base", - "model_module_version": "1.2.0", - "model_name": "LayoutModel", - "state": { - "_model_module": "@jupyter-widgets/base", - "_model_module_version": "1.2.0", - "_model_name": "LayoutModel", - "_view_count": null, - "_view_module": "@jupyter-widgets/base", - "_view_module_version": "1.2.0", - "_view_name": "LayoutView", - "align_content": null, - "align_items": null, - "align_self": null, - "border": null, - "bottom": null, - "display": null, - "flex": null, - "flex_flow": null, - "grid_area": null, - "grid_auto_columns": null, - "grid_auto_flow": null, - "grid_auto_rows": null, - "grid_column": null, - "grid_gap": null, - "grid_row": null, - "grid_template_areas": null, - "grid_template_columns": null, - "grid_template_rows": null, - "height": null, - "justify_content": null, - "justify_items": null, - "left": null, - "margin": null, - "max_height": null, - "max_width": null, - "min_height": null, - "min_width": null, - "object_fit": null, - "object_position": null, - "order": null, - "overflow": null, - "overflow_x": null, - "overflow_y": null, - "padding": null, - "right": null, - "top": null, - "visibility": null, - "width": null - } + { + "cell_type": "markdown", + "id": "0578e271", + "metadata": { + "id": "0578e271" + }, + "source": [ + "In general, the number of items allocated must equal the total available space.\n", + "\\begin{equation*}\n", + "n_1 + n_2 + \\dots + n_C = \\sum_{c}n_c = N \\\\\n", + "\\end{equation*}\n", + "Note that this, along with the constraint on the minimum number available means we don't have to specify an upper bound for each $n_c$." + ] }, - "c2a976389eb046728c1cf683b4920a85": { - "model_module": "@jupyter-widgets/base", - "model_module_version": "1.2.0", - "model_name": "LayoutModel", - "state": { - "_model_module": "@jupyter-widgets/base", - "_model_module_version": "1.2.0", - "_model_name": "LayoutModel", - "_view_count": null, - "_view_module": "@jupyter-widgets/base", - "_view_module_version": "1.2.0", - "_view_name": "LayoutView", - "align_content": null, - "align_items": null, - "align_self": null, - "border": null, - "bottom": null, - "display": null, - "flex": null, - "flex_flow": null, - "grid_area": null, - "grid_auto_columns": null, - "grid_auto_flow": null, - "grid_auto_rows": null, - "grid_column": null, - "grid_gap": null, - "grid_row": null, - "grid_template_areas": null, - "grid_template_columns": null, - "grid_template_rows": null, - "height": null, - "justify_content": null, - "justify_items": null, - "left": null, - "margin": null, - "max_height": null, - "max_width": null, - "min_height": null, - "min_width": null, - "object_fit": null, - "object_position": null, - "order": null, - "overflow": null, - "overflow_x": null, - "overflow_y": null, - "padding": null, - "right": null, - "top": null, - "visibility": null, - "width": null - } + { + "cell_type": "code", + "execution_count": 13, + "id": "140e739f", + "metadata": { + "id": "140e739f" + }, + "outputs": [], + "source": [ + "m.addConstr(n.sum() == N); #### remember we set N = 200 earlier" + ] }, - "e212c9758a0e4460950cb22f911671e1": { - "model_module": "@jupyter-widgets/controls", - "model_module_version": "1.5.0", - "model_name": "FloatSliderModel", - "state": { - "_dom_classes": [], - "_model_module": "@jupyter-widgets/controls", - "_model_module_version": "1.5.0", - "_model_name": "FloatSliderModel", - "_view_count": null, - "_view_module": "@jupyter-widgets/controls", - "_view_module_version": "1.5.0", - "_view_name": "FloatSliderView", - "continuous_update": true, - "description": "x", - "description_tooltip": null, - "disabled": false, - "layout": "IPY_MODEL_e58bd5527b16464da224d7e8125e2065", - "max": 0.08, - "min": 0, - "orientation": "horizontal", - "readout": true, - "readout_format": ".2f", - "step": 0.01, - "style": "IPY_MODEL_39b9c940f508451b9292f29639f5d128", - "value": 0.08 - } + { + "cell_type": "markdown", + "id": "b17a1e78", + "metadata": { + "id": "b17a1e78" + }, + "source": [ + "The last set of constraints are for price ordering. This requires the subsequent category to be cost between $50 and $100 less than the previous.\n", + "\\begin{equation*}\n", + "50 \\le p_c - p_{c+1} \\le 100\n", + "\\end{equation*}" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "1428be22", + "metadata": { + "id": "1428be22" + }, + "outputs": [], + "source": [ + "m.addConstr(p[1]-p[2] == [50,100]);" + ] }, - "e58bd5527b16464da224d7e8125e2065": { - "model_module": "@jupyter-widgets/base", - "model_module_version": "1.2.0", - "model_name": "LayoutModel", - "state": { - "_model_module": "@jupyter-widgets/base", - "_model_module_version": "1.2.0", - "_model_name": "LayoutModel", - "_view_count": null, - "_view_module": "@jupyter-widgets/base", - "_view_module_version": "1.2.0", - "_view_name": "LayoutView", - "align_content": null, - "align_items": null, - "align_self": null, - "border": null, - "bottom": null, - "display": null, - "flex": null, - "flex_flow": null, - "grid_area": null, - "grid_auto_columns": null, - "grid_auto_flow": null, - "grid_auto_rows": null, - "grid_column": null, - "grid_gap": null, - "grid_row": null, - "grid_template_areas": null, - "grid_template_columns": null, - "grid_template_rows": null, - "height": null, - "justify_content": null, - "justify_items": null, - "left": null, - "margin": null, - "max_height": null, - "max_width": null, - "min_height": null, - "min_width": null, - "object_fit": null, - "object_position": null, - "order": null, - "overflow": null, - "overflow_x": null, - "overflow_y": null, - "padding": null, - "right": null, - "top": null, - "visibility": null, - "width": null - } + { + "cell_type": "markdown", + "id": "51ee2204", + "metadata": { + "id": "51ee2204" + }, + "source": [ + "### Objective function" + ] + }, + { + "cell_type": "markdown", + "id": "6eb69e54", + "metadata": { + "id": "6eb69e54" + }, + "source": [ + "We want to maximize total revenue with the portion of total revenue coming from category $c$ being $p_cn_c$. This makes the total revenue $\\sum_{c} p_c n_c$. That is the first part of the objective. Earlier a price control parameter was introduced which is the second part of the objective. The lambda parameter captures the trade-off between the revenue and price-control pieces. This term penalizes the model from setting too high of prices since doing so could lose sales. Our model assumed we'll sell all of the items so having this penalty term can make this assumption more realistic. For reference, [here is a good source](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4565407).\n", + "\n", + "This term will be defined as $λ (\\sum_{c} p_c^2)$ for this problem. So, the complete objective is:\n", + "\\begin{equation*}\n", + "\\textrm{maximize} \\sum_{c} p_c n_c - λ (\\sum_{c} p_c^2)\n", + "\\end{equation*}" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "a2cfcabc", + "metadata": { + "id": "a2cfcabc" + }, + "outputs": [], + "source": [ + "revenue = gp.quicksum(p[c]*n[c] for c in products) #### you could also use the more simple p.prod(n)\n", + "penalty = l*(p[1]**2+p[1]**2) #### we used l as the lambda parameter earlier\n", + "m.setObjective(revenue - penalty, sense = GRB.MAXIMIZE)" + ] + }, + { + "cell_type": "markdown", + "id": "c707af60", + "metadata": { + "id": "c707af60" + }, + "source": [ + "### Integrate the ML model\n" + ] + }, + { + "cell_type": "markdown", + "id": "38d41039", + "metadata": { + "id": "38d41039" + }, + "source": [ + "Right now, if we were to run the optimization, the solution would be to set the price for Category 1 to $400, Category 2 to $300, and sell 150 and 50 of each item, respectively. That's because we have yet to add in the relationship between price and demand that was derived from the ML model. To integrate the machine learning model into the optimization model, we'll use the Gurobi Machine Learning package. The magic happens using `add_predictor_constr` function." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "b113eda2", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "b113eda2", + "outputId": "584d9f0e-ea69-41c9-e481-a45f0bc36dd6" + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Requirement already satisfied: gurobi-machinelearning in /usr/local/lib/python3.10/dist-packages (1.4.0)\n", + "Requirement already satisfied: numpy>=1.22.0 in /usr/local/lib/python3.10/dist-packages (from gurobi-machinelearning) (1.25.2)\n", + "Requirement already satisfied: gurobipy>=10.0.0 in /usr/local/lib/python3.10/dist-packages (from gurobi-machinelearning) (11.0.0)\n", + "Requirement already satisfied: scipy>=1.9.3 in /usr/local/lib/python3.10/dist-packages (from gurobi-machinelearning) (1.11.4)\n" + ] + } + ], + "source": [ + "#### install the package and load the required function\n", + "%pip install gurobi-machinelearning\n", + "from gurobi_ml import add_predictor_constr" + ] + }, + { + "cell_type": "markdown", + "id": "119a8f7e", + "metadata": { + "id": "119a8f7e" + }, + "source": [ + "This additional package is useful when we have **decision variables** that are also **features** of a machine learning model. First, we need a data frame that contains these decision variables. It is important to make sure the indices of the data frame have the **same name** as the training data for the machine learning model." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "0b1ac8cc", + "metadata": { + "id": "0b1ac8cc" + }, + "outputs": [], + "source": [ + "m_feats = pd.DataFrame({\"p[1]\":[p[1]],\"p[2]\":[p[2]]})" + ] + }, + { + "cell_type": "markdown", + "id": "3263a396", + "metadata": { + "id": "3263a396" + }, + "source": [ + "Adding the predictive model to the optimization model requires specifying the model we want to use `(m)`, regression object `(xgb_regressor)`, feature data frame `(m_feats)`, and the output decision variable `(n[1])`. Remember `n[2]` is **NOT** the output of the regression. We can then print the number of variables and constraints added to the model using `print_stats`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "1fe6c59f", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "1fe6c59f", + "outputId": "5aa94f80-e974-4e5b-de19-7c82525e5f75" + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Model for pipe:\n", + "90 variables\n", + "11 constraints\n", + "246 general constraints\n", + "Input has shape (1, 2)\n", + "Output has shape (1, 1)\n", + "\n", + "Pipeline has 1 steps:\n", + "\n", + "--------------------------------------------------------------------------------\n", + "Step Output Shape Variables Constraints \n", + " Linear Quadratic General\n", + "================================================================================\n", + "gbtree_reg (1, 1) 90 11 0 246\n", + "\n", + "--------------------------------------------------------------------------------\n" + ] + } + ], + "source": [ + "pred_constr = add_predictor_constr(m, xgb_regressor, m_feats, n[1])\n", + "pred_constr.print_stats()" + ] + }, + { + "cell_type": "markdown", + "id": "53ce50f6", + "metadata": { + "id": "53ce50f6" + }, + "source": [ + "### Solve the optimization and get the solution" + ] + }, + { + "cell_type": "markdown", + "id": "ec80341d", + "metadata": { + "id": "ec80341d" + }, + "source": [ + "Since this is a quadratic, non-convex problem we set the `NonConvex` parameter to 2. See the [documentation](https://www.gurobi.com/documentation/current/refman/nonconvex.html) for more information. We'll also print out the optimal solution." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "d04d5e5b", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "d04d5e5b", + "outputId": "7cbae0da-917f-45ec-9fe1-3af04c942d08" + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Set parameter NonConvex to value 2\n", + "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.3 LTS\")\n", + "\n", + "CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n", + "Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n", + "\n", + "Optimize a model with 17 rows, 97 columns and 102 nonzeros\n", + "Model fingerprint: 0x2255329b\n", + "Model has 2 quadratic objective terms\n", + "Model has 246 general constraints\n", + "Variable types: 17 continuous, 80 integer (80 binary)\n", + "Coefficient statistics:\n", + " Matrix range [1e-01, 1e+00]\n", + " Objective range [0e+00, 0e+00]\n", + " QObjective range [2e+00, 2e+00]\n", + " Bounds range [1e+00, 2e+02]\n", + " RHS range [1e+00, 4e+02]\n", + " GenCon rhs range [3e-01, 4e+02]\n", + " GenCon coe range [1e+00, 1e+00]\n", + "Presolve added 44 rows and 0 columns\n", + "Presolve removed 0 rows and 15 columns\n", + "Presolve time: 0.10s\n", + "Presolved: 66 rows, 85 columns, 383 nonzeros\n", + "Presolved model has 2 bilinear constraint(s)\n", + "\n", + "Solving non-convex MIQCP\n", + "\n", + "Variable types: 18 continuous, 67 integer (67 binary)\n", + "Found heuristic solution: objective 57353.074762\n", + "\n", + "Root relaxation: objective 6.842235e+04, 81 iterations, 0.00 seconds (0.00 work units)\n", + "\n", + " Nodes | Current Node | Objective Bounds | Work\n", + " Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n", + "\n", + " 0 0 68422.3454 0 21 57353.0748 68422.3454 19.3% - 0s\n", + " 0 0 68380.4264 0 21 57353.0748 68380.4264 19.2% - 0s\n", + "H 0 0 65051.531333 68380.4264 5.12% - 0s\n", + "H 0 0 65717.039440 68380.4264 4.05% - 0s\n", + " 0 0 67521.7136 0 30 65717.0394 67521.7136 2.75% - 0s\n", + " 0 0 67440.6730 0 34 65717.0394 67440.6730 2.62% - 0s\n", + " 0 0 67415.9059 0 33 65717.0394 67415.9059 2.59% - 0s\n", + " 0 0 67407.6626 0 32 65717.0394 67407.6626 2.57% - 0s\n", + " 0 0 67397.1289 0 32 65717.0394 67397.1289 2.56% - 0s\n", + " 0 0 67367.1673 0 26 65717.0394 67367.1673 2.51% - 0s\n", + "H 0 0 65919.660132 67367.1673 2.20% - 0s\n", + " 0 0 67325.6782 0 24 65919.6601 67325.6782 2.13% - 0s\n", + " 0 0 67027.1365 0 27 65919.6601 67027.1365 1.68% - 0s\n", + " 0 0 66962.1686 0 21 65919.6601 66962.1686 1.58% - 0s\n", + " 0 0 66754.1535 0 21 65919.6601 66754.1535 1.27% - 0s\n", + " 0 0 66724.4366 0 21 65919.6601 66724.4366 1.22% - 0s\n", + " 0 0 66357.6817 0 22 65919.6601 66357.6817 0.66% - 0s\n", + " 0 0 66314.3792 0 17 65919.6601 66314.3792 0.60% - 0s\n", + " 0 0 66212.7588 0 15 65919.6601 66212.7588 0.44% - 0s\n", + " 0 0 66096.5119 0 21 65919.6601 66096.5119 0.27% - 0s\n", + " 0 0 66096.5119 0 22 65919.6601 66096.5119 0.27% - 0s\n", + " 0 0 66096.5119 0 21 65919.6601 66096.5119 0.27% - 0s\n", + " 0 0 66096.5119 0 23 65919.6601 66096.5119 0.27% - 0s\n", + " 0 0 66092.0244 0 21 65919.6601 66092.0244 0.26% - 0s\n", + " 0 0 66092.0244 0 23 65919.6601 66092.0244 0.26% - 0s\n", + " 0 0 66092.0244 0 22 65919.6601 66092.0244 0.26% - 0s\n", + " 0 0 66092.0244 0 21 65919.6601 66092.0244 0.26% - 0s\n", + " 0 0 66092.0244 0 21 65919.6601 66092.0244 0.26% - 0s\n", + " 0 0 66092.0244 0 23 65919.6601 66092.0244 0.26% - 0s\n", + " 0 0 66092.0244 0 19 65919.6601 66092.0244 0.26% - 0s\n", + " 0 0 66087.8039 0 21 65919.6601 66087.8039 0.26% - 0s\n", + " 0 0 66078.8922 0 17 65919.6601 66078.8922 0.24% - 0s\n", + " 0 0 66074.6109 0 20 65919.6601 66074.6109 0.24% - 0s\n", + " 0 0 65978.7306 0 15 65919.6601 65978.7306 0.09% - 0s\n", + " 0 0 65955.6109 0 23 65919.6601 65955.6109 0.05% - 0s\n", + " 0 0 65941.1212 0 23 65919.6601 65941.1212 0.03% - 0s\n", + " 0 0 65933.9387 0 23 65919.6601 65933.9387 0.02% - 0s\n", + " 0 0 65933.9387 0 22 65919.6601 65933.9387 0.02% - 0s\n", + " 0 0 65931.8301 0 1 65919.6601 65931.8301 0.02% - 0s\n", + "H 0 0 65931.830113 65931.8301 0.00% - 0s\n", + " 0 0 65931.8301 0 1 65931.8301 65931.8301 0.00% - 0s\n", + "\n", + "Cutting planes:\n", + " Cover: 3\n", + " Implied bound: 2\n", + " Clique: 6\n", + " MIR: 2\n", + " GUB cover: 1\n", + " Relax-and-lift: 4\n", + "\n", + "Explored 1 nodes (377 simplex iterations) in 0.62 seconds (0.05 work units)\n", + "Thread count was 2 (of 2 available processors)\n", + "\n", + "Solution count 5: 65931.8 65919.7 65717 ... 57353.1\n", + "\n", + "Optimal solution found (tolerance 1.00e-04)\n", + "Best objective 6.593183011274e+04, best bound 6.593183011274e+04, gap 0.0000%\n" + ] + } + ], + "source": [ + "m.Params.NonConvex = 2\n", + "m.optimize()" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "ec56f1da", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "ec56f1da", + "outputId": "231a9f5b-2d76-48b2-899b-8f9db2232b94" + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "\n", + "Optimal price for the two categories:\n", + " 364.97 300.0\n", + "\n", + "Optimal number of space assigned to the two categories:\n", + " 91 109\n", + "\n", + "Total revenue:\n", + " 65931.83\n" + ] + } + ], + "source": [ + "print(\"\\nOptimal price for the two categories:\\n\",round(p[1].X,2),round(p[2].X,2))\n", + "print(\"\\nOptimal number of space assigned to the two categories:\\n\",round(n[1].X), round(n[2].X))\n", + "print(\"\\nTotal revenue:\\n\",round(revenue.getValue(),2))" + ] + }, + { + "cell_type": "markdown", + "id": "5e9e2bbd", + "metadata": { + "id": "5e9e2bbd" + }, + "source": [ + "## Follow-up questions\n", + "At this point the λ value is set to zero, meaning we are not penalizing the price. One way to see the sensitivity of the prices and number of items set to each category is solving the model for multiple values of λ. This can be done fairly easy in `gurobipy`. Check out the documentation for how Gurobi can handle [multiple scenarios]('https://www.gurobi.com/documentation/current/refman/multiple_scenarios.html').\n", + "\n", + "We also only used the trained `xgb_regressor` for the optimization to show the number of variables and constraints added to the model. Test the `linear_regressor` that was trained as well. You can do that by re-running the cells or by adding to the code to [remove the previous predictor constraints](https://gurobi-machinelearning.readthedocs.io/en/stable/auto_generated/gurobi_ml.sklearn.linear_regression.LinearRegressionConstr.html#gurobi_ml.sklearn.linear_regression.LinearRegressionConstr.remove) and adding new ones." + ] + }, + { + "cell_type": "markdown", + "id": "e9e00776", + "metadata": { + "id": "e9e00776" + }, + "source": [ + "The penalty portion of the objective can have a big impact on the solution. We can combine the parts of the model into a function to give us the ability to see what different values of $\\lambda$ do to the solution.\n" + ] + }, + { + "cell_type": "markdown", + "id": "e0888fb6", + "metadata": { + "id": "e0888fb6" + }, + "source": [ + "Copyright © 2024 Gurobi Optimization, LLC" + ] + } + ], + "metadata": { + "colab": { + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.6" } - } - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file