{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4ae5c2ae",
   "metadata": {},
   "source": [
    "# Factor Model as Objective\n",
    "\n",
    "In this notebook, we explore how to formulate a portfolio optimization problem where the risk estimator is not given explicitly but is instead represented by a factor model, as is common in US equity models [1]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ea2b25b7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-31T10:04:23.173377Z",
     "iopub.status.busy": "2025-01-31T10:04:23.173139Z",
     "iopub.status.idle": "2025-01-31T10:04:23.959695Z",
     "shell.execute_reply": "2025-01-31T10:04:23.959000Z"
    },
    "nbsphinx": "hidden"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: numpy in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (2.2.2)\r\n",
      "Requirement already satisfied: scipy in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (1.15.1)\r\n",
      "Requirement already satisfied: gurobipy in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (11.0.3)\r\n",
      "Requirement already satisfied: pandas in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (2.2.3)\r\n",
      "Requirement already satisfied: matplotlib in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (3.10.0)\r\n",
      "Requirement already satisfied: python-dateutil>=2.8.2 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from pandas) (2.9.0.post0)\r\n",
      "Requirement already satisfied: pytz>=2020.1 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from pandas) (2025.1)\r\n",
      "Requirement already satisfied: tzdata>=2022.7 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from pandas) (2025.1)\r\n",
      "Requirement already satisfied: contourpy>=1.0.1 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from matplotlib) (1.3.1)\r\n",
      "Requirement already satisfied: cycler>=0.10 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from matplotlib) (0.12.1)\r\n",
      "Requirement already satisfied: fonttools>=4.22.0 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from matplotlib) (4.55.8)\r\n",
      "Requirement already satisfied: kiwisolver>=1.3.1 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from matplotlib) (1.4.8)\r\n",
      "Requirement already satisfied: packaging>=20.0 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from matplotlib) (24.2)\r\n",
      "Requirement already satisfied: pillow>=8 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from matplotlib) (11.1.0)\r\n",
      "Requirement already satisfied: pyparsing>=2.3.1 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from matplotlib) (3.2.1)\r\n",
      "Requirement already satisfied: six>=1.5 in /opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages (from python-dateutil>=2.8.2->pandas) (1.17.0)\r\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "# Install dependencies\n",
    "%pip install numpy scipy gurobipy pandas matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8c44ae69",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-31T10:04:23.961699Z",
     "iopub.status.busy": "2025-01-31T10:04:23.961509Z",
     "iopub.status.idle": "2025-01-31T10:04:24.605076Z",
     "shell.execute_reply": "2025-01-31T10:04:24.604409Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import gurobipy as gp\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "65f2022b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-31T10:04:24.607459Z",
     "iopub.status.busy": "2025-01-31T10:04:24.607051Z",
     "iopub.status.idle": "2025-01-31T10:04:24.615490Z",
     "shell.execute_reply": "2025-01-31T10:04:24.614864Z"
    },
    "nbsphinx": "hidden"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter WLSAccessID\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter WLSSecret\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter LicenseID to value 2443533\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WLS license 2443533 - registered to Gurobi GmbH\n"
     ]
    }
   ],
   "source": [
    "# Hidden cell to avoid licensing messages\n",
    "# when docs are generated.\n",
    "with gp.Model():\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a651963d",
   "metadata": {},
   "source": [
    "## Background\n",
    "\n",
    "The *standard mean-variance (Markowitz) portfolio selection model* determines optimal investments that balance risk and expected return. In this notebook, we maximize *utility*, which is defined as a weighted linear function of return and risk that can be adjusted by varying the risk-aversion coefficient $\\gamma$.\n",
    "\n",
    "The standard formulation of this problem on $n$ assets reads\n",
    "\n",
    "\\begin{alignat*}{2}\n",
    "\\tag{1}\n",
    "\\max_{x \\in \\mathbf{R}^{n}} \\quad        & \\mu^\\top x - \\tfrac12 \\gamma\\ x^\\top\\Sigma x \\\\\n",
    "\\mbox{s.t.} \\quad & 1_n^\\top x = 1 \\\\\n",
    "& 0 \\le x \\le 1\\\\\n",
    "\\end{alignat*}\n",
    "\n",
    "where $\\mu \\in \\mathbf{R}^n$ is the return estimator, $\\Sigma \\in \\mathbf{R}^{n,n}$ is used as an estimator for the variance, and $1_n$ is the vector of $n$ ones. While it is certainly possible that the matrix $\\Sigma$ is given explicitly (e.g., as the covariance matrix of a time series), it is often expressed implicitly through a *factor model*.  In that case, the matrix $\\Sigma$ takes the form\n",
    "\n",
    "$$\n",
    "\\Sigma = X \\Sigma_0 X^T + D\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $X \\in \\mathbf{R}^{n,k}$ is the factor exposure matrix,\n",
    "- $\\Sigma_0 \\in \\mathbf{R}^{k,k}$ is the symmetric positive definite (SPD) factor covariance matrix, and\n",
    "- $D \\in \\mathbf{R}^{n,n}$ is the *diagonal* matrix of specific risks with $d_{ii} > 0$ for all $i$.\n",
    "\n",
    "Effectively, this splits the risk into two sources:  One that arises from common risk factors (macroeconomic conditions, market trends, etc.) and a specific risk that is uncorrelated among the assets.  The number of factors, $k$, is typically much smaller than the number $n$ of assets.\n",
    "\n",
    "The important observation from a computational point of view is the following: The factor model data is of size $nk + n$ (with $k \\ll n$), while the covariance matrix $\\Sigma$ is on the order of $n^2$. As we will see, optimization models using the factor model admit a smaller representation and typically offer improved computational performance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b20a6ff",
   "metadata": {},
   "source": [
    "## A synthetic factor model\n",
    "\n",
    "Leading industry factor models are commercial products that are not necessary for our demonstration purposes. We are more interested in a qualitative comparison of solver performance using *some* factor model than in the actual economic meaning.  For this reason, our first step will be to create a simple, synthetic multivariate factor model that we can use in the study that follows.\n",
    "\n",
    "The following function uses an uncorrelated factor covariance matrix $\\Sigma_0$ that is used for sampling a multivariate normal distribution on `num_factors` factors along `time_steps` sampling points for `num_assets` assets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "115db05b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-31T10:04:24.617646Z",
     "iopub.status.busy": "2025-01-31T10:04:24.617454Z",
     "iopub.status.idle": "2025-01-31T10:04:24.622199Z",
     "shell.execute_reply": "2025-01-31T10:04:24.621565Z"
    }
   },
   "outputs": [],
   "source": [
    "def factor_model(num_assets, num_factors, timesteps):\n",
    "    # Generate random factor model, risk is X * sigma0 * X.T + cov(u)\n",
    "    sigma0 = np.diag(1 + np.arange(num_factors) + np.random.rand(num_factors))\n",
    "    X = np.random.normal(size=(num_assets, num_factors))\n",
    "    alpha = np.random.normal(loc=1, size=(num_assets, 1))\n",
    "    u = np.random.multivariate_normal(\n",
    "        np.zeros(num_assets), np.eye(num_assets), timesteps\n",
    "    ).T\n",
    "    d = np.diag(np.cov(u))  # NOTE: This is the _diagonal_ of D!\n",
    "\n",
    "    # Time series in factor space\n",
    "    TS_factor = np.random.multivariate_normal(\n",
    "        np.zeros(num_factors), sigma0, timesteps\n",
    "    ).T\n",
    "\n",
    "    # Estimate mu from time series in full space\n",
    "    mu = np.mean(alpha + X @ TS_factor + u, axis=1)\n",
    "\n",
    "    return X, sigma0, d, mu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0169b6e2",
   "metadata": {},
   "source": [
    "We skip the details of this statistical procedure for constructing $X$, $\\Sigma_0$, $D$, and $\\mu$; more background and details can be found in [2, Sect. 6.6].  The essential point is that we now have a synthetic factor model that acts similarly to a commercial factor model regarding solver performance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "993b9b7b",
   "metadata": {},
   "source": [
    "## Taking advantage of the factor model structure\n",
    "\n",
    "Ignoring the constant factor $\\frac{1}{2}\\gamma$ for a moment, the objective function we want to maximize includes\n",
    "\n",
    "$$\n",
    "x^T \\Sigma x = x^T (X \\Sigma_0 X + D) x = x^T X \\Sigma_0 X^T x + x^T D x\n",
    "$$\n",
    "\n",
    "The first term, $x^T X \\Sigma_0 X^T x$, would result in a quadratic function having $(n+1)\\frac n2$ terms.  Since $\\Sigma_0$ is SPD, it admits a *Cholesky factorization* $\\Sigma_0 = LL^T$ where $L \\in \\mathbf{R}^{k,k}$ is a triangular matrix.  This allows us to rewrite the first term as\n",
    "\n",
    "$$\n",
    "x^T X \\Sigma_0 X^T x = x^T \\underbrace{X L}_{B} \\underbrace{L^T X^T}_{B^T} x = \\underbrace{x^T B}_{y^T} \\underbrace{B^T x}_{y} = y^T y\n",
    "$$\n",
    "\n",
    "where we have substituted $y = B^T x$ in the last step. Also, note that $B = X L \\in \\mathbf{R}^{n, k}$ comprises only $nk$ elements.\n",
    "\n",
    "The second term, $x^T D x$, comprises only $n$ terms (i.e., $\\sum_i d_{ii} x_{i}^2$) and can be used as is.  Putting this together, the standard optimization model (1) can be rewritten as\n",
    "\n",
    "\\begin{alignat}{2}\n",
    "\\tag{2}\n",
    "\\max_{x \\in \\mathbf{R}^{n}, y \\in \\mathbf{R}^k} \\quad        & \\mu^\\top x - \\tfrac12 \\gamma\\ y^T y + \\tfrac12 \\gamma\\ \\sum_{i=1}^n d_{ii} x_i^2\\\\\n",
    "\\mbox{s.t.} \\quad & 1_n^\\top x = 1 \\\\\n",
    "& y = B^T x.\\\\\n",
    "& 0 \\le x \\le 1\n",
    "\\end{alignat}\n",
    "\n",
    "Note that the $k$ variables $y$ do not have bound constraints.  Model (2) contains the much smaller matrix $B$ instead of $\\Sigma$ in model (1), at the expense of $k$ additional optimization variables.  Generally, form (2) is advantageous for solver performance; we shall compare both in the next step."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a1060ef",
   "metadata": {},
   "source": [
    "## Comparing the two optimization models\n",
    "\n",
    "In order to compare the solution times for models (1) and (2) above, we will define two auxiliary functions to build an optimization model in either form.  Let's start with a function for the traditional form (1):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8ab6fffd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-31T10:04:24.624108Z",
     "iopub.status.busy": "2025-01-31T10:04:24.623926Z",
     "iopub.status.idle": "2025-01-31T10:04:24.627524Z",
     "shell.execute_reply": "2025-01-31T10:04:24.626873Z"
    }
   },
   "outputs": [],
   "source": [
    "def build_sigma_model(model, gamma, sigma, mu):\n",
    "    x = model.addMVar(mu.shape, lb=0.0, ub=1.0)\n",
    "    model.addConstr(x.sum() == 1)\n",
    "    model.setObjective(x @ mu - x @ (gamma / 2.0 * sigma) @ x, gp.GRB.MAXIMIZE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44ed7bb2",
   "metadata": {},
   "source": [
    "The second function build the equivalent model (2):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ded8ee55",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-31T10:04:24.629410Z",
     "iopub.status.busy": "2025-01-31T10:04:24.629228Z",
     "iopub.status.idle": "2025-01-31T10:04:24.633262Z",
     "shell.execute_reply": "2025-01-31T10:04:24.632806Z"
    }
   },
   "outputs": [],
   "source": [
    "def build_factor_model(model, gamma, B, d, mu):\n",
    "    x = model.addMVar(B.shape[0], lb=0.0, ub=1.0)  # n variables\n",
    "    y = model.addMVar(B.shape[1], lb=-float(\"inf\"), ub=float(\"inf\"))  # k variables\n",
    "    model.addConstr(x.sum() == 1)\n",
    "    model.addConstr(B.T @ x - y == 0)\n",
    "\n",
    "    model.setObjective(\n",
    "        x @ mu - gamma / 2.0 * y @ y - ((gamma / 2.0 * d) * x**2).sum(),\n",
    "        gp.GRB.MAXIMIZE,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38aee780",
   "metadata": {},
   "source": [
    "We are now ready to run a small benchmark.  We will keep a fixed risk aversion coefficient $\\gamma$, and solve both models over a range of data with an increasing number of assets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d463c359",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-31T10:04:24.635077Z",
     "iopub.status.busy": "2025-01-31T10:04:24.634870Z",
     "iopub.status.idle": "2025-01-31T10:04:31.882815Z",
     "shell.execute_reply": "2025-01-31T10:04:31.882208Z"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(0xACAC)  # Fix seed for reproducibility\n",
    "\n",
    "num_factors = 72  # USE4 has 72 factors, too (see [1])\n",
    "timesteps = 700\n",
    "gamma = 0.025\n",
    "\n",
    "problem_dims = np.linspace(100, 750, 16, endpoint=True, dtype=int)\n",
    "time_sigma = []\n",
    "time_factor = []\n",
    "\n",
    "for num_assets in problem_dims:\n",
    "    X, sigma0, d, mu = factor_model(num_assets, num_factors, timesteps)\n",
    "\n",
    "    with gp.Model() as m:\n",
    "        m.Params.OutputFlag = 0\n",
    "        sigma = X @ sigma0 @ X.T + np.diag(d)\n",
    "        build_sigma_model(m, gamma, sigma, mu)\n",
    "        m.optimize()\n",
    "        time_sigma.append(m.RunTime)\n",
    "\n",
    "    with gp.Model() as m:\n",
    "        m.Params.OutputFlag = 0\n",
    "        L = np.linalg.cholesky(sigma0)\n",
    "        B = X @ L\n",
    "        build_factor_model(m, gamma, B, d, mu)\n",
    "        m.optimize()\n",
    "        time_factor.append(m.RunTime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "926124d9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-31T10:04:31.885628Z",
     "iopub.status.busy": "2025-01-31T10:04:31.885304Z",
     "iopub.status.idle": "2025-01-31T10:04:32.027019Z",
     "shell.execute_reply": "2025-01-31T10:04:32.026292Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYG9JREFUeJzt3Xl4E1Xbx/FvGrpQStmhpVQK4gKyiCCIgBsFVEQQV0QBF3hcUDYRUHZlFwSXR1RE9FUExYKoiFZkFQRkURBkEZStUBahhUJbknn/mKeBkNCmadq0ze9zXbk6mTmZuXMImTtnzpxjMQzDQERERKSYCPJ3ACIiIiK+pORGREREihUlNyIiIlKsKLkRERGRYkXJjYiIiBQrSm5ERESkWFFyIyIiIsVKCX8HUNDsdjsHDx6kdOnSWCwWf4cjIiIiHjAMg9TUVKpWrUpQUPZtMwGX3Bw8eJDY2Fh/hyEiIiJe2LdvH9WqVcu2TMAlN6VLlwbMyomMjPTpvjMzM/nhhx9o06YNwcHBPt13UaU6cU/14kp14kp14p7qxVUg1ElKSgqxsbGO83h2Ai65yboUFRkZmS/JTXh4OJGRkcX2w5VbqhP3VC+uVCeuVCfuqV5cBVKdeNKlRB2KRUREpFhRciMiIiLFipIbERERKVYCrs+Np2w2G5mZmbl6TWZmJiVKlODs2bPYbLZ8iqxouVSdBAcHY7Va/RiZiIgUV0puLmIYBocOHeLEiRNevTYqKop9+/ZpDJ3/ya5OypYtS1RUlOpKRER8SsnNRbISm8qVKxMeHp6rE6/dbufUqVNERETkOMBQoHBXJ4ZhkJaWRnJyMgDR0dH+DFFERIoZJTcXsNlsjsSmQoUKuX693W4nIyODsLAwJTf/c6k6KVmyJADJyclUrlxZl6hERMRndAa+QFYfm/DwcD9HEhiy6jm3fZtERESyo+TGDfUBKRiqZxERyQ+6LCUiIiI+YbPBihWQlATR0dCyJfij14GSGxEREcmzhATo3Rv27z+/rlo1mDoVOnUq2Fh0WSqf2GywdCl89pn511/D3nTv3p2OHTv65+AiIhIQEhLgvvucExuAAwfM9QkJBRuPWm7yQUIC9O1bOLLXqVOnYhhGwR5UREQChs1mtti4O9UYBlgs0KcPdOhQcJeo1HLjY19/HcwDD1gKTfZapkwZypYtW7AHFRGRgLFihWuLzYUMA/btM8sVFCU3OTAMOH3as0dKCgwcWPKS2SuY2W1Kimf7y02Dy9y5c6lXrx4lS5akQoUKxMfHc/r0aZfLUqmpqXTp0oVSpUoRHR3N66+/zi233EKfPn0cZeLi4nj11Vfp2rUrERERVK9enQULFnDkyBE6dOhAREQE9evX59dff3W85tixY3Tu3JmYmBjCw8OpV68en332WS5rW0REipqkJN+W8wUlNzlIS4OICM8e5coFkZQUBLi/xdkwzOy2TBnP9peW5lmMSUlJdO7cmccff5xt27axdOlSOnXq5PZyVL9+/fj5559ZsGABiYmJrFixgg0bNriUe/3112nevDkbN26kXbt2PProo3Tt2pVHHnmEDRs2cPnll9O1a1fHMc6ePUujRo349ttv2bJlCz179uTRRx9l7dq1Hte1iIgUPZ4OMl+Qg9Grz00xkJSUxLlz5+jUqRPVq1cHoF69ei7lUlNT+eijj5g1axatWrUC4MMPP6Rq1aouZe+8807+85//ADBs2DDeeecdrr/+eu6//34ABg4cSLNmzTh8+DBRUVHExMTwwgsvOF7/3HPP8f333/PFF18wdOhQn79nEREpHFq2hMqV4X8z6riwWMx+py1bFlxMSm5yEB4Op055VnbZMjvt2uXcGLZwIdx0k2fH9kSDBg1o1aoV9erVo23btrRp04b77ruPcuXKOZXbvXs3mZmZNGnSxLGuTJkyXHXVVS77rF+/vmO5SpUqgHPClLUuOTmZqKgobDYbY8aM4fPPP+fAgQNkZGSQnp7umGZBRESKp3PnIDTU/bassVqnTCnY8W50WSoHFguUKuXZo3VrqFrVjsXivrOMxQKxsdCmjWf783QAX6vVSmJiIt999x116tThzTff5KqrrmLPnj1ev+/g4OAL4rZccp3dbgdg4sSJTJ06lYEDB7JkyRI2bdpE27ZtycjI8DoGEREp/MaMMTsMR0bCxRcCqlWDuXM1zk2RZrXCuHFnANfEJL+zV4vFQvPmzRk5ciQbN24kJCSEefPmOZWpWbMmwcHBrFu3zrHu5MmT7NixI8/H//nnn+nQoQOPPPIIDRo0oGbNmj7Zr4iIFF6//24mNwDvvw9798KSJTBrlvl3z56CT2xAl6V8rn37TD7/3KBvX4vLODdTpuTPP/KaNWtYvHgxbdq0oXLlyqxZs4YjR45Qu3Ztfv/9d0e50qVL061bNwYMGED58uWpXLkyw4cPJygoKM/zPF1xxRXMnTuXVatWUa5cOSZPnszhw4epXbt2Xt+eiIgUQufOweOPm387doT77zd/yN9yi78jU3KTLzp1gnvuKbj5NSIjI1m+fDlTpkwhJSWF6tWrM2nSJO644w7mzJnjVHby5Mk89dRT3HXXXURGRvLiiy+yb98+wsLC8hTDkCFD2L17N23btiU8PJyePXvSsWNHTpw4kaf9iohI4TR5MqxfD2XLwn//63lXioKg5CafWK0Fl73Wrl2bRYsWud02c+ZMp+elS5fm008/dTw/ffo0I0eOpGfPno51f//9t8t+Lr6tPC4uzmld+fLlmT9/vsvr7HY7KSkpHrwLEREpKrZvh2HDzOXJkwv2Nm9PKLkJMBs3buTPP/+kSZMmnDx5klGjRgHQoUMHP0cmIiJFgd0OTz4J6enmDTLdu/s7IldKbgLQa6+9xvbt2wkJCaFRo0asWLGCihUr+jssEREpAv77X1i50ryr9913C9flqCxKbgJMw4YNWb9+vb/DEBGRIujvv2HQIHN53DiIi/NnNJemW8FFREQkR4YBPXuacx+2aAHPPOPviC5NyY2IiIjkaOZMSEyEsDD44AMIKsQZRCEOTURERAqDgwehXz9zeeRIuPJK/8aTEyU3IiIickmGYV6COnECGjU6n+QUZkpuRERE5JI+/xy++gpKlIAZM8y/hZ2Sm2LCMAx69uxJ+fLlsVgsbNq0yd8hiYhIEXf0KDz3nLn80ktQv75/4/GUkpv8YrPB0qXw2WfmX5stXw+3aNEiZs6cyTfffENSUhJ169bN0/5GjBjBtdde65vgRESkSOrdG44cgWuugZdf9nc0nisCjUtFUEIC9O2Ly8yZU6fm2/Sof/31F9HR0dx44435sn9vZWRk+DsEERHxwjffmLN7BwWZl6NCQvwdkefUcuNjwV9/jeWBB5wTG4ADB+C++8zEx8e6d+/Oc889x969e7FYLMTFxbFo0SJatGhB2bJlqVChAnfddRd//fWX0+v2799P586dKV++PKVKlaJx48asWbOGmTNnMnLkSH777TcsFgsWi8UxR9XevXvp0KEDERERREZG8sADD3D48GHHPrNafKZPn06NGjUIDw/3+fsVEZH8dfIkPPWUudyvHzRp4t94ckstNzkxDEhL86xsZiYlBw40X+NuPxaL2cYXH+/ZFOHh4R6Naz116lQuv/xy3nvvPdatW4fVamX58uX069eP+vXrc+rUKYYNG8Y999zDpk2bCAoK4tSpU9x8883ExMSwYMECoqKi2LBhA3a7nQcffJAtW7awaNEifvzxRwDKlCmD3W53JDbLli3j3LlzPPvsszz44IMsXbrUEc+uXbv48ssvSUhIwFIYx+UWEZFsDRhg/iavVcu89buoUXKTk7Q0iIjwqGiOzWCGYbbolCnj2bFPnTIn78hBmTJlKF26NFarlaioKADuvfdepzIzZsygUqVKbN26lbp16zJr1iyOHDnCunXrKF++PAC1atVylI+IiKBEiRKO/QEkJiayefNm9uzZQ2xsLAAff/wx11xzDevWreP6668HzEtRH3/8MZUqVdKs4CIiRczixfD+++byBx+Yv7OLGl2WKqZ27txJ586dqVmzJpGRkcT9bwKQvXv3ArBp0yYaNmzoSGw8sW3bNmJjYx2JDUCdOnUoW7Ys27Ztc6yrXr06lSpV8s0bERGRAnP6NPToYS4//TTcdJN/4/GWWm5yEh5utqB4wL5sGUHt2uVccOFCzz4xeUiX27dvT/Xq1Xn//fepWrUqdrudunXrOjr4lixZ0ut956SUB61NIiJS+Lz8MuzZA7Gx5sSYRZWSm5xYLB5dGgKgdWvsVatiSUrC4q7fjcVi3jXVpo1nfW68dOzYMbZv3877779Py5YtAVi5cqVTmfr16zN9+nSOHz/utvUmJCQE20W3r9euXZt9+/axb98+R+vN1q1bOXHiBHXq1MmndyMiIgVh1Sp44w1z+b33IDLSv/HkhS5L+ZLVypmsVPfijrRZz6dMydfEBqBcuXJUqFCB9957j127dvHTTz/R76Lxsjt37kxUVBQdO3bk559/Zvfu3Xz55ZesXr0agLi4OPbs2cOmTZs4evQo6enpxMfHU69ePbp06cKGDRtYu3YtXbt25eabb6Zx48b5+p5ERCT/nD0LTzxhdg3t2hVuv93fEeWNkhsfy2zfHuPzzyEmxnlDtWowd26+jXNzoaCgIGbPns369eupW7cuffv2ZeLEiU5lQkJC+OGHH6hcuTJ33nkn9erVY9y4cVj/l3jde++93H777dx6661UqlSJzz77DIvFwldffUW5cuW46aabiI+Pp2bNmsyZMyff35OIiOSfV16BP/+EKlXg9df9HU3e6bJUfujUCe65B1asgKQkiI6Gli3ztcWmT58+9OnTx/E8Pj6erVu3OpUxLrpUVr16debOnet2f6GhoW63XXbZZXz11VeXjGPEiBGMGDHC88BFRMSvNm6E8ePN5f/+F3Jxn0mhpeQmv1itcMst/o5CRETkkjIz4fHHzRmC7ruvQC4uFAhdlhIREQlQEybApk1ma81bb/k7Gt9RciMiIhKAtm6FUaPM5SlTzP42xYWSGxERkQBjs5l3R2VkwB13wCOP+Dsi31Jy48bFHW8lf6ieRUT848034ZdfoHRpePddj6YxLFKU3FwgODgYgDRPJ8qUPMmq56x6FxGR/Ld7tzkSMcDEieZoxMWN7pa6gNVqpWzZsiQnJwMQHh6eq1mt7XY7GRkZnD17lqAg5Y3gvk4MwyAtLY3k5GTKli3rGFtHRETyl2GYc0elpZk39GbNI1XcKLm5SNYs2FkJTm4YhsGZM2coWbJkrpKi4iy7OilbtqzTrOMiIpK/pk+Hn36CkiXN5eL6O1zJzUUsFgvR0dFUrlyZzMzMXL02MzOT5cuXc9NNN+lSy/9cqk6Cg4PVYiMiUoD274f+/c3lV1+Fyy/3bzz5ScnNJVit1lyffK1WK+fOnSMsLEzJzf+oTkRE/M8w4KmnIDUVmjaF3r39HVH+KqYNUiIiIpJl1iz49lsIDoYPPsj3+Zv9TsmNiIhIMZacfL6lZuhQuOYa/8ZTEJTciIiIFGPPPQfHjkH9+jBokL+jKRjqcyMiIlKM2GywYgUkJcFff8Hnn5uXoWbMMC9LBQIlNyIiIsXEvHkW+vc374y60N13Q6NG/onJH/x+Wertt98mLi6OsLAwmjZtytq1a7MtP2XKFK666ipKlixJbGwsffv25ezZswUUrYiISOG0enU0Dz1kdUlsAObPh4SEAg/Jb/ya3MyZM4d+/foxfPhwNmzYQIMGDWjbtu0lB9CbNWsWgwYNYvjw4Wzbto0PPviAOXPm8NJLLxVw5CIiIoWHzQbTp9cjuyn7+vQxywUCvyY3kydPpkePHjz22GPUqVOHadOmER4ezowZM9yWX7VqFc2bN+fhhx8mLi6ONm3a0Llz5xxbe0RERIqzlSstHDtWEnA/Or5hwL59Zl+cQOC3PjcZGRmsX7+ewYMHO9YFBQURHx/P6tWr3b7mxhtv5JNPPmHt2rU0adKE3bt3s3DhQh599NFLHic9PZ309HTH85SUFMAcOTe3IxDnJGt/vt5vUaY6cU/14kp14kp14p7qxdX+/XY8OaXv23eOzMxsmncKsdz8e/stuTl69Cg2m40qVao4ra9SpQp//vmn29c8/PDDHD16lBYtWmAYBufOneOpp57K9rLU2LFjGTlypMv6H374gfDw8Ly9iUtITEzMl/0WZaoT91QvrlQnrlQn7qleztu/vwLQIsdy//zzCwsXHsv/gPJBWlqax2UthpHdFbr8c/DgQWJiYli1ahXNmjVzrH/xxRdZtmwZa9ascXnN0qVLeeihh3j11Vdp2rQpu3btonfv3vTo0YOhQ4e6PY67lpvY2FiOHj1KZGSkT99TZmYmiYmJtG7dWlMN/I/qxD3ViyvViSvViXuqF1dnz2ZSo4aFY8fCcHdpymIxiImBnTvPFdnRiVNSUqhYsSInT57M8fztt5abihUrYrVaOXz4sNP6w4cPX3Km6KFDh/Loo4/y5JNPAlCvXj1Onz5Nz549efnllwlyM71paGgooaGhLuuDg4Pz7T9Ffu67qFKduKd6caU6caU6cU/14uyJJzYyYcL1LustFgALU6dCWFjRra/c/Fv7rUNxSEgIjRo1YvHixY51drudxYsXO7XkXCgtLc0lgcma3NJPDVAiIiKFQunSmbhrtalWDebOhU6dCj4mf/HrIH79+vWjW7duNG7cmCZNmjBlyhROnz7NY489BkDXrl2JiYlh7NixALRv357JkyfTsGFDx2WpoUOH0r59+1zP4C0iIlKcfPHFlQA8/TQ88IA5QnF0NLRsWfwnyryYX5ObBx98kCNHjjBs2DAOHTrEtddey6JFixydjPfu3evUUjNkyBAsFgtDhgzhwIEDVKpUifbt2zN69Gh/vQURERG/W7vWwu+/V6JECYOBAy1Ur+7viPzL79Mv9OrVi169erndtnTpUqfnJUqUYPjw4QwfPrwAIhMRESkaxo83GwI6dzaoXt39WDeBxO/TL4iIiIj3Nm+Gr78OwmIxGDAgQIYgzoGSGxERkSJs3Djzb7NmB7n6av/GUlgouRERESmidu2C2bPN5fvu2+nfYAoRJTciIiJF1IQJYLfD7bfbqVnzpL/DKTSU3IiIiBRBBw7AzJnm8qBBdr/GUtgouRERESmCXnsNMjPhppvgxhs1kO2FlNyIiIgUMUeOwHvvmcvZzB0dsJTciIiIFDFTp0JaGjRqBG3a+DuawkfJjYiISBFy8iS89Za5/NJLWRNjyoWU3IiIiBQh77xjJji1a0PHjv6OpnBSciMiIlJEpKXB5Mnm8uDBEKSzuFuqFhERkSLigw/MzsRxcfDQQ/6OpvBSciMiIlIEZGSYg/YBDBwIwcH+jacwU3IjIiJSBHzyCezfD9HR0L27v6Mp3JTciIiIFHI22/kJMvv3h7Aw/8ZT2Cm5ERERKeTmzoWdO6F8efjPf/wdTeGn5EZERKQQMwwYM8Zc7t0bIiL8G09RoORGRESkEPv2W/j9dzOp6dXL39EUDUpuRERECinDgNGjzeWnnzYvS0nOlNyIiIgUUsuWwS+/QGgo9Ovn72iKDiU3IiIihVRWq80TT0BUlH9jKUqU3IiIiBRCa9fCjz9CiRLw4ov+jqZoUXIjIiJSCGXdIdWlC1Sv7t9YiholNyIiIoXMli3w1VdgscCgQf6OpuhRciMiIlLIZI1GfO+9cPXV/o2lKFJyIyIiUoj89Rd89pm5/NJL/o2lqFJyIyIiUohMmAB2O9xxBzRs6O9oiiYlNyIiIoXEgQMwc6a5rFYb7ym5ERERKSQmTYKMDLjpJmjRwt/RFF1KbkRERAqBo0fh3XfNZbXa5I2SGxERkUJg6lRIS4NGjaBNG39HU7QpuREREfGzlBR4801z+aWXzPFtxHtKbkRERPzsv/+Fkyehdm3o2NHf0RR9Sm5ERET8KC0NJk82lwcPhiCdmfNMVSgiIuJHM2bAkSMQFwcPPeTvaIoHJTciIiJ+kpFhDtoH5szfwcH+jae4UHIjIiLiJ59+Cvv2QVQUPPaYv6MpPpTciIiI+IHNBmPHmsv9+0NYmH/jKU6U3IiIiPjBl1/Czp1Qrhw89ZS/oylelNyIiIgUMMOAMWPM5d69ISLCv/EUN0puRERECtjChfDbb2ZS89xz/o6m+FFyIyIiUoAMA0aPNpeffhrKl/dvPMWRkhsREZECtGwZrF4NoaHQt6+/oymelNyIiIgUoKy+Nk88AdHR/o2luFJyIyIiUkDWrYPERLBaYcAAf0dTfCm5ERERKSBZrTaPPGJOtyD5Q8mNiIhIAfjjD5g/HywWGDjQ39EUb0puRERECkDWaMSdOkHt2v6NpbgrkdsXbNu2jdmzZ7NixQr++ecf0tLSqFSpEg0bNqRt27bce++9hIaG5kesIiIiRdLu3fDZZ+bySy/5N5ZA4HHLzYYNG4iPj6dhw4asXLmSpk2b0qdPH1555RUeeeQRDMPg5ZdfpmrVqowfP5709PT8jFtERKRQs9lg6VIzqenTB+x2uP12uO46f0dW/HnccnPvvfcyYMAA5s6dS9myZS9ZbvXq1UydOpVJkybxktJTEREJQAkJ5rQK+/c7r2/Rwj/xBBqPk5sdO3YQHBycY7lmzZrRrFkzMjMz8xSYiIhIUZSQAPfdZ45EfLGhQ83+Np06FXxcgcTjy1KeJDZ5KS8iIlLU2Wxmi427xCZLnz5mOck/Pr1b6tdff2X58uW+3KWIiEiRsWKF66WoCxkG7NtnlpP8k+u7pbLz6KOPsmPHDmxKSUVEJAAlJfm2nHjHp8nN4sWL1ddGREQC0tGjMG+eZ2U1p1T+8mlyU7VqVV/uTkREpNBLTYXJk2HSJHM5OxYLVKsGLVsWTGyByqs+NwsXLuT77793Wf/999/z3Xff5TkoERGRwu7sWTOpqVkTRowwE5uGDWHIEDOJsVicy2c9nzLFnDhT8o9Xyc2gQYPc9qsxDINBgwblOSgREZHC6tw5mD4drrgC+vc3L0ddeSV8/jn8+iu88grMnQsxMc6vq1bNXK/bwPOfV5eldu7cSZ06dVzWX3311ezatSvPQYmIiBQ2druZnAwdCjt2mOuqVTNbbbp1gxIXnFE7dYIOHcy7opKSzD42LVuqxaageJXclClTht27dxN30Xztu3btolSpUr6IS0REpFAwDFi0CF5+GTZuNNdVrGjOEfX00xAW5v51VivcckuBhSkX8OqyVIcOHejTpw9//fWXY92uXbvo378/d999t8+CExER8aeff4abb4Y77zQTm9KlYeRIcyLMvn0vndiIf3mV3EyYMIFSpUpx9dVXU6NGDWrUqEHt2rWpUKECr732Wq729fbbbxMXF0dYWBhNmzZl7dq12ZY/ceIEzz77LNHR0YSGhnLllVeycOFCb96GiIiIW5s2Qbt25lxQK1ZAaKjZv2b3bhg2zExypPDy+rLUqlWrSExM5LfffqNkyZLUr1+fm266KVf7mTNnDv369WPatGk0bdqUKVOm0LZtW7Zv307lypVdymdkZNC6dWsqV67M3LlziYmJ4Z9//sl2Ik8RERFP7dxpJi+zZ5vPrVZ4/HFzXbVq/o1NPOf1ODcWi4U2bdpw0003ERoaiuXie948MHnyZHr06MFjjz0GwLRp0/j222+ZMWOG27uuZsyYwfHjx1m1apVj7qqL+/2IiIjk1oEDMGoUfPDB+XmfHnrIvAR15ZX+jU1yz6vkxm63M3r0aKZNm8bhw4fZsWMHNWvWZOjQocTFxfHEE0/kuI+MjAzWr1/P4MGDHeuCgoKIj49n9erVbl+zYMECmjVrxrPPPstXX31FpUqVePjhhxk4cCDWS3RBT09PJz093fE8JSUFgMzMTJ+Pppy1P43SfJ7qxD3ViyvViSvViXve1IvNBitXWhx3LrVoYWC1wrFjMHFiEP/9bxBnz5o/0u+4w87IkTauvTbreL5+B74XCJ+V3Lw3r5KbV199lY8++ogJEybQo0cPx/q6desyZcoUj5Kbo0ePYrPZqFKlitP6KlWq8Oeff7p9ze7du/npp5/o0qULCxcuZNeuXTzzzDNkZmYyfPhwt68ZO3YsI0eOdFn/ww8/EB4enmOc3khMTMyX/RZlqhP3VC+uVCeuVCfueVovq1dHM316PY4dK+lYV778GerUOcb69VU4c8b8cVynzlEeeWQbdeoc5+BBOHgwX8LOV8X5s5KWluZxWYthZDcxu3u1atXi3XffpVWrVpQuXZrffvuNmjVr8ueff9KsWTP+/fffHPdx8OBBYmJiWLVqFc2aNXOsf/HFF1m2bBlr1qxxec2VV17J2bNn2bNnj6OlZvLkyUycOJGkS8xC5q7lJjY2lqNHjxIZGZnbt56tzMxMEhMTad26teOyWaBTnbinenGlOnGlOnEvN/Uyb56Fhx6yYp7pLuw+YTieN2hg8MorNtq2NVxGFS4qAuGzkpKSQsWKFTl58mSO52+vWm4OHDhArVq1XNbb7XaPm40qVqyI1Wrl8OHDTusPHz5MVFSU29dER0cTHBzsdAmqdu3aHDp0iIyMDEJCQlxeExoaSmhoqMv64ODgfPsA5Oe+iyrViXuqF1eqE1eqE/dyqhebzbzDyf1PeDOLqVAB1q2zEBzs06kW/aY4f1Zy8768uhW8Tp06rFixwmX93LlzadiwoUf7CAkJoVGjRixevNixzm63s3jxYqeWnAs1b96cXbt2YbfbHet27NhBdHS028RGREQC14oVsH9/9mWOHTPHspHixatUddiwYXTr1o0DBw5gt9tJSEhg+/btfPzxx3zzzTce76dfv35069aNxo0b06RJE6ZMmcLp06cdd0917dqVmJgYxo4dC8DTTz/NW2+9Re/evXnuuefYuXMnY8aM4fnnn/fmbYiISDF2id4KXpeTosOr5KZDhw58/fXXjBo1ilKlSjFs2DCuu+46vv76a1q3bu3xfh588EGOHDnCsGHDOHToENdeey2LFi1ydDLeu3cvQUHnG5diY2P5/vvv6du3L/Xr1ycmJobevXszcOBAb96GiIgUY9HRvi0nRYfXFxlbtmzpk17ZvXr1olevXm63LV261GVds2bN+OWXX/J8XBERKd5atoSoKDh0yP12i8UcmK9ly4KNS/KfV31u9u3bx/4LLmSuXbuWPn368N577/ksMBERkbywWuGyy9xvy7orasoUzdRdHHmV3Dz88MMsWbIEgEOHDhEfH8/atWt5+eWXGTVqlE8DFBER8cbKlbB2rZnIXDSkGtWqwdy50KmTf2KT/OVVcrNlyxaaNGkCwOeff069evVYtWoVn376KTNnzvRlfCIiIrlmGDBggLnco4c5vcKSJTBrlvl3zx4lNsWZV31uMjMzHWPH/Pjjj9x9990AXH311ZccTE9ERKSgfPkl/PILhIfDiBHmpadbbvF3VFJQvGq5ueaaa5g2bRorVqwgMTGR22+/HTBHHa5QoYJPAxQREcmNjAzImrbwhRd0N1Qg8iq5GT9+PO+++y633HILnTt3pkGDBoA5sWXW5SoRERF/eO892LXL7Gfzwgv+jkb8wavLUrfccgtHjx4lJSWFcuXKOdb37Nkz3yajFBERyUlKCmTNlTxiBJQu7ddwxE+8HufGarU6JTYAcXFxeY1HRETEa+PHw9GjcNVV8MQT/o5G/MXjy1K33367R4PnpaamMn78eN5+++08BSYiIpIb+/fD5Mnm8rhxUEznjxQPeNxyc//993PvvfdSpkwZ2rdvT+PGjalatSphYWH8+++/bN26lZUrV7Jw4ULatWvHxIkT8zNuERERJ8OHw9mz0KIFdOjg72jEnzxObp544gkeeeQRvvjiC+bMmcN7773HyZMnAbBYLNSpU4e2bduybt06ateunW8Bi4iIXGzzZsgaZm3ixPMjEEtgylWfm9DQUB555BEeeeQRAE6ePMmZM2eoUKECwWr/ExERPxk4EOx2uO8+uOEGf0cj/uZ1h2KAMmXKUKZMGV/FIiIikmuLF8N330GJEjBmjL+jkcLAq3FuRERECgO7HV580Vx++mm44gr/xiOFg5IbEREpsmbPhg0bzPFshg71dzRSWCi5ERGRIunsWXjpJXN50CCoVMm/8UjhoeRGRESKpHfeCeKffyAmBvr08Xc0Uph4ndycOHGC6dOnM3jwYI4fPw7Ahg0bOHDggM+CExERcSc1NZixY81T2CuvmLN/i2Tx6m6p33//nfj4eMqUKcPff/9Njx49KF++PAkJCezdu5ePP/7Y13GKiIg4fPnlFZw4YaFuXeja1d/RSGHjVctNv3796N69Ozt37iQsLMyx/s4772T58uU+C05ERORif/8N33xTE4AJE8Bq9W88Uvh4ldysW7eO//znPy7rY2JiOHToUJ6DEhERuZThw62cO2fl1lvt3H67v6ORwsir5CY0NJSUlBSX9Tt27KCSuquLiEg+2bABPvvMPHWNG2fTNAvillfJzd13382oUaPIzMwEzLml9u7dy8CBA7n33nt9GqCIiAiAYcCAAebyzTfvo2FD/8YjhZdXyc2kSZM4deoUlStX5syZM9x8883UqlWL0qVLM3r0aF/HKCIiwqJF8NNPEBJi8PDD2/wdjhRiXt0tVaZMGRITE1m5ciW///47p06d4rrrriM+Pt7X8YmIiGCznZ9m4dln7VSpcsa/AUmhlqeJM1u0aEGLFi18FYuIiIhbH38MW7ZA2bIwcKCdX37xd0RSmHmd3Kxbt44lS5aQnJyM3W532jZ58uQ8ByYiIgKQlnZ+3qghQ6B8ef/GI4WfV8nNmDFjGDJkCFdddRVVqlTBckF3dYu6rouIiA9NmQIHDkD16vDss/6ORooCr5KbqVOnMmPGDLp37+7jcERERM47cgTGjTOXR4+GsDD43426Ipfk1d1SQUFBNG/e3NexiIiIOHnlFUhNheuug86d/R2NFBVeJTd9+/bl7bff9nUsIiIiDrt2wTvvmMsTJ0KQ11M9S6Dx6rLUCy+8QLt27bj88supU6cOwcHBTtsTEhJ8EpyIiASul16Cc+fgjjvgttv8HY0UJV4lN88//zxLlizh1ltvpUKFCupELCIiPvXLL/DFF2CxwPjx/o5GihqvkpuPPvqIL7/8knbt2vk6HhERCXAXTrPQvTvUq+fXcKQI8uoKZvny5bn88st9HYuIiAgLFsDKlVCyJIwa5e9opCjyKrkZMWIEw4cPJy0tzdfxiIhIADt3DgYONJf79oVq1fwbjxRNXl2WeuONN/jrr7+oUqUKcXFxLh2KN2zY4JPgREQksEyfDtu3Q8WK5+eSEsktr5Kbjh07+jgMEREJdKmpMHy4uTxsGJQp4994pOjyKrkZnvXpExER8ZFJkyA5GWrVgv/8x9/RSFGmIZFERMTvkpLgtdfM5bFjISTEv/FI0eZxy0358uXZsWMHFStWpFy5ctmObXP8+HGfBCciIoFhxAg4fRqaNoV77/V3NFLUeZzcvP7665QuXdqxrIH7RETEF7ZtMzsSg9l6o9OL5JXHyU23bt0cy5oNXEREfGXQILDboWNHaNHC39FIceBVnxur1UpycrLL+mPHjmG1WvMclIiIBIbly81B+6xWs6+NiC94ldwYhuF2fXp6OiHqBSYiIh64cJqFHj3g6qv9G48UH7m6FfyNN94AwGKxMH36dCIiIhzbbDYby5cv52p9OkVExANffAFr10KpUufHtxHxhVwlN6+//jpgttxMmzbN6RJUSEgIcXFxTJs2zbcRiohIsZORAYMHm8sDBkBUlH/jkeIlV8nNnj17ALj11ltJSEigXLly+RKUiIgUPzYbrFhhjmmzahXs3m0mNf37+zsyKW68GqF4yZIlvo5DRESKsYQE6N0b9u93Xt+hA1zQw0HEJzRCsYiI5KuEBLjvPtfEBuC998ztIr6k5EZERPKNzWa22FziJlsA+vQxy4n4ipIbERHJNytWuG+xyWIYsG+fWU7EV5TciIhIvklK8m05EU941aEY4MSJE6xdu5bk5GTsdrvTtq5du+Y5MBERKfo8vcU7Ojp/45DA4lVy8/XXX9OlSxdOnTpFZGSk0ySaFotFyY2IiLB3L4wenX0ZiwWqVYOWLQsmJgkMXl2W6t+/P48//jinTp3ixIkT/Pvvv47H8ePHfR2jiIgUIYZhzvJdty4sXgzBweb6i2f7zno+ZYo5t5SIr3iV3Bw4cIDnn3+e8PBwX8cjIiJF2P79cMcd5lxRqanQrBls3gxffgkxMc5lq1WDuXOhUyf/xCrFl1eXpdq2bcuvv/5KzZo1fR2PiIgUQYYBH31k3tZ98iSEhsKrr0LfvmarzFVXmQP2ZY1QHB1tXopSi43kB6+Sm3bt2jFgwAC2bt1KvXr1CM5qc/yfu+++2yfBiYhI4XfwIPTsCd9+az5v0gRmzoTatZ3LWa1wyy0FHZ0EIq+Smx49egAwatQol20WiwWbRmMSESn2DAM++QSefx5OnICQEBg1ypwrqoTX9+KK5J1XH7+Lb/0WEZHAcugQ/Oc/sGCB+bxxY7O15ppr/BqWCKBB/EREJBcMAz77zExiFiww74QaPRpWr1ZiI4WH18nNsmXLaN++PbVq1aJWrVrcfffdrPBy/Oy3336buLg4wsLCaNq0KWvXrvXodbNnz8ZisdCxY0evjisiIp5LTjYnwHz4YTh+HBo2hPXr4aWXdBlKChevkptPPvmE+Ph4wsPDef7553n++ecpWbIkrVq1YtasWbna15w5c+jXrx/Dhw9nw4YNNGjQgLZt25KcnJzt6/7++29eeOEFWmrkJxGRfPf552bLTEKCmciMHAlr1kC9ev6OTMSVV8nN6NGjmTBhAnPmzHEkN3PmzGHcuHG88sorudrX5MmT6dGjB4899hh16tRh2rRphIeHM2PGjEu+xmaz0aVLF0aOHKnb0UVE8tGRI/DAA/Dgg3D0KDRoAOvWwbBh5wfnEylsvGpI3L17N+3bt3dZf/fdd/PSSy95vJ+MjAzWr1/P4MGDHeuCgoKIj49n9erVl3zdqFGjqFy5Mk888USOl8LS09NJT093PE9JSQEgMzOTzMxMj2P1RNb+fL3fokx14p7qxZXqxJW/6yQhwcJzz1k5csSC1WowaJCdwYPthISAP/+Z/F0vhVEg1Elu3ptXyU1sbCyLFy+mVq1aTut//PFHYmNjPd7P0aNHsdlsVKlSxWl9lSpV+PPPP92+ZuXKlXzwwQds2rTJo2OMHTuWkSNHuqz/4Ycf8m2E5cTExHzZb1GmOnFP9eJKdeKqoOskJSWY99+vz4oV1QC47LIUevfewOWXn+THHws0lGzps+KqONdJWlqax2W9Sm769+/P888/z6ZNm7jxxhsB+Pnnn5k5cyZTp071ZpceSU1N5dFHH+X999+nYsWKHr1m8ODB9OvXz/E8JSWF2NhY2rRpQ2RkpE/jy8zMJDExkdatW7sMbBioVCfuqV5cqU5c+aNOFiywMGCAlcOHzdaaAQPsvPxySUJDmxfI8T2hz4qrQKiTrCsvnvAquXn66aeJiopi0qRJfP755wDUrl2bOXPm0KFDB4/3U7FiRaxWK4cPH3Zaf/jwYaKiolzK//XXX/z9999Ol8SyxtwpUaIE27dv5/LLL3d6TWhoKKGhoS77Cg4OzrcPQH7uu6hSnbinenGlOnHl6zqx2VynQTh5Enr3NgflA6hTB2bOtHD99VagcM6RoM+Kq+JcJ7l5X17fvHfPPfdwzz33ePtyAEJCQmjUqBGLFy923M5tt9tZvHgxvXr1cil/9dVXs3nzZqd1Q4YMITU1lalTp+bqkpiISCBKSDCTmP37z6+rUMFMeE6cgKAgGDAARoyAsDB/RSmSN34fmaBfv35069aNxo0b06RJE6ZMmcLp06d57LHHAOjatSsxMTGMHTuWsLAw6tat6/T6smXLArisFxERZwkJ5jg1huG8/tgx82/VqmaZpk0LPjYRX/I4uSlfvjw7duygYsWKlCtXDovFcsmyx48f9ziABx98kCNHjjBs2DAOHTrEtddey6JFixydjPfu3UtQkAZSFhHJC5vNbLG5OLG5UFCQOY2CSFHncXLz+uuvU7p0acdydslNbvXq1cvtZSiApUuXZvvamTNn+iwOEZHiasUK50tR7uzfb5bTzN1S1Hmc3HTr1s2x3L179/yIRURE8klSkm/LiRRmXl3vsVqtbqdHOHbsGFZr4exVLyISyKKjfVtOpDDzKrkxLnHRNj09nZCQkDwFJCIivlejBmT329NigdhY87ZwkaIuV3dLvfHGGwBYLBamT59ORESEY5vNZmP58uVcffXVvo1QRETy5NgxuPNOs1MxmInMhb9Rs7pQTpmSfQIkUlTkKrl5/fXXAbPlZtq0aU6XoEJCQoiLi2PatGm+jVBERLx26pSZ2GzdCjExMHQovPqqc+fiatXMxKZTJ7+FKeJTuUpu9uzZA8Ctt95KQkIC5cqVy5egREQk79LTzYRl7VooXx5++MEcefjJJ11HKFaLjRQnXg3it2TJEl/HISIiPmSzwaOPQmIilCoF331nJjZgJjK63VuKM6+Sm8cffzzb7TNmzPAqGBERyTvDgGefhS++gOBgmD8fmjTxd1QiBcer5Obff/91ep6ZmcmWLVs4ceIEt912m08CExER7wwZAu++a3YUnjUL4uP9HZFIwfIquZk3b57LOrvdztNPP+0yK7eIiBScyZNhzBhzedo0cy4pkUDjs0mbgoKC6Nevn+OOKhERKVgffQT9+5vLY8ZAz57+jUfEX3w6I+Vff/3FuXPnfLlLERHxwIIF8MQT5nL//jBokH/jEfEnry5L9evXz+m5YRgkJSXx7bffOs1BJSIi+W/ZMnjgAfMOqe7dYeLE8wPziQQir5KbjRs3Oj0PCgqiUqVKTJo0Kcc7qURExHc2boS77zbHtLn7bnj/fSU2IhrnRkSkiNqxA9q2hZQUuPlmmD0bSnj1rS5SvPi0z42IiBSMAwegTRs4cgQaNoSvvoKSJf0dlUjh4HGO37BhQywetnVu2LDB64BERCR7x4+bic0//8AVV8CiRVCmjL+jEik8PE5uOnbsmI9hiIiIJy6eCDMxESpX9ndUIoWLx8nN8OHD8zMOERHJQXo63HsvrFlzfiLM6tX9HZVI4ZOnrmfr169n27ZtAFxzzTU0bNjQJ0GJiIgzmw26djUTmlKlYOHC8xNhiogzr5Kb5ORkHnroIZYuXUrZsmUBOHHiBLfeeiuzZ8+mUqVKvoxRRCSgZU2E+fnn5kSYCQnQtKm/oxIpvLy6W+q5554jNTWVP/74g+PHj3P8+HG2bNlCSkoKzz//vK9jFBEJaCNGBDkmwvzkE7MzsYhcmlctN4sWLeLHH3+kdu3ajnV16tTh7bffpo3+14mI+MyCBTWZMcMKwDvvmCMRi0j2vGq5sdvtBAcHu6wPDg7GbrfnOSgREYH/+z8LM2bUA2D0aPjPf/wckEgR4VVyc9ttt9G7d28OHjzoWHfgwAH69u1Lq1atfBaciEig+vpr6NnTbLHp3dvG4MF+DkikCPEquXnrrbdISUkhLi6Oyy+/nMsvv5waNWqQkpLCm2++6esYRUQCyvLlWRNhWrj11r2MH2/XfFEiueBVn5vY2Fg2bNjAjz/+yJ9//glA7dq1iY+P92lwIiKBZtMmaN8ezp6Fdu3sPPbYJoKCov0dlkiR4vU4NxaLhdatW9O6dWtfxiMiEjBsNlixApKSIDoaoqLOT4R5000wa5aNJUsMf4cpUuTk6rLU6tWr+eabb5zWffzxx9SoUYPKlSvTs2dP0tPTfRqgiEhxlJAAcXFw663w8MPm37p1ITkZrr0WFizQRJgi3spVcjNq1Cj++OMPx/PNmzfzxBNPEB8fz6BBg/j6668ZO3asz4MUESlOEhLgvvtg/37n9Tab+fe55zQRpkhe5Cq52bRpk9PdULNnz6Zp06a8//779OvXjzfeeIPPP//c50GKiBQXNhv07m2OOuyOxQIjRpxPdEQk93KV3Pz7779UqVLF8XzZsmXccccdjufXX389+/bt8110IiLFzIoVri02FzIM2LfPLCci3slVclOlShX27NkDQEZGBhs2bOCGG25wbE9NTXU7uJ+IiJiSknxbTkRc5Sq5ufPOOxk0aBArVqxg8ODBhIeH07JlS8f233//ncsvv9znQYqIFAenT5uzeXsiWnd/i3gtV7eCv/LKK3Tq1Imbb76ZiIgIPvroI0JCQhzbZ8yYobmlREQuYhjmjN4vvJD9JSkw+9xUqwYtW4JmsxHxTq6Sm4oVK7J8+XJOnjxJREQEVqvVafsXX3xBRESETwMUESnKfv8dnn8eli0zn8fFmaMPT5xoPr+wY3HWKMRTpoDVquRGxFteTb9QpkwZl8QGoHz58k4tOSIiger4cejVCxo2NBObsDAYORK2boXx42HuXIiJcX5NtWrm+k6d/BOzSHHh9QjFIiLiymaD6dPh5Zfh2DFz3X33wWuvQfXq58t16gQdOjiPUNyypdliIyJ5o+RGRMRHVq40B+DbtMl8fs018MYbcNtt7stbrXDLLQUVnUjg8OqylIiInHfgADzyiNnysmkTlC1rJjWbNl06sRGR/KOWGxERL6Wnw+uvw6uvmrd5Wyzw5JMwejRUquTv6EQCl5IbEREvfPst9OkDu3aZz5s1gzffhEaN/BqWiKDLUiIiubJjB7RrB3fdZSY2UVHw8cdmfxslNiKFg5IbEREPpKbCoEFQt645ynBwMAwYYCY7jz4KQfo2FSk0dFlKRCQbhgGffgovvnh+vqfbbzcH2rvqKr+GJiKXoORGRAKezeZ+vJkNG8xbu1etMsvVrGkmNXfddX40YREpfJTciEhAS0iA3r2d53yqWtUco+bHH82Wm/BwGDIE+vY1RxoWkcJNyY2IBKyEBHP04AvndwI4eNB8AHTuDBMmmFMjiEjRoORGRAKSzWa22Fyc2FyocmX4v//TlAgiRY3694tIQFqxwvlSlDvJyWY5ESlalNyISEDKuvPJV+VEpPBQciMiAenIEc/KRUfnbxwi4ntKbkQkoBgGvPaaeedTdiwWiI01bwsXkaJFyY2IBIwTJ6BTJ3NkYbvdTFwsFtcxa7KeT5mizsQiRZGSGxEJCJs2QePGMH8+hITAf/8Ly5bB3LkQE+Nctlo1c32nTv6IVETySreCi0ixN2MGPPssnD0L1avDF1/A9deb2zp1gg4d3I9QLCJFk5IbESm20tKgVy/48EPz+Z13muPWlC/vXM5qhVtuKfDwRCSf6LKUiBRLu3ZBs2ZmYhMUBKNHw9dfuyY2IlL8qOVGRIqdefOge3dISTFHGf7sM7jtNn9HJSIFRS03IlJsZGbCCy+Y/WhSUqBFC3NmbyU2IoFFLTciUiwcPAgPPggrV5rP+/eHsWMhONi/cYlIwVNyIyJF3k8/mbN3JydDZKTZz0a3cYsErkJxWertt98mLi6OsLAwmjZtytq1ay9Z9v3336dly5aUK1eOcuXKER8fn215ESm+7HYYMwZatzYTmwYNYP16JTYigc7vyc2cOXPo168fw4cPZ8OGDTRo0IC2bduSnJzstvzSpUvp3LkzS5YsYfXq1cTGxtKmTRsOHDhQwJGLiD8dPw7t28PLL5tJzmOPwerVUKuWvyMTEX/z+2WpyZMn06NHDx577DEApk2bxrfffsuMGTMYNGiQS/lPP/3U6fn06dP58ssvWbx4MV27dnUpn56eTnp6uuN5SkoKAJmZmWRmZvryrTj25+v9FmWqE/dUL65yUyfr11t46CEr//xjISzM4I03bHTvbvzv9fkaZoHS58Q91YurQKiT3Lw3i2EYRj7Gkq2MjAzCw8OZO3cuHTt2dKzv1q0bJ06c4KuvvspxH6mpqVSuXJkvvviCu+66y2X7iBEjGDlypMv6WbNmER4enqf4RaRgGQYsWhTHBx/U5dw5K1FRp3jxxXXUrJni79BEJJ+lpaXx8MMPc/LkSSIjI7Mt69eWm6NHj2Kz2ahSpYrT+ipVqvDnn396tI+BAwdStWpV4uPj3W4fPHgw/fr1czxPSUlxXMrKqXJyKzMzk8TERFq3bk2wbtEAVCeXonpxlVOdnD4Nzzxj5bPPzKvpd99t54MPQilTpkVBh1pg9DlxT/XiKhDqJOvKiyf8flkqL8aNG8fs2bNZunQpYWFhbsuEhoYSGhrqsj44ODjfPgD5ue+iSnXinurFlbs6+fNPuPde2LrVnCph/Hjo1y8Ii8Xv3QYLhD4n7qleXBXnOsnN+/JrclOxYkWsViuHDx92Wn/48GGioqKyfe1rr73GuHHj+PHHH6lfv35+hikiBcRmg2XLLBw5cn4Cyy++gCefNFtuoqNhzhxzvYjIpfg1uQkJCaFRo0YsXrzY0efGbrezePFievXqdcnXTZgwgdGjR/P999/TuHHjAopWRPLTvHkWnnmmDceOnf9aioiAU6fM5VtvNadRuOgqtoiIC79flurXrx/dunWjcePGNGnShClTpnD69GnH3VNdu3YlJiaGsWPHAjB+/HiGDRvGrFmziIuL49ChQwBEREQQERHht/chIt5LSICHHrJiGFan9VmJTadOZotNCb9/Y4lIUeD3r4oHH3yQI0eOMGzYMA4dOsS1117LokWLHJ2M9+7dS1DQ+evq77zzDhkZGdx3331O+xk+fDgjRowoyNBFxAdsNujd27wTCixuy6xbBxb3m0REXPg9uQHo1avXJS9DLV261On533//nf8BiUiBWbEC9u+HSyU2APv2meVuuaWgohKRoiwwbjUQkUJr2TLPyiUl5W8cIlJ8FIqWGxEJPOvWwfDh8N13npWPjs7feESk+FDLjYgUqI0b4e67oUkTM7EJCoJSpcBicT9YusUCsbG6/VtEPKfkRkQKxObN5kB8110HX39tJjXdusGOHfDxx1mlnBOcrE7EU6aYg/eJiHhCyY2I5Ktt2+Chh6BBA/OWb4sFHn7YXD9zJlx+uXmr9+zZNipUOOv02mrVYO5cc7uIiKfU50ZE8sXOnTBqFMyaBXa7ue7++2HECKhTx7X8PfcYlCjxA5GR7ThypIRjhGK12IhIbim5ERGf2r0bXn3VvNRks5nrOnaEkSMhp5lSrFa4+WaDYjo1jogUECU3IuIT//wDo0fDhx/CuXPmurvuMltqGjXya2giEmCU3IhInhw4AGPGwPvvQ2amua5tW7OlpmlT/8YmIoFJyY2IeOXQIRg3DqZNg/R0c91tt5n9bJo3929sIhLYlNyIiAubzZzuICkJl469R47AhAnw9ttw5oy5rmVLM6nR9AgiUhgouRERJwkJ5kSW5nxPpmrVzE7C27fDG2/A6dPm+htugFdegVatNLGliBQeSm5ExCEhAe67L2uG7vP274fu3c8/b9zYbKm5/XYlNSJS+Ci5ERHAvBTVu7drYnOh4GCYM8e8tVtJjYgUVhqhWEQAWL7c+VKUO5mZUK6cEhsRKdzUciMSwDIyzKTmm2/gs888e01SUv7GJCKSV0puRALM4cOwcKGZ0CQmQmpq7l4fHZ0/cYmI+IqSG5FizjBg0yYzmfnmG1i71nl7VBS0awd33GH2uTl40H2/G4vFvGuqZcsCCVtExGtKbkSKodOnYfFiM5n59lszYblQ48ZmQnPXXXDddRD0v953Fot5t5TF4pzgZPWxmTJFE1mKSOGn5EakEMtuML2L/fPP+WTmp5/OjxoMUKoUtG5tJjN33nnpS0udOsHcue7HuZkyxdwuIlLYKbkRKaQuNZje1KlmknHuHPzyi5nMfPMNbNni/Pq4OGjf3kxobr4ZQkM9O26nTtChg+dJlYhIYaPkRqQQutRgegcOwL33msnGH3/A8ePnt1mt5pxOWZebatf2/pZtq1VTKYhI0aXkRqSQyW4wvax1K1aYf8uVMzsC33WXORN3+fIFF6eISGGl5EakkFmxIufB9MC8PPXMM1BC/4tFRJxohGKRQsQw4PvvPStbqZISGxERd5TciBQCdjt89RU0aQLjxnn2Gg2mJyLinn73ifiRzWbeej1mzPm7ncLCzBaZ06c1mJ6IiDfUciPiBxkZ8OOPl1GvXgkefthMbEqXhsGDzfFqPvrILHfx3U4aTE9EJGdKbkQK0Jkz8PbbUKdOCd56qyG7dlkoXx5GjTKTmjFjoHLl84PpxcQ4v75aNXO9BtMTEbk0XZYSKQCnTsG0aTBpEhw6BGChXLmzDBoUzDPPWImIcH2NBtMTEfGOkhuRfHTiBLz5pnkZKWvAvcsugxdesBEVlUjHjrcTHHzpbEWD6YmI5J6SG5F8kJwMr79uXoJKTTXXXXGF2aemSxewWOwsXGj3b5AiIsWUkhsRH9q/H157Dd57z+xfA1CvHrz0Etx///lLSpmZ/otRRKS4U3Ij4oGcZufevRvGj4cPPzyfuFx/PQwZYk6NEKSu+yIiBUbJjUgOspudu3ZtGDsWZs0yEyCAm24yk5r4eO8nrhQREe8puRHJRk6zc1/o9tvh5ZehRYuCi09ERFwpuRG5BE9m5wbo2NFMaho3LrDQREQkG+oJIHIJns7O3bu3EhsRkcJELTciF0lPh+XLYfJkz8onJeVvPCIikjtKbkSAgwdh4UL49lv48UdzRGFPaXZuEZHCRcmNBCSbDdatM5OZb7+FjRudt0dHmx2EFywwRxbW7NwiIkWHkhsJGCdOwPffm8nMd9/B0aPnt1ks0KQJtGtnPq691hybJutuKYvFOcHR7NwiIoWXkhspUnIaTO9ChgFbt55vnfn55/Nj0QCUKQNt25rJzO23m7NxXyxrdm5349xMmaLZuUVECiMlN1JkZDeYXlaSceYMLFlyPqH55x/nfdSpc7515sYbITg45+Nqdm4RkaJFyY0UCdkNpnffffDkk2an4J9+Oj+nE0BoKNx2m5nM3Hkn1Kjh3fE1O7eISNGh5EYKPU8G03v//fPrYmPPt87cdhuEhxdMnCIiUjgouZFC7/vvPRtMr0cPeO45qFtXczqJiAQyJTdS6Ozda3b+zXr89ptnr7v1VqhXL39jExGRwk/JjfjVuXPw++/OyYwnrTTuaDA9EREBJTfihdzcjn2xkyfhl1/OJzJr1sDp085lrFZo2BCaNzcfN9xg3tl04IAG0xMRkZwpuZFc8eR27CyGAX//DcuWVeO774JYvRo2b3ZNUMqUgWbNziczTZpAqVLOZaZO1WB6IiLiGSU34rGcbseePRtq1nS+xHTwYDDQyKl8jRrnE5nmzeGaa8zRgLOjwfRERMRTSm7EI57cjv3QQ67bS5QwqFHjBHfcEUnLllaaN/e+b4wG0xMR8UBe+g4U5WNfQMmN5Mhmg08/zbmjr2GYl5NuvtlskbnxRmjY8BxLly7nzjvvJDg47x9wDaYnUoQF2knXD8e0zJsH/ft71nfA13LTbyGfKbkp4nz9fyc52bx7afPm83//+APOnvXs9e++C126nH+emel9LCLFnr9O9v44rj9PfP44th+OGb16NdYJEy7dd2Du3Px9v9n1W8jPY7uh5KYIy8v/nTNnzEklL05kkpPdlw8JgYwMczkIGy1ZQTRJJBHNClpix/xijInxwRuTwGWzYVm2DI4c0cm+OB3Xnyc+fxzbH8e02ag3ffql+w5YLOYop9dfbz7PzDTH4vDF34wMGD8++2P36WP2KyigljolN0VU1v8di2Hj5gsSjZX7W3LffVbH/x27HfbscU5gNm+GnTvNbRezWODyy6F+fXNAvKy/1aub65vsT2AKvYnl/BfjPqrRh6msi+1UfG/HDqRf2H46rmXePNo88wwljh07v1In+6J/3Jw67OXlxJd1kr74kXXSPXsWnnkm+86CTz8NZcuev6vBMM4/Lnx+qeWLn9ts0LNn9sd88klztNJz58zEwAePEqmpBB8/nn1dHTwIl12Wqyr2CcOAffvM75QC6leg5MZHbBk2fn9zOak/beH3XRFc+9wtWEPy52SQ9V3R0UhgqptEo7cxla5dOzF+vHlJ6eJxZLJUrOicwNSvb86affFt2Fm+6JxAk4n3Ac7/aWM4wBfcx9qH5mK15nOzo35hF8/jJiRgfeghrDrZmyf75583x0fIzCTs6FFzevugIPMXic3m+nC3/uJ1GRmenXSTks7HmbWPvPzdty/7DntZJ77rrjO/fLISk2weJc6do31GBkE2W97/PZKToVWrvO8nN/79F/r29ekuPZ5xxmIxm+KDg6FEiUv/zW7bxX/37YNly3I+dtZnqwBYDMPdJ734SklJoUyZMpw8eZLIyEif7POXFxO4bHJvqtrO/wc+aK3G3n5TuWFC3r8Y09LMVvrkZPOxfDnsnJDAXMxE48K7qO3/+4jfx1zmYR47NNRMWi5OZKpUycUcTDYbxMVh7N/v9j+RgQVLbDWzmeiCZCMzM5OFCxf+r0NxsFfv36Ew/cLOqjgvT3451ks+HTdH/jhuDp8txyiNF322PHapX/dnz5o93w8duvRrK1SA1183y6ene//IyDi/nJIC2f3ClrwJCjp/4jUM8ws0J1WrmgNugfl5y/q8Zy1f/Dy75ePHYffunI95ww1wxRVmouGDx7kNGyjx9NM5H3fJEt+3nixdas5/k8/Hzs35W8lNHv3y4vnWDHdJxtoBc10SnIyM88nKhUnLhcvJyfBvcianktMg7TSlOE04aZTiNBGkMosuVODYJRINSCGSddf9h4b17ZQrfY4g+znzF1F2j6zmXHePkyfNX485adsWrr7a/KIoU4ZzERFs+OsvrrvlFkpUqOBYT5kyEBbmeXblx5PuJX91envSNQwy09NZ9M033N66NcFBQeaxzp07/wu7adNL/8qxWKByZfjmm/O/5g0j738zM81f9hdeFrpY+fIwbtz5113cSpDTX3fr9u2DBQtyrreGDSEiwvXSQ04Pd9dfiwgjKAjDYsFSogQWq9X8nGU9goKcn+e07d9/YdeunA96/fXm5z5rH3n9+/ffMGNGzscdPhyuvfZ8YpLNI9Mw+GnFCm5r25bg8HDn7RcOmlVAJ10n/jgmkHn2LOeqVSPs+HEs2Q3l7u2PhOxkfVfmNIx8Ho+t5CYbvkxubBk2DofHEWXb75TYZLEDJynH/Gte5lzqGWypaXD6NNaMNKdk5cLlC/+GUMxvNQoOdk52LvUoXRpeein7k27lyuYogjbb+Q5uGRnnly/+68m2Awdg5cqc30e1ambzWNbJOytJuXg563kRPtkWG1knX09u56tb1+x0Fhrq/AgJcV2X3SOr/ObNZhKZk59+IrNFC9+1fPrppJsfJz6PW4QL6KTr92Ni1snGoUO5fsIE80evu6HcC6Ljdj4eOzfnb/W5yYPN/13BtbZLX0sOAsrxL4/98UKejmMEBUGpUlhKlYLwcIzMTCz79uX4OvvtdxBUr67rddTsHpcqt2WLZ9eIe/Qwm/JPnoSTJ7GfOMG/e/ZQ3mrFkpJirk9JOd9ScPSo+cir5GS47ba878cb3s70eSkXzzFxKWXLmn0ULBbzRH3xX3frsit77Jhnv+yvuw5iY93/Ss/NL/qs5b174aOPcj7ukCHQoIHz9f6cHu7KlShhHtfTk/2bb/r2ZH/99TBqVM4nv5tu8m0i3LKlud+CnqTNavXf/Cn+OLYf329Ss2bYZs+mhLtxbvJ7KPdCNox8oUhu3n77bSZOnMihQ4do0KABb775Jk2aNLlk+S+++IKhQ4fy999/c8UVVzB+/HjuvPPOAozYlPaXZ52jdlS+kRL16hBSrhRh5cMJr1iKsIqlCCoVbp6cwi/6e9E6S0iI0+Ubi4dfykEDX/Tdl/Ktt8KkSTl/Mb7zjtN/WltmJisv/oVlt8OpU44EKNvHiRPmPeubNuUcY3S0mVhlXYcODnb/N7ttF5bZs8f8T5mTKVOgcWPzfZco4Xwp4MLnFyxn2u388NNPtLnjDoLDws5vCwoyO+Z5ctKdN8+3J11PT/aTJvn+l/3ixRgHDmTfnD5ihG9PCEXhZO/L5MafSYY/T3z+OLYf369xzz1w773+ucuyMA0jb/jZ7NmzjZCQEGPGjBnGH3/8YfTo0cMoW7ascfjwYbflf/75Z8NqtRoTJkwwtm7dagwZMsQIDg42Nm/e7NHxTp48aQDGyZMn8xz7xteXXHgz4CUfG19fkudjOTl3zjCqVTPs5qnA5WHHYhixsWY5X/ryS8OwWMzHhcfMWvflly4vycjIMObPn29kZGR4f9wlSzyqZ2PJEu+P4c7/6tnl/V74vr2s52zrJR+Pmy1/HdcwDOPLLw27xWLY3R3zEp8tXx03t59pnx67WjXn48bGOh3TJ/9/vDhuvjl3zvx/OmuW+dfLz5JX9eKjY+dKAR4zXz4rhUxuzt9+73PTtGlTrr/+et566y0A7HY7sbGxPPfccwwaNMil/IMPPsjp06f55ptvHOtuuOEGrr32WqZNm+ZSPj09nfT0dMfzlJQUYmNjOXr0qE/63ByPvIIo+wGCcK1GOxaSrNWocHKHz28Lt8ybh/Whh8zlC/4Jjf/9ArPNnm1m8D5mmTcPa79+WA4cOH/MatWwTZrk9niZmZkkJibSunVr7/sM2GyUqFULDh50+8vesFggJoZzO3f6/BdCftVzTvXi139fPxwXwD53LvbnnqPkBX2rsvts+UpuP9M+ZbNhWbnS8SvXaNHC5W7DPP//8eK4hV2+1UsRFgh1kpKSQsWKFQt/h+KMjAzCw8OZO3cuHTt2dKzv1q0bJ06c4KuvvnJ5zWWXXUa/fv3o06ePY93w4cOZP38+v/32m0v5ESNGMHLkSJf1s2bNIjw8PM/vIWXmn3SZPxjAKcHJulvq045jiex+dZ6P40706tXUmz7d6WSQVrEiW554gqRmzfLlmADYbFTYupWwf//lbLlyHKtTJ9+/GKNXr+b68eMB5/Ecsmp83cCB+fae/VXPgXZcwC+fLb8eV0Q8lpaWxsMPP1z4k5uDBw8SExPDqlWraHbBl+aLL77IsmXLWLNmjctrQkJC+Oijj+jcubNj3X//+19GjhzJ4cOHXcrnZ8tNlrWD5lNjSj+q2s9fWz1gjeXv3pNoMq6jT45xSYX8F5gvf00U5l/YueVxvfjr39cPxw2EX565pTpxT/XiKhDqJDctN4WiQ3F+Cg0NJTQ01GV9cHCwzz4AzSfdj21sJ9a/uZRtP/1C7dtu4NrnbiEmn0YodhIcDPHx+X+cPPJJfT/wgEtHOUvLlpQoiJN9PtVzjvXir39fP36ufPl/s7hQnbinenFVnOskN+/Lr8lNxYoVsVqtLi0uhw8fJioqyu1roqKiclW+oFhDrNR//ib21zpF/TtvwhpceFpPihWrtcDmJhERkaLJ3dhzBSYkJIRGjRqxePFixzq73c7ixYudLlNdqFmzZk7lARITEy9ZXkRERAKL3y9L9evXj27dutG4cWOaNGnClClTOH36NI899hgAXbt2JSYmhrFjxwLQu3dvbr75ZiZNmkS7du2YPXs2v/76K++9954/34aIiIgUEn5Pbh588EGOHDnCsGHDOHToENdeey2LFi2iSpUqAOzdu5egC+YKufHGG5k1axZDhgzhpZde4oorrmD+/PnUrVvXX29BREREChG/JzcAvXr1olevXm63LV261GXd/fffz/3335/PUYmIiEhR5Nc+NyIiIiK+puRGREREihUlNyIiIlKsKLkRERGRYkXJjYiIiBQrheJuqYKUNZVWSkqKz/edmZlJWloaKSkpxXb469xSnbinenGlOnGlOnFP9eIqEOok67ztyZSYAZfcpKamAhAbG+vnSERERCS3UlNTKVOmTLZl/DoruD/Y7XYOHjxI6dKlsVgsPt131ozj+/bt89mM40Wd6sQ91Ysr1Ykr1Yl7qhdXgVAnhmGQmppK1apVnQb3dSfgWm6CgoKoVq1avh4jMjKy2H64vKU6cU/14kp14kp14p7qxVVxr5OcWmyyqEOxiIiIFCtKbkRERKRYUXLjQ6GhoQwfPpzQ0FB/h1JoqE7cU724Up24Up24p3pxpTpxFnAdikVERKR4U8uNiIiIFCtKbkRERKRYUXIjIiIixYqSGxERESlWlNzkYPny5bRv356qVatisViYP3++03bDMBg2bBjR0dGULFmS+Ph4du7c6VTm+PHjdOnShcjISMqWLcsTTzzBqVOnCvBd+NbYsWO5/vrrKV26NJUrV6Zjx45s377dqczZs2d59tlnqVChAhEREdx7770cPnzYqczevXtp164d4eHhVK5cmQEDBnDu3LmCfCs+9c4771C/fn3HIFrNmjXju+++c2wPxDq50Lhx47BYLPTp08exLhDrZMSIEVgsFqfH1Vdf7dgeiHUCcODAAR555BEqVKhAyZIlqVevHr/++qtjeyB+18bFxbl8ViwWC88++ywQuJ8VjxiSrYULFxovv/yykZCQYADGvHnznLaPGzfOKFOmjDF//nzjt99+M+6++26jRo0axpkzZxxlbr/9dqNBgwbGL7/8YqxYscKoVauW0blz5wJ+J77Ttm1b48MPPzS2bNlibNq0ybjzzjuNyy67zDh16pSjzFNPPWXExsYaixcvNn799VfjhhtuMG688UbH9nPnzhl169Y14uPjjY0bNxoLFy40KlasaAwePNgfb8knFixYYHz77bfGjh07jO3btxsvvfSSERwcbGzZssUwjMCskyxr16414uLijPr16xu9e/d2rA/EOhk+fLhxzTXXGElJSY7HkSNHHNsDsU6OHz9uVK9e3ejevbuxZs0aY/fu3cb3339v7Nq1y1EmEL9rk5OTnT4niYmJBmAsWbLEMIzA/Kx4SslNLlyc3NjtdiMqKsqYOHGiY92JEyeM0NBQ47PPPjMMwzC2bt1qAMa6descZb777jvDYrEYBw4cKLDY81NycrIBGMuWLTMMw6yD4OBg44svvnCU2bZtmwEYq1evNgzDTBqDgoKMQ4cOOcq88847RmRkpJGenl6wbyAflStXzpg+fXpA10lqaqpxxRVXGImJicbNN9/sSG4CtU6GDx9uNGjQwO22QK2TgQMHGi1atLjkdn3Xmnr37m1cfvnlht1uD9jPiqd0WSoP9uzZw6FDh4iPj3esK1OmDE2bNmX16tUArF69mrJly9K4cWNHmfj4eIKCglizZk2Bx5wfTp48CUD58uUBWL9+PZmZmU71cvXVV3PZZZc51Uu9evWoUqWKo0zbtm1JSUnhjz/+KMDo84fNZmP27NmcPn2aZs2aBXSdPPvss7Rr187pvUNgf0527txJ1apVqVmzJl26dGHv3r1A4NbJggULaNy4Mffffz+VK1emYcOGvP/++47t+q6FjIwMPvnkEx5//HEsFkvAflY8peQmDw4dOgTg9MHJep617dChQ1SuXNlpe4kSJShfvryjTFFmt9vp06cPzZs3p27duoD5nkNCQihbtqxT2YvrxV29ZW0rqjZv3kxERAShoaE89dRTzJs3jzp16gRsncyePZsNGzYwduxYl22BWidNmzZl5syZLFq0iHfeeYc9e/bQsmVLUlNTA7ZOdu/ezTvvvMMVV1zB999/z9NPP83zzz/PRx99BOi7FmD+/PmcOHGC7t27A4H7/8dTATcruPjWs88+y5YtW1i5cqW/QykUrrrqKjZt2sTJkyeZO3cu3bp1Y9myZf4Oyy/27dtH7969SUxMJCwszN/hFBp33HGHY7l+/fo0bdqU6tWr8/nnn1OyZEk/RuY/drudxo0bM2bMGAAaNmzIli1bmDZtGt26dfNzdIXDBx98wB133EHVqlX9HUqRoJabPIiKigJw6Z1++PBhx7aoqCiSk5Odtp87d47jx487yhRVvXr14ptvvmHJkiVUq1bNsT4qKoqMjAxOnDjhVP7ienFXb1nbiqqQkBBq1apFo0aNGDt2LA0aNGDq1KkBWSfr168nOTmZ6667jhIlSlCiRAmWLVvGG2+8QYkSJahSpUrA1Yk7ZcuW5corr2TXrl0B+TkBiI6Opk6dOk7rateu7bhcF+jftf/88w8//vgjTz75pGNdoH5WPKXkJg9q1KhBVFQUixcvdqxLSUlhzZo1NGvWDIBmzZpx4sQJ1q9f7yjz008/Ybfbadq0aYHH7AuGYdCrVy/mzZvHTz/9RI0aNZy2N2rUiODgYKd62b59O3v37nWql82bNzt9GSUmJhIZGenyJVeU2e120tPTA7JOWrVqxebNm9m0aZPj0bhxY7p06eJYDrQ6cefUqVP89ddfREdHB+TnBKB58+Yuw0ns2LGD6tWrA4H7XZvlww8/pHLlyrRr186xLlA/Kx7zd4/mwi41NdXYuHGjsXHjRgMwJk+ebGzcuNH4559/DMMwb08sW7as8dVXXxm///670aFDB7e3JzZs2NBYs2aNsXLlSuOKK64o0rcnPv3000aZMmWMpUuXOt2mmJaW5ijz1FNPGZdddpnx008/Gb/++qvRrFkzo1mzZo7tWbcotmnTxti0aZOxaNEio1KlSkX6FsVBgwYZy5YtM/bs2WP8/vvvxqBBgwyLxWL88MMPhmEEZp1c7MK7pQwjMOukf//+xtKlS409e/YYP//8sxEfH29UrFjRSE5ONgwjMOtk7dq1RokSJYzRo0cbO3fuND799FMjPDzc+OSTTxxlAvG71jAMw2azGZdddpkxcOBAl22B+FnxlJKbHCxZssQAXB7dunUzDMO8RXHo0KFGlSpVjNDQUKNVq1bG9u3bnfZx7Ngxo3PnzkZERIQRGRlpPPbYY0Zqaqof3o1vuKsPwPjwww8dZc6cOWM888wzRrly5Yzw8HDjnnvuMZKSkpz28/fffxt33HGHUbJkSaNixYpG//79jczMzAJ+N77z+OOPG9WrVzdCQkKMSpUqGa1atXIkNoYRmHVysYuTm0CskwcffNCIjo42QkJCjJiYGOPBBx90Gs8lEOvEMAzj66+/NurWrWuEhoYaV199tfHee+85bQ/E71rDMIzvv//eAFzeq2EE7mfFExbDMAy/NBmJiIiI5AP1uREREZFiRcmNiIiIFCtKbkRERKRYUXIjIiIixYqSGxERESlWlNyIiIhIsaLkRkRERIoVJTciIiJSrCi5ERERkWJFyY2I+ET37t2xWCyMGzfOaf38+fOxWCx+iip/WCwW5s+f7+8wROQSlNyIiM+EhYUxfvx4/v33X3+HIiIBTMmNiPhMfHw8UVFRjB07NttyX375Jddccw2hoaHExcUxadKkXB3n2LFjdO7cmZiYGMLDw6lXrx6fffaZU5m5c+dSr149SpYsSYUKFYiPj+f06dMALF26lCZNmlCqVCnKli1L8+bN+eeffxyv/eqrr7juuusICwujZs2ajBw5knPnzgEQFxcHwD333IPFYnE8/+2337j11lspXbo0kZGRNGrUiF9//TVX70tEfEPJjYj4jNVqZcyYMbz55pvs37/fbZn169fzwAMP8NBDD7F582ZGjBjB0KFDmTlzpsfHOXv2LI0aNeLbb79ly5Yt9OzZk0cffZS1a9cCkJSUROfOnXn88cfZtm0bS5cupVOnThiGwblz5+jYsSM333wzv//+O6tXr6Znz56OS2crVqyga9eu9O7dm61bt/Luu+8yc+ZMRo8eDcC6desA+PDDD0lKSnI879KlC9WqVWPdunWsX7+eQYMGERwc7G1Vikhe+HlWchEpJrp162Z06NDBMAzDuOGGG4zHH3/cMAzDmDdvnnHhV83DDz9stG7d2um1AwYMMOrUqZOn47dr187o37+/YRiGsX79egMw/v77b5dyx44dMwBj6dKlbvfTqlUrY8yYMU7r/u///s+Ijo52PAeMefPmOZUpXbq0MXPmzDy9BxHxDbXciIjPjR8/no8++oht27a5bNu2bRvNmzd3Wte8eXN27tyJzWbzaP82m41XXnmFevXqUb58eSIiIvj+++/Zu3cvAA0aNKBVq1bUq1eP+++/n/fff9/RD6h8+fJ0796dtm3b0r59e6ZOnUpSUpJj37/99hujRo0iIiLC8ejRowdJSUmkpaVdMqZ+/frx5JNPEh8fz7hx4/jrr788ei8i4ntKbkTE52666Sbatm3L4MGD82X/EydOZOrUqQwcOJAlS5awadMm2rZtS0ZGBmBeHktMTOS7776jTp06vPnmm1x11VXs2bMHMC8prV69mhtvvJE5c+Zw5ZVX8ssvvwBw6tQpRo4cyaZNmxyPzZs3s3PnTsLCwi4Z04gRI/jjjz9o164dP/30E3Xq1GHevHn58v5FJHtKbkQkX4wbN46vv/6a1atXO62vXbs2P//8s9O6n3/+mSuvvBKr1erRvn/++Wc6dOjAI488QoMGDahZsyY7duxwKmOxWGjevDkjR45k48aNhISEOCUbDRs2ZPDgwaxatYq6desya9YsAK677jq2b99OrVq1XB5BQeZXZnBwsNtWpiuvvJK+ffvyww8/0KlTJz788EOP3o+I+JaSGxHJF/Xq1aNLly688cYbTuv79+/P4sWLeeWVV9ixYwcfffQRb731Fi+88IKjTKtWrXjrrbcuue8rrriCxMREVq1axbZt2/jPf/7D4cOHHdvXrFnDmDFj+PXXX9m7dy8JCQkcOXKE2rVrs2fPHgYPHszq1av5559/+OGHH9i5cye1a9cGYNiwYXz88ceMHDmSP/74g23btjF79myGDBni2H9cXByLFy/m0KFD/Pvvv5w5c4ZevXqxdOlS/vnnH37++WfWrVvn2KeIFDB/d/oRkeLhwg7FWfbs2WOEhIQYF3/VzJ0716hTp44RHBxsXHbZZcbEiROdtlevXt0YPnz4JY917Ngxo0OHDkZERIRRuXJlY8iQIUbXrl0dx9+6davRtm1bo1KlSkZoaKhx5ZVXGm+++aZhGIZx6NAho2PHjkZ0dLQREhJiVK9e3Rg2bJhhs9kc+1+0aJFx4403GiVLljQiIyONJk2aGO+9955j+4IFC4xatWoZJUqUMKpXr26kp6cbDz30kBEbG2uEhIQYVatWNXr16mWcOXPGi5oUkbyyGIZh+DvBEhEREfEVXZYSERGRYkXJjYiIiBQrSm5ERESkWFFyIyIiIsWKkhsREREpVpTciIiISLGi5EZERESKFSU3IiIiUqwouREREZFiRcmNiIiIFCtKbkRERKRY+X8j9yqy7uH8JAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fix, ax = plt.subplots()\n",
    "ax.plot(problem_dims, time_sigma, c=\"blue\", marker=\"o\", label=\"sigma\")\n",
    "ax.plot(problem_dims, time_factor, c=\"red\", marker=\"o\", label=\"factor\")\n",
    "ax.legend()\n",
    "plt.xlabel(\"No. assets\")\n",
    "plt.ylabel(\"Solution time (sec.)\")\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca4030a2",
   "metadata": {},
   "source": [
    "From the data, it is clear that the solve times of (2) are much smaller both in absolute and asymptotic terms.  The data may appear somewhat nonsmooth over the problem dimensions.  These are artifacts from multi-threading and different code paths in the computational kernels used for dense matrix multiplication."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab923db0",
   "metadata": {},
   "source": [
    "## Takeaways\n",
    "\n",
    "- Incorporating explicit factors to model risk in a gurobipy model is straightforward.\n",
    "- Using factor risk models instead of the derived full covariance matrix can greatly improve solution time."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79432dfa",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] Menchero, J., Orr, D., Wang, J.: The Barra US equity model (USE4), methodology notes. English,\n",
    "MSCI (May (2011)\n",
    "\n",
    "[2] Gérard Cornuéjols, Javier Peña, and Reha Tütüncü. Optimization Methods in Finance. Cambridge University Press, 2 edition, 2018. doi:10.1017/9781107297340."
   ]
  }
 ],
 "metadata": {
  "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}