{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 4: tuning parameters for iDDN\n",
    "\n",
    "This tutorial shows how to tune the two hyperparameters using data driven cross validation based methods.\n",
    "Our approach is to first search for $\\lambda_1$, then based on the chosen value of $\\lambda_1$ to search for $\\lambda_2$.\n",
    "\n",
    "Although we also provide functions to jointly search $\\lambda_1$ and $\\lambda_2$ on a grid, due to computational reason, it is often not efficient when the network is large, so we will not discuss this option here. However, as you will see, it is straightforward to do that if desired. In addition, we also provide methods purely based on statistic and probability theory, which is very easy to calculate, but may not always be accurate enough. We will briefly talk about it.\n",
    "\n",
    "Again, it is important to note that the best choice of $\\lambda_1$ and $\\lambda_2$ often relies on domain knowledge. \n",
    "Those provided by cross validation may not reflect the underlying biological process. \n",
    "Therefore, we advise user to choose a range of parameters and plot networks to choose good parameters.\n",
    "\n",
    "Compared with the similar tutorial in the DDN3.0 package, now we 1) support iDDN (multi-omics with constraints), 2) support parallel computing for cross validation, and 3) include some tricks to handle large problems.\n",
    "\n",
    "We first import the necessary modules.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:09:33.045515Z",
     "start_time": "2024-08-26T16:09:31.675905Z"
    }
   },
   "source": [
    "import numpy as np\n",
    "from iddn_data import load_data  # Load example data and images\n",
    "from iddn import parameter_tuning_iddn as parameter_tuning  # parameter tuning tools\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ],
   "outputs": [],
   "execution_count": 1
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "We load the data."
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "We use the same example as before. The parameter tuning supports parallel computing.\n",
    "As our example is small, one core (`cores=1`) is good enough. Using multiple core will not have clear benefits here. "
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:09:33.125325Z",
     "start_time": "2024-08-26T16:09:33.050501Z"
    }
   },
   "cell_type": "code",
   "source": [
    "example = load_data.load_example()\n",
    "dat1 = example[\"dat1\"]\n",
    "dat2 = example[\"dat2\"]\n",
    "dep_mat = example[\"dep_mat\"]\n",
    "\n",
    "cores = 1"
   ],
   "outputs": [],
   "execution_count": 2
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Step 1: search $\\lambda_1$ using cross validation\n",
    "\n",
    "We first keep $\\lambda_2$ fixed and scan $\\lambda_1$ to find the optimal value using cross validation.\n",
    "This allows us to efficiently work on larger networks (compared to using grid search for two parameters).\n",
    "\n",
    "Specifically, for each $\\lambda_1$, certain part (`ratio_val`) of the data (training set) is used to estimate the networks, and then for each node, the regression coefficients are calculated to use the values of all other nodes to estimate the value of this node.\n",
    "These regression models are then applied on the remaining part of the data (test set) to estimate each node using all other nodes, and we calculate the unexplained variance as the error measure.\n",
    "\n",
    "Therefore, we have an error value for each $\\lambda_1$. Since we repeat the cross validation for several time (`n_cv`), we have an associated standard error. This allows us to pick the optimal $\\lambda_1$, and optionally considers the uncertainty of estimation.\n",
    "\n",
    "Here we keep the $\\lambda_2$ value as 0.05, which is a typical value used in many cases.\n",
    "We search $\\lambda_1$ from 0.02 to 0.4, with a step size of 0.02.\n"
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:09:34.170582Z",
     "start_time": "2024-08-26T16:09:33.396653Z"
    }
   },
   "cell_type": "code",
   "source": "l1_list = np.arange(0.02, 0.401, 0.02)",
   "outputs": [],
   "execution_count": 3
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "Here we call the cross validation function. This function actually allows doing grid search if we input a list of $\\lambda_1$ as well as a list of $\\lambda_2$. We need to specify \n",
    "- `dat1` and `dat2`: the input data\n",
    "- `n_cv`: the repeats for cross validation\n",
    "- `ratio_val`: the ratio of data left out for validation (0.5 means 50%)\n",
    "- `lambda1_lst` and `lambda2_lst`: the list of parameters to search\n",
    "- `dep_mat`: the constraint matrix of iDDN\n",
    "- `cores`: the number of cores to use"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:09:57.320359Z",
     "start_time": "2024-08-26T16:09:34.774688Z"
    }
   },
   "source": [
    "val_err_l1 = parameter_tuning.cv_2d(\n",
    "    dat1,\n",
    "    dat2,\n",
    "    n_cv=5,\n",
    "    ratio_val=0.5, \n",
    "    lambda1_lst=l1_list,\n",
    "    lambda2_lst=[0.05],\n",
    "    dep_mat=dep_mat,\n",
    "    cores=cores,\n",
    ")"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.02 0.05 s0\n",
      "0 0.04 0.05 s0\n",
      "0 0.06 0.05 s0\n",
      "0 0.08 0.05 s0\n",
      "0 0.1 0.05 s0\n",
      "0 0.12000000000000001 0.05 s0\n",
      "0 0.13999999999999999 0.05 s0\n",
      "0 0.16 0.05 s0\n",
      "0 0.18 0.05 s0\n",
      "0 0.19999999999999998 0.05 s0\n",
      "0 0.22 0.05 s0\n",
      "0 0.24 0.05 s0\n",
      "0 0.26 0.05 s0\n",
      "0 0.28 0.05 s0\n",
      "0 0.30000000000000004 0.05 s0\n",
      "0 0.32 0.05 s0\n",
      "0 0.34 0.05 s0\n",
      "0 0.36000000000000004 0.05 s0\n",
      "0 0.38 0.05 s0\n",
      "0 0.4 0.05 s0\n",
      "1 0.02 0.05 s0\n",
      "1 0.04 0.05 s0\n",
      "1 0.06 0.05 s0\n",
      "1 0.08 0.05 s0\n",
      "1 0.1 0.05 s0\n",
      "1 0.12000000000000001 0.05 s0\n",
      "1 0.13999999999999999 0.05 s0\n",
      "1 0.16 0.05 s0\n",
      "1 0.18 0.05 s0\n",
      "1 0.19999999999999998 0.05 s0\n",
      "1 0.22 0.05 s0\n",
      "1 0.24 0.05 s0\n",
      "1 0.26 0.05 s0\n",
      "1 0.28 0.05 s0\n",
      "1 0.30000000000000004 0.05 s0\n",
      "1 0.32 0.05 s0\n",
      "1 0.34 0.05 s0\n",
      "1 0.36000000000000004 0.05 s0\n",
      "1 0.38 0.05 s0\n",
      "1 0.4 0.05 s0\n",
      "2 0.02 0.05 s0\n",
      "2 0.04 0.05 s0\n",
      "2 0.06 0.05 s0\n",
      "2 0.08 0.05 s0\n",
      "2 0.1 0.05 s0\n",
      "2 0.12000000000000001 0.05 s0\n",
      "2 0.13999999999999999 0.05 s0\n",
      "2 0.16 0.05 s0\n",
      "2 0.18 0.05 s0\n",
      "2 0.19999999999999998 0.05 s0\n",
      "2 0.22 0.05 s0\n",
      "2 0.24 0.05 s0\n",
      "2 0.26 0.05 s0\n",
      "2 0.28 0.05 s0\n",
      "2 0.30000000000000004 0.05 s0\n",
      "2 0.32 0.05 s0\n",
      "2 0.34 0.05 s0\n",
      "2 0.36000000000000004 0.05 s0\n",
      "2 0.38 0.05 s0\n",
      "2 0.4 0.05 s0\n",
      "3 0.02 0.05 s0\n",
      "3 0.04 0.05 s0\n",
      "3 0.06 0.05 s0\n",
      "3 0.08 0.05 s0\n",
      "3 0.1 0.05 s0\n",
      "3 0.12000000000000001 0.05 s0\n",
      "3 0.13999999999999999 0.05 s0\n",
      "3 0.16 0.05 s0\n",
      "3 0.18 0.05 s0\n",
      "3 0.19999999999999998 0.05 s0\n",
      "3 0.22 0.05 s0\n",
      "3 0.24 0.05 s0\n",
      "3 0.26 0.05 s0\n",
      "3 0.28 0.05 s0\n",
      "3 0.30000000000000004 0.05 s0\n",
      "3 0.32 0.05 s0\n",
      "3 0.34 0.05 s0\n",
      "3 0.36000000000000004 0.05 s0\n",
      "3 0.38 0.05 s0\n",
      "3 0.4 0.05 s0\n",
      "4 0.02 0.05 s0\n",
      "4 0.04 0.05 s0\n",
      "4 0.06 0.05 s0\n",
      "4 0.08 0.05 s0\n",
      "4 0.1 0.05 s0\n",
      "4 0.12000000000000001 0.05 s0\n",
      "4 0.13999999999999999 0.05 s0\n",
      "4 0.16 0.05 s0\n",
      "4 0.18 0.05 s0\n",
      "4 0.19999999999999998 0.05 s0\n",
      "4 0.22 0.05 s0\n",
      "4 0.24 0.05 s0\n",
      "4 0.26 0.05 s0\n",
      "4 0.28 0.05 s0\n",
      "4 0.30000000000000004 0.05 s0\n",
      "4 0.32 0.05 s0\n",
      "4 0.34 0.05 s0\n",
      "4 0.36000000000000004 0.05 s0\n",
      "4 0.38 0.05 s0\n",
      "4 0.4 0.05 s0\n"
     ]
    }
   ],
   "execution_count": 4
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "The output (val_err_l1) is of shape (repeats of CV, number of lambda1, number of lambda2), and each element is the unexplained variance in the validation set.\n",
    "We can plot the curve with respect to the $\\lambda_1$ values. Here we need to `squeeze()` the `val_err_l1`, as the drawing function need the input to have shape (repeats of CV, number of lambda1 or lambda2)."
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:09:57.603966Z",
     "start_time": "2024-08-26T16:09:57.384187Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# val_avg = np.mean(val_err_l1_1, axis=0)\n",
    "# val_std = np.std(val_err_l1_1, axis=0)\n",
    "parameter_tuning.plot_error_1d(val_err_l1.squeeze(), lambda_lst=l1_list)"
   ],
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABN40lEQVR4nO3deVhU9f4H8PfMwDCsg+yLKIsLruBKansUaqV57aZmLrRqWreoa1qmZbe0bnk1tTK7pVdzTbPtRhk3K3dFwAVxAQRUdpSBwQFmzvn9gYzxE5fBmTmzvF/PM8+jh+85fL4Oy9tzvotMFEURRERERDZMLnUBRERERNfDwEJEREQ2j4GFiIiIbB4DCxEREdk8BhYiIiKyeQwsREREZPMYWIiIiMjmMbAQERGRzXORugBzEQQB586dg7e3N2QymdTlEBER0Q0QRRE1NTUICwuDXH71+ygOE1jOnTuHiIgIqcsgIiKiNigqKkL79u2v+nGHCSze3t4Amjrs4+Nj9uu//UM21u0rwuO3RiHl3i5mvz4REZEz0mg0iIiIMP4evxqHCSzNj4F8fHwsEli6dQiGPKsS57SwyPWJiIic2fWGc3DQ7Q2KDvQCAORVaCWuhIiIyPkwsNygmKCmwFJQqYXeIEhcDRERkXNhYLlBoT4qqFzlaDSIKDp/UepyiIiInAoDyw2Sy2WICrj0WKi8VuJqiIiInAsDiwmiAz0BAHnlHMdCRERkTQwsJoi5NPA2l3dYiIiIrIqBxQQxvMNCREQkCQYWE0Q3j2Gp4B0WIiIia2JgMUHUpTssFbUNqK5rlLgaIiIi58HAYgIvNxeE+KgAALm8y0JERGQ1DCwm4kwhIiIi62NgMdHlwMI7LERERNbCwGIi48Bb3mEhIiKyGgYWEzXvKcS1WIiIiKyHgcVE0QFNj4QKKutgEESJqyEiInIODCwmCvd1h5uLHA0GAWfO10ldDhERkVNgYDFR0yaInClERERkTQwsbcA9hYiIiKyLgaUNmqc25/IOCxERkVUwsLQB12IhIiKyLgaWNri8CSLvsBAREVkDA0sbNN9hKa+ph0bHTRCJiIgsjYGlDbxVrgjydgPAmUJERETWYHJg+f333/Hggw8iLCwMMpkMW7duve4527dvR9++feHm5oZOnTph5cqVV7RZtmwZIiMjoVKpkJCQgH379plamlVxHAsREZH1mBxYtFot4uLisGzZshtqn5+fj/vvvx933XUXMjMz8cILL+DJJ5/ETz/9ZGyzYcMGpKSkYO7cuTh48CDi4uKQlJSEsrIyU8uzmuhA7ilERERkLS6mnjBs2DAMGzbshtt/8skniIqKwgcffAAA6NatG3bs2IF//etfSEpKAgAsXLgQTz31FJKTk43n/PDDD/j8888xc+ZMU0u0Cq7FQkREZD0WH8Oye/duJCYmtjiWlJSE3bt3AwAaGhqQnp7eoo1cLkdiYqKxTWvq6+uh0WhavKzp8iMh3mEhIiKyNIsHlpKSEgQHB7c4FhwcDI1Gg4sXL6KiogIGg6HVNiUlJVe97vz586FWq42viIgIi9R/NTGXpjbnV2q5CSIREZGF2e0soVmzZqG6utr4KioqsurnD2/nDqWLHA16AecuXLTq5yYiInI2Jo9hMVVISAhKS0tbHCstLYWPjw/c3d2hUCigUChabRMSEnLV67q5ucHNzc0iNd8IhVyGKH9PHC+twanyWkT4eUhWCxERkaOz+B2WQYMGIS0trcWxbdu2YdCgQQAApVKJfv36tWgjCALS0tKMbWwVx7EQERFZh8mBpba2FpmZmcjMzATQNG05MzMThYWFAJoe1UycONHYfsqUKcjLy8OMGTOQk5ODjz76CBs3bsSLL75obJOSkoIVK1Zg1apVOHbsGKZOnQqtVmucNWSruBYLERGRdZj8SOjAgQO46667jH9PSUkBAEyaNAkrV65EcXGxMbwAQFRUFH744Qe8+OKLWLx4Mdq3b4/PPvvMOKUZAMaMGYPy8nLMmTMHJSUliI+PR2pq6hUDcW2NcU8h3mEhIiKyKJkoig4xxUWj0UCtVqO6uho+Pj5W+ZyZRRfw0LKdCPJ2w77XEq9/AhEREbVwo7+/7XaWkC1ofiRUVlOPGm6CSEREZDEMLDfBR+WKAK+mmUr5FXwsREREZCkMLDeJM4WIiIgsj4HlJnFPISIiIstjYLlJMbzDQkREZHEMLDep+ZEQ77AQERFZDgPLTWpeiyW/QguBmyASERFZBAPLTWrfzh1KhRz1egFnuQkiERGRRTCw3CQXhRwd/Zs2Pszj1GYiIiKLYGAxA+4pREREZFkMLGYQHcg9hYiIiCyJgcUMuBYLERGRZTGwmAFXuyUiIrIsBhYziLk0tblEo4O2Xi9xNURERI6HgcUM1B6u8PdUAuAmiERERJbAwGImHMdCRERkOQwsZnJ5iX7eYSEiIjI3BhYz4VosRERElsPAYibNewpxphAREZH5MbCYSUzQpcBSUctNEImIiMyMgcVMItq5w1Uhg65RQLFGJ3U5REREDoWBxUxcFHJ08Lu0CSLHsRAREZkVA4sZcU8hIiIiy2BgMSOuxUJERGQZDCxmxD2FiIiILIOBxYxiuBYLERGRRTCwmFHzWiznqnWoa+AmiERERObCwGJG7TyV8Lu0CSIfCxEREZkPA4uZRQdceizEXZuJiIjMhoHFzLinEBERkfkxsJgZ12IhIiIyPwYWM+NaLERERObHwGJmzY+E8iu0EEVugkhERGQODCxm1sHPAy5yGeoaDCjhJohERERmwcBiZq4tNkHkOBYiIiJzYGCxgGiOYyEiIjIrBhYLiOGeQkRERGbFwGIBzQNveYeFiIjIPBhYLIBrsRAREZkXA4sFNK/FcvbCRVxsMEhcDRERkf1jYLEAP08lfD1cATStx0JEREQ3h4HFQi5vgshxLERERDeLgcVCOI6FiIjIfBhYLIR7ChEREZkPA4uFRHMtFiIiIrNhYLGQy4vH1XITRCIiopvEwGIhHfw8oZDLoG0woKymXupyiIiIbs66dZJ+egYWC1G6XN4EMbeM41iIiMg+6RoN2HSgCK//WiRpHS6SfnYHFx3gifwKLXIrtBjcKUDqcoiIiG5YqUaHNXsKsHZvISq1DYBfDzxarEG3UB9J6mFgsaDoQE+k5TSNYyEiIrJ1oijiYOEFrNx1Gj8eLoZeaBqDGaZWYcKJ7QhV3ytZbQwsFsS1WIiIyB7U6w347+FirNx5Gllnqo3HB0b6YfKQSNzXPRguoxYDHkrJamRgsSCuxUJERLasrEaHtXsLsWZPISpqmyaIKF3kGBkXhkmDI9EzXH258bhxElXZhIHFgprXYjl74SJ0jQaoXBUSV0RERARkFTU99vn+0Dk0Gpoe+wT7uGHCLR0xbmAH+Hu5XXkSA4vj8vdUwkflAo1Oj9OVWsSGSDNQiYiIqNEg4McjJVi5Mx8HCy8Yj/ft4IvJQ6IwrGcIXBW2O3mYgcWCZDIZogO9kFl0AXnlDCxERGR9lbX1WLevEKv3FKBU0/TYx1Uhw4O9mx77xEX4SlvgDWpTlFq2bBkiIyOhUqmQkJCAffv2XbVtY2Mj5s2bh5iYGKhUKsTFxSE1NbVFG4PBgNdffx1RUVFwd3dHTEwM3nrrLYdYIdY4joVrsRARkRUdOVuNlzdlYdCC/+H9n0+gVFOPAC83vJDYGTtn3o2FY+LtJqwAbbjDsmHDBqSkpOCTTz5BQkICFi1ahKSkJBw/fhxBQUFXtJ89ezbWrFmDFStWIDY2Fj/99BNGjRqFXbt2oU+fPgCAd999Fx9//DFWrVqFHj164MCBA0hOToZarcbzzz9/872UkHFPoQrOFCIiIgtatw76R8bg5+xSfLEzH/tPnzd+KK69GslDojC8VyiULrb72OdaZKKJtzESEhIwYMAALF26FAAgCAIiIiLw3HPPYebMmVe0DwsLw2uvvYZp06YZj40ePRru7u5Ys2YNAOCBBx5AcHAw/v3vf1+1zfVoNBqo1WpUV1fDx8d2Hr2kHinGlDUHEddejW+m3yp1OURE5ICq6xrx5bS3sKbLHThXrQMAuMhlGN4rFJOHRKJPhC9kMpnEVbbuRn9/m3SHpaGhAenp6Zg1a5bxmFwuR2JiInbv3t3qOfX19VCpVC2Oubu7Y8eOHca/Dx48GJ9++ilOnDiBLl26ICsrCzt27MDChQuvWkt9fT3q6y/v0aPRaEzpitX8eS0WURRt9guGiIjsj67RgJW7TuOjX09BE5wAVOvg76nEowkd8NgtHRHso7r+ReyESYGloqICBoMBwcHBLY4HBwcjJyen1XOSkpKwcOFC3H777YiJiUFaWhq2bNkCg8FgbDNz5kxoNBrExsZCoVDAYDDg7bffxvjx469ay/z58/Hmm2+aUr4kOvp7QC4Daur1KK+pR5ADffEQEZE09AYBX6WfwaJfTqJE03RHpYuuCk9PuAsP9A51yGU0LP4ga/HixejcuTNiY2OhVCoxffp0JCcnQy6//Kk3btyIL7/8EmvXrsXBgwexatUqvP/++1i1atVVrztr1ixUV1cbX0VF0m7KdDVuLgpENG+CyBVviYjoJoiiiJ+OliBp0e+YueUwSjQ6hPu644O/xuHHvK/wcL/2DhlWABPvsAQEBEChUKC0tLTF8dLSUoSEhLR6TmBgILZu3QqdTofKykqEhYVh5syZiI6ONrb5+9//jpkzZ2Ls2LEAgF69eqGgoADz58/HpEmTWr2um5sb3NxaWdjGBkUHeKKgsg55FbUYFOMvdTlERGSH9uZV4t3UHOMaKr4erph+Vyc8dkvHppAybqy0BVqYSXdYlEol+vXrh7S0NOMxQRCQlpaGQYMGXfNclUqF8PBw6PV6bN68GSNHjjR+rK6ursUdFwBQKBQQBMGU8mwW9xQiIqK2yinR4PGV+zHm0z04WHgBKlc5pt/VCb/PuAtP3hZ9+Y6KxCvRWprJ05pTUlIwadIk9O/fHwMHDsSiRYug1WqRnJwMAJg4cSLCw8Mxf/58AMDevXtx9uxZxMfH4+zZs3jjjTcgCAJmzJhhvOaDDz6It99+Gx06dECPHj2QkZGBhQsX4vHHHzdTN6XVPLWZewoREdGNOnO+Dgu3ncDXGWchioBCLsPYARH42z2dnXI8pMmBZcyYMSgvL8ecOXNQUlKC+Ph4pKamGgfiFhYWtrhbotPpMHv2bOTl5cHLywvDhw/H6tWr4evra2yzZMkSvP7663j22WdRVlaGsLAwPPPMM5gzZ87N99AGxPAOCxER3aAqbQOW/XoKq3cXoMHQ9KTh/l6heOm+LsY79s7I5HVYbJWtrsMCNO2GOfDtNMhlwLG3hsLNxTEHRBERUdvVNejx7z/y8enveaip1wMABsf445WhsXa1Iq2pLLIOC7VNoJcbvN1cUFOvR0FlHboEe0tdEhER2YhGg4D1+4vwYdpJlNc0rS/WPdQHM4fF4rbOAVy/6xIGFito2gTRE1lnqpFbVsvAQkREEAQR/z1SjPd/Oo7TlXUAgA5+Hnjpvi54sHcY5HIGlT9jYLGSmEAvZJ2p5p5CRESEnacqsODHHBw+Ww0ACPBS4rm7O2PcwA52u9ePpTGwWAlnChER0ZGz1Xg3NQd/nKwAAHgqFXj69hg8cVsUvNz4K/la+K9jJVyLhYjIeZXV6PCP74/h26xzAABXhQzjEzpi+t2dEOBlH4ugSo2BxUr+fIeFmyASETmPXbkVeH5dJipqmwbUjowPw0v3dkUHfw+JK7MvfFBmJZH+npDJgBqdHhW1DVKXQ0REFiYIIpakncRjn+1FRW09ugZ74/u8r7B4bB+GlTbgHRYrUbkq0L6dO4qqLiKvvBaB3rwFSETkqCpr6/HChkzjWJVH+rfHmyN6wj3krMSV2S8GFiuKDvBqCiwVWiREcxNEIiJHtP90FZ5bm4ESjQ4qVzn+8VAvPNyvfdMHHXy/H0tiYLGi6EBP/HaiHLllnClERORoBEHEp3/k4Z8/HYdBEBET6ImPxvdD1xCuvWUODCxWZNxTiGuxEBE5lPPaBry0KQv/yykDADwUH4a3R/WCJ6cqmw3/Ja2oeaZQHtdiISJyGBmF5zF9bQbOXrgIpYscb47ogbEDIjgb1MwYWKyo+Q5L0fmLqNcbuAkiEZEdE0URn+88jQU/HkOjQUSkvweWje+LHmFqqUtzSAwsVhTk7QZPpQLaBgMKK+vQmXsKERHZpeqLjZjxVRZ+OloKALi/VygWjO4Fb5WrxJU5LgYWK5LJZIgJ8sKhM9XILdcysBAR2aHDZ6rx7Np0FFVdhKtChtcf6I4Jt3TkIyALY2CxsugATxw6U428Co5jISKyJ6IoYs2eArz1/TE0GAS0b+eOj8b3Re/2vlKX5hQYWKyMewoREdmf2no9Zm4+hO8PFQMA7u0ejPcfjoPag4+ArIWBxcq4azMRkX05VqzBtC8PIq9CCxe5DDOHxeKJW6P4CMjKGFisLOZPd1i4CSIRke0SRREb9hdh7rdHUa8XEKZWYcmjfdGvYzupS3NKDCxWFhXQtAli9cVGVGkb4M9txYmIbE5dgx6zvz6CLRlNe//c1TUQCx+JRztPpcSVOS8GFitTuSoQpnbH2QtNewoxsBAR2ZaTpTV49suDOFlWC4Vchpfu64Ipt8dALucdcSnJpS7AGRnHsXBPISIim7I5/QxGLN2Jk2W1CPJ2w9onE/DsnZ0YVmwA77BIICbQC3+crOCeQkRENqLRIOD1rUewfn8RAODWTgFYNDYeAbwLbjMYWCQQwz2FiIhsRqNBwPPrMvDjkRLIRBEv3NsV0+/uBAXvqtgUPhKSANdiISKyDY0GAX9b3xRWlAo5VhT9hL8ldmZYsUG8wyKB5jEsBVV1aNALULowNxIRWVujQcAL6zPx38NNYeWTCX1xd8YFqcuiq+BvSgmE+KjgoVTAIIgorKqTuhwiIqejvxRWfjhcDFeFDB8/1hd3xwYD48ZJXRpdBQOLBGQymfEuC8exEBFZl94g4G9/CiufPNYP93QLlrosug4GFolEB1wax8KZQkREVqM3CPjbhj/dWRnPsGIvGFgkwrVYiIisS28Q8MKGTPxwqCmsfDS+HxK7M6zYCwYWiRj3FOIdFiIii9MbBLy4MQvfXworyx7ti3sZVuwKA4tEOIaFiMg69AYBKRuz8F3WObjIm8LKfT1CpC6LTMTAIpGogKbAcr6uaRNEIiIyP71BwEubsvDtpbDy0XiGFXvFwCIRD6ULwtQqAEDeuq3SFkNE5IAMgoiXNmXhm8xLd1YYVuwaA4uEYoIujWP5fb/ElRARORaDIOKljZnGsLL00b5IYlixawwsEoq+9FgoV+krbSFERA7EIIh4eVMWthrDSh8M7cmwYu8YWCRk3FNIqZa4EiIix2AQRPx9Uxa+zjgLhVyGJeP6YGjPUKnLIjNgYJGQcS0Wt3YSV0JEZP8Mgoi/f5WFLZfCytJxfTCsF8OKo2BgkVC3UB8AQJ6bL2cKERHdBIMgYsZXh7Dl4OU7KwwrjoWBRUIBXm6IDfEGAOzOrZS4GiIi+2QQRLyy+RA2HzwDhVyGD8f2wXCGFYfDwCKxwTEBAIAdpyokroSIyP40h5Wv0i+Hlft7M6w4IgYWid3a2R8AsJOBhYjIJIIgYuafwsrisfEMKw6MgUViA6P84SIaUFhVh6KqOqnLISKyC8KlOyubLoWVRWPi8UDvMKnLIgtiYJGYl5sL+niJAHiXhYjoRgiCiJlbmsKKXAYsGhOPB+MYVhwdA4sNGHxLNwAcx0JEdD2CIGLWlsPYeOBSWBnbh2HFSTCw2IBbOzcNvN2VWwlBECWuhojINgmCiFe/PowNB4oglwH/GhOPEQwrToOBxQbER/jCU6lAlbYBx0o0UpdDRGRzRFHEnG+PYP3+y2FlZHy41GWRFTGw2ABXhRwJ0U2zhXad4nosRET/3+c7T2PNnkLIIWLhIwwrzoiBxUYMjmkKLBzHQkTU0m8nyvH2D9kAgFfD6vFQH4YVZ+QidQHUpHkcy778KtTrDXBzUUhcERGR9HLLazF97UEIIvBI//Z4YnRvqUsiifAOi43oGuyNAC8lLjYakFF4QepyiIgkV13XiKdWHUCNTo/+HdvhrYd6QiaTSV0WSYSBxUbIZDIM6XRpthAfCxGRk9MbBExfdxB5FVqE+7rjkwn9eOfZyTGw2JAh3FeIiAgA8PZ/j+GPkxVwd1Xg04n9EODlJnVJJDEGFhsy5NI4lqwz1dDoGiWuhohIGhv2F+KLnacBAP8aE4ceYWppCyKb0KbAsmzZMkRGRkKlUiEhIQH79u27atvGxkbMmzcPMTExUKlUiIuLQ2pq6hXtzp49i8ceewz+/v5wd3dHr169cODAgbaUZ7fCfd0RFeAJgyBib16V1OUQEVnd/tNVmL31CADgxcQuGNqTmxlSE5MDy4YNG5CSkoK5c+fi4MGDiIuLQ1JSEsrKylptP3v2bCxfvhxLlixBdnY2pkyZglGjRiEjI8PY5vz58xgyZAhcXV3x448/Ijs7Gx988AHatWvX9p7ZqSGduHszETmnM+frMGV1OhoNIu7vHYrn7+kkdUlkQ2SiKJq0FnxCQgIGDBiApUuXAgAEQUBERASee+45zJw584r2YWFheO211zBt2jTjsdGjR8Pd3R1r1qwBAMycORM7d+7EH3/80eaOaDQaqNVqVFdXw8fHp83XkVrqkWJMWXMQnYO8sC3lDqnLISKyCm29HqM/3oWckhr0DPfBpmcGw13JQbbO4EZ/f5t0h6WhoQHp6elITEy8fAG5HImJidi9e3er59TX10OlUrU45u7ujh07dhj//u2336J///7461//iqCgIPTp0wcrVqwwpTSHcUu0P2Qy4GRZLUo1OqnLISKyOEEQ8eKGTOSU1CDAyw2fTujPsEJXMCmwVFRUwGAwIDg4uMXx4OBglJSUtHpOUlISFi5ciJMnT0IQBGzbtg1btmxBcXGxsU1eXh4+/vhjdO7cGT/99BOmTp2K559/HqtWrbpqLfX19dBoNC1ejsDXQ4le4U0DzPhYiIicwb9+OYGfs0uhVMjx6cR+CPN1l7okskEWnyW0ePFidO7cGbGxsVAqlZg+fTqSk5Mhl1/+1IIgoG/fvnjnnXfQp08fPP3003jqqafwySefXPW68+fPh1qtNr4iIiIs3RWraV6PhdObicjRfZd1Dkv+dwoAMP8vvdC3g/ONXaQbY1JgCQgIgEKhQGlpaYvjpaWlCAkJafWcwMBAbN26FVqtFgUFBcjJyYGXlxeio6ONbUJDQ9G9e/cW53Xr1g2FhYVXrWXWrFmorq42voqKikzpik271biAXCVMHGJERGQ3Dp+pxsubsgAAT98ejdH92ktcEdkykwKLUqlEv379kJaWZjwmCALS0tIwaNCga56rUqkQHh4OvV6PzZs3Y+TIkcaPDRkyBMePH2/R/sSJE+jYseNVr+fm5gYfH58WL0fRr2M7KF3kKNHokFuulbocIiKzK9Po8NR/DqBeL+CuroF4ZWis1CWRjTP5kVBKSgpWrFiBVatW4dixY5g6dSq0Wi2Sk5MBABMnTsSsWbOM7ffu3YstW7YgLy8Pf/zxB4YOHQpBEDBjxgxjmxdffBF79uzBO++8g1OnTmHt2rX49NNPW8wsciYqVwUGRDbdFuU4FiJyNLpGA55enY4SjQ6dgrzw4bg+UMi5RxBdm8m7NY8ZMwbl5eWYM2cOSkpKEB8fj9TUVONA3MLCwhbjU3Q6HWbPno28vDx4eXlh+PDhWL16NXx9fY1tBgwYgK+//hqzZs3CvHnzEBUVhUWLFmH8+PE330M7NaRTAHaeqsSOUxWYNDhS6nKIiMxCFEXM2nIYmUUXoHZ3xWcT+8Nb5Sp1WWQHTF6HxVY5yjoszQ6duYARS3fCW+WCjNfvhYuCuygQkf375LdcLPgxBwq5DKsfH4jBl8bskfOyyDosZD09wtTwUbmgRqfH4bPVUpdDRHTT0o6V4t3UHADA3Ae7M6yQSRhYbJRCLsPgS7s3cxwLEdm7E6U1+Nv6TIgi8GhCB0y45eqTKohaw8Biw5p3b+Z6LERkz85rG/DkqgOordfjlmg/vDmiB2QyDrIl0zCw2LDm9VgOFlzAxQaDxNUQEZmu0SDg2S8PorCqDhF+7vhofD+4ckwetQG/amxYpL8HwtQqNBgE7D9dJXU5REQmm/ddNnbnVcJTqcC/Jw2An6dS6pLITjGw2DCZTGZcpp/jWIjI3qzeU4DVewogkwGLx/ZBl2BvqUsiO8bAYuNu5TgWIrJDu3Ir8Ma3RwEAf0/qisTuwdc5g+jaGFhsXPNMoexiDaq0DRJXQ0R0fQWVWjz75UEYBBEPxYdh6h0xUpdEDoCBxcYFeruha7A3RBHYnVspdTlERNdUo2vEk6sO4EJdI+IifLFgdG/OCCKzYGCxA83jWPhYiIhsmUEQ8cL6TJwsq0WwjxtWTOgHlatC6rLIQTCw2IFbO/sD4MBbIrJtH6adRFpOGdxc5FgxsT+CfFRSl0QOhIHFDgyM8oeLXIbCqjoUVdVJXQ4R0RWyz2mw7NdTAIB3R/dG7/a+0hZEDoeBxQ54ubmgTwdfALzLQkS2R28Q8MrmQ9ALIoapG/FQn3CpSyIHxMBiJ5pnC3EcCxHZmn/vyMfhs9XwUbngzelDpS6HHBQDi51oXo9lV24lBEGUuBoioianK7RYuO0EAGD2A90R5M1xK2QZDCx2Ij7CF55KBaq0DcgpqZG6HCIiiKKImVsOoV4v4NZOAfhrv/ZSl0QOjIHFTrgq5EiI5mwhIrId6/cXYU9eFdxdFXhnVC+ut0IWxcBiRwbHNAUWjmMhIqmVVOvwzg/HAAAv3dcFHfw9JK6IHB0Dix1pHseyL78KDXpB4mqIyFmJoojZW4+gpl6PuAhfJA+JkrokcgIMLHaka7A3AryUuNhoQEbheanLISIn9cPhYvxyrBSuChneG90bCjkfBZHlMbDYEZlMZlymn+NYiEgK57UNxl2Yp97ZCV1DvCWuiJwFA4udGcL1WIhIQv/44RgqahvQOcgL0+7iLsxkPQwsdmbIpXEsWWeqodE1SlwNETmT306UY/PBM5DJgAWje8PNhRsbkvUwsNiZcF93RAV4wiCI2JtXJXU5ROQktPV6vLrlMABg8uBI9OvYTuKKyNkwsNihIZ24HgsRWdc/fzqOsxcuItzXHS/f11XqcsgJMbDYoeZxLAwsRGQN6QXnsWr3aQDA/L/0gqebi7QFkVNiYLFDg2L8IZMBJ8tqUarRSV0OETmwer0Br2w+BFEERvdtj9u7BEpdEjkpBhY75OuhRK9wNQDeZSEiy1r2ay5OldUiwEuJ1x/oJnU55MQYWOzU5fVYKiWuhIgcVU6JBh/9egoA8OaInvD1UEpcETkzBhY79edxLKIoSlwNETkagyDila8OQS+IuLd7MIb3CpG6JHJyDCx2qn9kOyhd5CjR6JBbrpW6HCJyMF/szEfWmWp4q1zwj4d6cidmkhwDi51SuSowILJpHQSOYyEicyqsrMP7Px8HALw6vBuCfVQSV0TEwGLXuK8QEZmbKIqY9fUh6BoF3BLth7EDIqQuiQgAA4tdax7HsjuvEnqDIHE1ROQINh04g52nKuHmIseCv/TmoyCyGQwsdqxnuBo+KhfU6PQ4fLZa6nKIyM6VaXT4xw/ZAICUe7sgMsBT4oqILmNgsWMKuQyDueotEZnJnG+OQqPTo1e4Gk/cGiV1OUQtMLDYuebdm7keCxHdjB8PFyP1aAlc5DK8O7o3XBT89UC2hV+Rdu7WSwNv0wvO42KDQeJqiMgeVdc1Ys63RwEAU+6IQfcwH4krIroSA4udi/T3QJhahQaDgP2nq6Quh4js0Nv/zUZ5TT2iAz0x/e5OUpdD1CoGFjsnk8k4vZmI2mzHyQpsPHAGMhnw3ujeULkqpC6JqFUMLA7g1uZxLLkMLER04+oa9Jj19SEAwIRbOqJ/pJ/EFRFdHQOLA2ieKXT0nAZV2gaJqyEie7Hw5xMoqrqIMLUKM4bGSl0O0TUxsDiAQG83dA32higCu3M5W4iIri+z6AI+35kPAHh7VC94ublIXBHRtTGwOIjmcSw7OI6FiK6jQS/gla8OQRCBh+LDcFdskNQlEV0XA4uDuLWzPwBgF8exENF1fLw9F8dLa+DnqcScB3tIXQ7RDWFgcRADo/zhIpehoLIORVV1UpdDRDbqZGkNlv56EgAw98Hu8PNUSlwR0Y1hYHEQXm4uiI/wBcDpzUTUOoMgYsbmQ2g0iLgnNggj4sKkLonohjGwOBCOYyGia1n6v1PIKLwAL0MD/jGqJ3diJrvCwOJAmtdj2Z1bCUEQJa6GiGzJtuxS/OuXEwCAOSW7EKp2l7giItMwsDiQ+AhfeCoVqNQ2IKekRupyiMhGnCqrxYsbMgE0LRD3yPC+0hZE1AYMLA7EVSHHwKimlSo5joWIAECja8TTqw+gtl6PgZF+eP2B7sC4cVKXRWQyBhYHw3EsRNRMEES8uD4TeeVahKpVWDa+L5Qu/LFP9olfuQ6meRzLvvwqNOgFiashIikt+uUE0nLKoHSRY/mEfgj0dpO6JKI2Y2BxMF2DvRHgpcTFRgMyCs9LXQ4RSST1SDE+/N8pAMD8Ub3Qu72vtAUR3SQGFgcjk8mMmyFyHAuRczpeUoOUjVkAgOQhkRjdr73EFRHdvDYFlmXLliEyMhIqlQoJCQnYt2/fVds2NjZi3rx5iImJgUqlQlxcHFJTU6/afsGCBZDJZHjhhRfaUhoBuJXjWIicVnVd0yDbugYDBkX749Xh3aQuicgsTA4sGzZsQEpKCubOnYuDBw8iLi4OSUlJKCsra7X97NmzsXz5cixZsgTZ2dmYMmUKRo0ahYyMjCva7t+/H8uXL0fv3r1N7wkZDbk0jiXrTDVqdI0SV0NE1mIQRDy/PgMFlXUI93XH0kf7wFXBG+nkGEz+Sl64cCGeeuopJCcno3v37vjkk0/g4eGBzz//vNX2q1evxquvvorhw4cjOjoaU6dOxfDhw/HBBx+0aFdbW4vx48djxYoVaNeuXdt6QwCAcF93RAV4wiCI2JtXJXU5RGQl7/98HL+dKIfKtWmQrb8XB9mS4zApsDQ0NCA9PR2JiYmXLyCXIzExEbt37271nPr6eqhUqhbH3N3dsWPHjhbHpk2bhvvvv7/Fta+lvr4eGo2mxYsuG9KpafdmPhYicg7fHzqHj7fnAgDeHd0bPcPVEldEZF4mBZaKigoYDAYEBwe3OB4cHIySkpJWz0lKSsLChQtx8uRJCIKAbdu2YcuWLSguLja2Wb9+PQ4ePIj58+ffcC3z58+HWq02viIiIkzpisMb8ueBt+vWSVwNEVlS9jkN/r7pEADgmdujMTI+XOKKiMzP4g83Fy9ejM6dOyM2NhZKpRLTp09HcnIy5PKmT11UVIS//e1v+PLLL6+4E3Mts2bNQnV1tfFVVFRkqS7YpUEx/pDJgJNltSjd9I3U5RCRhZzXNuDp1QdwsdGA2zoHYMbQWKlLIrIIkwJLQEAAFAoFSktLWxwvLS1FSEhIq+cEBgZi69at0Gq1KCgoQE5ODry8vBAdHQ0ASE9PR1lZGfr27QsXFxe4uLjgt99+w4cffggXFxcYDIZWr+vm5gYfH58WL7rM10OJXpduCe/y5P+2iByR3iDguXUZOHP+Ijr4eWDJuD5QyLkDMzkmkwKLUqlEv379kJaWZjwmCALS0tIwaNCga56rUqkQHh4OvV6PzZs3Y+TIkQCAe+65B4cPH0ZmZqbx1b9/f4wfPx6ZmZlQKBRt6BYBl5fp/8Wro8SVEJElvJuagx2nKuChVODTif3g66GUuiQii3Ex9YSUlBRMmjQJ/fv3x8CBA7Fo0SJotVokJycDACZOnIjw8HDjeJS9e/fi7NmziI+Px9mzZ/HGG29AEATMmDEDAODt7Y2ePXu2+Byenp7w9/e/4jiZ5v5eofh4ey5+8olEqUaHYJ8bf+RGRLZta8ZZrPgjHwDw/l/jEBvCu8zk2EwOLGPGjEF5eTnmzJmDkpISxMfHIzU11TgQt7Cw0Dg+BQB0Oh1mz56NvLw8eHl5Yfjw4Vi9ejV8fX3N1glqXc9wNfp3bIcDBefx5d5CpNzbReqSiMgMjpytxiubmwbZTrsrBsN7hUpcEZHlyURRFKUuwhw0Gg3UajWqq6s5nuVPvj90DtPXZiDAS4mdM++GmwsfsRHZs4raeoxYsgPnqnW4q2sgPps0gONWyK7d6O9vLoHo4JJ6hCCksRYVtQ344VDx9U8gIpvVaBAw7cuDOFetQ1SAJxaN5SBbch4MLA7OVSHHhIimJ39f7DwNB7mhRuSU3v7hGPbmV8HLzQUrJvaD2t1V6pKIrIaBxQmMffIBKF3kOHy2GgcLL0hdDhG1waYDRVi56zQAYOEjcegU5C1tQURWxsDiBPy93DAyLgwAsOrSDzwish+ZRRfw2tYjAIAXEjvjvh6tr3tF5MgYWJzEpMGRAID/Hi5GqUYnbTFEdMPKanSYsjodDXoB93YPxvN3d5a6JCJJMLA4iZ7hagyIbAe9IOLLPQVSl0NEN6BBL+DZNQdRotEhJtATCx+Jg5yDbMlJMbA4kcmDowAAX+4tRL2+9S0PiMh2zPv+KA4UnIe3mwtWTOwPbxUH2ZLzYmBxIvf1CEaoWoVKbQO+z+IUZyJbtm5fIdbsKYRMBiweF4/oQC+pSyKSFAOLE3FVyPHYLU37Cq3cxSnORLYqvaAKc75pGmT78n1dcXdssMQVEUmPgcXJjBvY4U9TnM9LXQ4R/T+lGh2mrNiJRoOI4b1C8OydMVKXRGQTGFicjJ+nEg/FN01x/mLnaWmLIaIWdI0GPLM6HeV6OboGe+OfD8dBJuMgWyKAgcUpNU9xTj1SgpJqTnEmsgUGQcQL6zORWXQBaoMOn07sB083k/enJXJYDCxOqEeYGgOj/JqmOO/lFGciqYmiiLe+z0bq0RIoFXIsL/oZHf09pS6LyKYwsDipyZfusqzdWwhdI6c4E0np09/zLi+7PyYOt9RxFh/R/8fA4qTu6/6nKc7cxZlIMt9knsX8H3MAALPv74YHeocB48ZJXBWR7WFgcVIuCjkmDGqa4vzFznxOcSaSwK5TFXh5UxYA4Ilbo/DkbdFNH2BgIboCA4sTGzugA9xc5Dh6ToP0Ak5xJrKmY8UaPLM6HY0GEff3DsVrw7tJXRKRTWNgcWJNU5zDAQBfcBdnIqs5d+Eikr/Yj5p6PQZG+eGDv3KPIKLrYWBxcn+e4lxcfVHaYoicQPXFRkz+Yh9KNDp0DvLCign9oXJVSF0Wkc1jYHFy3cN8kBDlB4Mg4ss9hVKXQ+TQ6vUGPP2fAzhRWotgHzesfHwg1B7c0JDoRjCwEJKHRAIA1u7jFGciSxEEES9tzMLe/Cp4ubngi8kDEe7rLnVZRHaDgYWQ2C0Y4b7uqNI24Lusc1KXQ+SQ5v94DN8fKoarQoblE/qhe5iP1CUR2RUGFoILd3EmsqjPd+RjxR/5AIB/PhyHIZ0CJK6IyP4wsBAAYOyACOMU5wOc4kxkNv89XIy3fsgGALwyNBYP9QmXuCIi+8TAQgCAdp5KjLr0g3Qld3EmMot9+VV4YUMmRBGYcEtHTLkjWuqSiOwWAwsZGac4Hy3BuQuc4kx0M06W1uDJVfvRoBdwX/dgvDGiB2QyrrVC1FYMLGTULdQHt0Q3TXFes4e7OBO1ValGh8lf7IdGp0ffDr74cFwfKLgwHNFNYWChFiYPjgIArOMUZ6I2qdE1YvIX+3H2wkVEB3ji35MGcGE4IjNgYKEWErsFIdzXHefrGvEtpzgTmaRBL2DqmoM4VqxBgJcbVj0+EO08lVKXReQQGFioBReFHBMv7eK8cienOBPdKFEUMXPzIew4VQEPpQJfTB6ACD8PqcsichgMLHSFMQMioHKVI7tYg/2nOcWZ6Eb886fj2JJxFgq5DMvG90Wv9mqpSyJyKAwsdAVfjz9Ncd6VL3E1RLZv9Z4CfLQ9FwAw/y+9cFfXIIkrInI8DCzUquYpzj8dLcVZTnEmuqqfj5Zg7jdHAAAvJnbBI/0jJK6IyDExsFCrYkN8MCjan1Ocia4hveA8nluXAUEExg2MwPP3dJK6JCKHxcBCVzX50i7OnOJMdKW88lo8uWo/6vUC7o4Nwlsje3JhOCILYmChq2rexflCXSO+zeQUZ6Jm5f9Zj0lf7MP5ukb0bq/G0kf7wEXBH6dElsTvMLoqhVyGSYObpjh/wV2ciQAA2no9Hk/XoajqIjr4eeDzyQPgoXSRuiwih8fAQtc0pn8HuLsqcKxYg335VVKXQySper0BU788iMPugfDzVGLV4wMR4OUmdVlEToGBha5J7eGKUX2bpziflrYYIgk1GgQ8tzYDv58oh7vQiH9P6o+oAE+pyyJyGgwsdF2TBkUCAH46WsIpzuSUDIKIlzZm4efsUihd5FhR9BP6dGgndVlEToWBha6ra4g3Bsf4QxCB1bs5xZmciyA0Lbn/bdY5uMhl+OSxvrhVe1bqsoicDgML3ZDJlxaSW7+/EBcbOMWZnIMoinjju6PYlH4Gchnw4bg+uDs2GBg3TurSiJwOAwvdkHu6BaN9u6Ypzt9k8n+X5PhEUcT8H3Pwn90FkMmADx6Jw/BeoU0fZGAhsjoGFrohCrnMOJZlJac4kxNY9MtJfPp7HgDg7Yd6YVSf9hJXROTcGFjohj3SPwLurgrklNRgL6c4kwP75LdcLE47CQCY80B3PJrQQeKKiIiBhW6Y2sMVf2me4rzztLTFEFnIql2nseDHHADA35O64vFboySuiIgABhYyUfPg25+zS3DmfJ20xRCZ2Yb9hZj77VEAwHN3d8K0u7iZIZGtYGAhk3QO9satnQKapjhzF2dyIN9knsXMLYcBAE/eGoWUe7tIXBER/RkDC5lsUvMU531FnOJMDiH1SAlSNmZBFIHxCR3w2v3duPMykY1hYCGT3R0bhAg/d1RfbMTWT7+Wuhyim/Lr8TI8t+4gDIKI0X3b462RPRlWiGwQAwuZrMUU51N1nOJMdmvXqQpMWZ2ORoOI+3uH4t3RvSCXM6wQ2SIGFmqTv16a4nxc5Y9fj5dJXQ6RyQ6crsKT/zmAer2AxG7BWDQmHi4K/kgkslX87qQ2Ubu7YvyltSlmfHUY5TX1EldEdOMOnbmA5C/2o67BgNs6B2Dpo33gyrBCZNP4HUpt9nJSV3TVVaKith4pGzMhCHw0RLYvp0SDiZ/vQ029HgOj/PDphP5QuSqkLouIroOBhdpM5arA0jO/QOUqxx8nK/DJ77lSl0R0TbnltXjss724UNeI+AhffD55ANyVDCtE9qBNgWXZsmWIjIyESqVCQkIC9u3bd9W2jY2NmDdvHmJiYqBSqRAXF4fU1NQWbebPn48BAwbA29sbQUFBeOihh3D8+PG2lEZW1nn0MLw5ogcA4IOfTyC9gEv2k20qrKzD+BV7UVHbgO6hPliVPBBebi5Sl0VEN8jkwLJhwwakpKRg7ty5OHjwIOLi4pCUlISystYHXs6ePRvLly/HkiVLkJ2djSlTpmDUqFHIyMgwtvntt98wbdo07NmzB9u2bUNjYyPuu+8+aLXatveMrGPcODzSPwIj4sJgEEQ8vy4TF+oapK6KqIVzFy7i0c/2oESjQ+cgL6x+YiDUHq5Sl0VEJpCJJs5JTUhIwIABA7B06VIAgCAIiIiIwHPPPYeZM2de0T4sLAyvvfYapk2bZjw2evRouLu7Y82aNa1+jvLycgQFBeG3337D7bfffkN1aTQaqNVqVFdXw8fHx5QukRnU6BrxwJIdKKisw33dg7F8Qj+uZUE2oaxGh7HL9yCvQotIfw9sfGYQgnxUUpdFRJfc6O9vk+6wNDQ0ID09HYmJiZcvIJcjMTERu3fvbvWc+vp6qFQtfzi4u7tjx44dV/081dXVAAA/P7+rtqmvr4dGo2nxIul4q1yxdFxfuCpk+Dm7lMv2k02o0jZgwmf7kFehRbivO7586haGFSI7ZVJgqaiogMFgQHBwcIvjwcHBKCkpafWcpKQkLFy4ECdPnoQgCNi2bRu2bNmC4uLiVtsLgoAXXngBQ4YMQc+ePa9ay/z586FWq42viIgIU7pCFtCrvRqzhnUDAPzj+2M4eq5a4orImVVfbMTEz/fieGkNgrzdsPapBIT7uktdFhG1kcVnCS1evBidO3dGbGwslEolpk+fjuTkZMjlrX/qadOm4ciRI1i/fv01rztr1ixUV1cbX0VFRZYon0yUPCQSid2C0GAQ8NzaDGjr9VKXRE6otl6PyV/sw5GzGvh7KrH2qQR09PeUuiwiugkmBZaAgAAoFAqUlpa2OF5aWoqQkJBWzwkMDMTWrVuh1WpRUFCAnJwceHl5ITo6+oq206dPx/fff49ff/0V7du3v2Ytbm5u8PHxafEi6clkMvzz4TiEqlXIq9Di9a1HpC6JnIyu0YAn3/sBGYUXoHZ3xeonEtApyFvqsojoJpkUWJRKJfr164e0tDTjMUEQkJaWhkGDBl3zXJVKhfDwcOj1emzevBkjR440fkwURUyfPh1ff/01/ve//yEqKsrEbpAtaeepxOKxfSCXAVsyzuKr9DNSl0ROoqRah7Gf7sEerQu83Fzwn8cHonsY/zND5AhMfiSUkpKCFStWYNWqVTh27BimTp0KrVaL5ORkAMDEiRMxa9YsY/u9e/diy5YtyMvLwx9//IGhQ4dCEATMmDHD2GbatGlYs2YN1q5dC29vb5SUlKCkpAQXL140QxdJCgOj/PBiYhcAwOtbj+BUWa3EFZGjSy+owoNLdyCz6ALUBh2+SB6AuAhfqcsiIjMxObCMGTMG77//PubMmYP4+HhkZmYiNTXVOBC3sLCwxYBanU6H2bNno3v37hg1ahTCw8OxY8cO+Pr6Gtt8/PHHqK6uxp133onQ0FDja8OGDTffQ5LMs3d1wuAYf1xsNGD62oPQNRqkLokc1Lp9hRj76R6U19Sja7A3vs37GgMirz7LkIjsj8nrsNgqrsNim8o0Ogxb/AcqtQ2YcEtHvPXQ1Wd+EZmqQS9g3vdHsWZPIQBgWM8QvP/XOHhu2QSMGydxdUR0IyyyDguRqYJ8VFg4Jh4AsHpPAX483Pp0diJTldfUY/xne7BmTyFkMuDl+7rgo/F94enmwrBC5IAYWMji7ugSiCl3xAAAZmw+hKKqOokrInt36MwFjFi6A/tPn4e3mws+m9gf0+/uzNWViRwYAwtZxUv3dUGfDr6o0enx/PoMNBoEqUsiO7U5/Qwe/mQ3iqt1iA70xNbpQ3BPt+Drn0hEdo2BhazCVSHHh2P7wEflgozCC3j/Z+7GTabRGwTM+y4bL23KQoNewD2xQdg6bQhiAr2kLo2IrICBhawmws8D7z3cGwCw/Lc8bD/e+g7fRP9flbYBEz/fh8935gMAnr+7E1ZM7A8fFXdcJnIWDCxkVUN7hmLCLR0BAC9tzEKZRidxRWTrss9pMGLpDuzKrYSHUoFPHuuLlPu6Qi7neBUiZ8LAQlb32v3d0C3UB5XaBrywIRMGwSFm1pMFfJd1Dn/5eCfOnL+IDn4e+PrZIRjaM1TqsohIAgwsZHUqVwWWPtoHHkoFduVW4qNfT0ldEtkYgyDi3dQcPLcuA7pGAbd1DsC304egawj3BCJyVgwsJImYQC+8NbJpEbl//XIC+/KrJK6IbEX1xUY8sWo/Pt6eCwB45vZorEweCF8PpcSVEZGUGFhIMqP7tcdf+oZDEIG/rc/AeW2D1CWRxE6W1uChZTux/Xg5VK5yLB4bj1nDu0HB8SpETo+BhST11sieiA7wRHG1Di9vyoKD7BRBbfDT0RI8tGwn8iu0CPd1x1dTBmNkfLjUZRGRjWBgIUl5urlgyaN9oHSRIy2nDJ/vPC11SWRlgiDiX9tO4JnV6dA2GHBLtB++nT4EPcPVUpdGRDaEgYUk1yNMjdn3dwMALPjxGA6fqZa4IrKWGl0jnlmTjsVpJwEAkwdHYvUTCfD3cpO4MiKyNQwsZBMm3NIRST2C0WgQMX3dQdToGqUuiSwsv0KLUR/twrbsUihd5Pjnw73xxogecFXwxxIRXYk/GcgmyGQyvDc6DuG+7iiorMNrXx/heBYH9uvxMoxYugOnymoR7OOGjc8Mwl/7R0hdFhHZMAYWshlqD1d8OK4PFHIZvs06h40HiqQuicwsveA8nly1H8lf7EeNTo9+Hdvhu+duRXyEr9SlEZGNY2Ahm9KvYzu8fF9XAMDcb4/ixBcbJK6IbpYoith+vAyPLN+N0R/vwi/HyiCTAY9VHcXapxIQ5K2SukQisgMuUhdA9P89c3s0duVW4I+TFZh6qBErymsRzR157Y5BEPHfw8X4eHsusos1AABXhQx/6dMeT98RjZhfNICLQuIqicheyEQHGSig0WigVqtRXV0NHx8fqcuhm1ReU4/7P/wDZTX1ULnKMSMpFpMHR3LDOzugazRgy8GzWP57Lgoq6wAAHkoFHh3YAU/cFoVQtbvEFRKRLbnR398MLGSzzl64iFdm/wc7vNoDAG6J9sM/H45DhJ+HxJVRa2p0jVi7txCf7chHeU09AKCdhysmD47CpMEdubQ+EbWKgYUcgjhiBNa8ugTv/HAMFxsN8FQq8Nr93TFuYARkMt5tsQUVtfVYufM0/rP7NDQ6PQAgVK3CU7dFY+zACHgo+eSZiK7uRn9/8ycJ2TTZuHGYcEtH3N45AH/fdAj7Tlfh1a8PI/VoCd4d3YuPFyRUVFWHFX/kYcP+ItTrBQBATKAnptwRg5Hx4VC6cEw/EZkP77CQ3TAIIr7YmY/3fjqOBr0Ab5UL3hzRA6P6hPNuixUdL6nBJ7/l4tusczAITT8+4iJ88eydMbi3WzDHGRGRSfhIiBzWqbIavLQxC1mXlvC/t3sw3hnVC4HeXM7dktILqvDx9lz8cqzMeOy2zgGYemcMBkX7MzQSUZswsJBD0xsELP89D4t+OYFGg4h2Hq74x0O9cH/vUKlLcyiiKGL7iXJ8vD0X+/KrAAAyGTCsZwim3tEJvdpzg0IiujkMLOQUss9p8NKmLBy7tM7Hg3FhmDeiB9p5ckbKzdAbBPx3+WZ83Bhi/LdtXkPlmTuiuS4OEZkNAws5jQa9gCX/O4mPtufCIIgI9HbDgr/0wj3dgqUuza6Iooij5zT4Luscvss6h3PVOgBNa6iMT+iAJ26NRoiaq9ISkXkxsJDTySq6gJc2ZeFUWS0A4OF+7THnwe7wUblKXJltO1la0xRSDhUjv0JrPN5OfxHJw+IwcRDXUCEiy2FgIaekazRg4bYTWPFHHkQRCFOr8N7Dcbi1c4DUpdmUgkotvj9UjO+yziGnpMZ4XOUqxz2xwXgwLgx3vjYVqm++lrBKInIGXIeFnJLKVYFXh3fDvd2D8fKmLBRU1uGxf+/FY7d0wKxh3eDp5rxf8sXVF/HDpZDSPMMKaBqbckeXQDwYF4Z7ugXDq/nfaOwjElVKRHQl3mEhh1XXoMeCH3Pwn90FAIAOfh54/69xGBjlJ3Fl1lNRW48fj5Tgu8xz2He6ynhcLgOGdArAg73DkNQjBGoPPjYjImnwkRDRJTtPVeDvm7JwrloHmQx4fEgU/p7UFSpXx9wpuPpiI346WoLvss5hV26lcXE3ABgY6YcH40IxrFcoAry4bg0RSY+BhehPNLpG/OP7bGw8cAYAEB3oiYWPxCM+wlfawsxEW6/HL8dK8V1WMX4/UY4Gg2D8WFx7NR6MC8PwXqEI8+VWBkRkWxhYiFrxv5xSvLL5MMpr6iGXAff6NKL7gB6IDfVGbIg3Itp52M3S8rpGA7YfL8d3h84h7VgpdI2XQ0psiDcejAvDA71D0dHfU8IqiYiujYGF6Cou1DVg7rdH8U3muSs+5qFUoEtwU3iJDfFG1xAfxIZ4S7oQXW29HgWVWhRW1qGgqg4FlXUorNLiUF4FaoTL4SrS3wMj4sLwQFwYugR7S1YvEZEpGFiIriO9oArpbyxEzrCHkVNcg1NltS0epfxZsI+bMbw0BRlvdArygpvLzY+DEUURldoGFFTWoaBSeymQNP25sKoOFbUNVz03TK3CA3FhGBEXhh5hPtzPh4jsDqc1E11Hv45+6HdfD+CReABNy9GfrtTiWHENjpfUIKekBjklGpw5fxGlmnqUasrx+4ly4/kKuQzRAZ7oGuKNbqE+6BrsjdhQb4T7ul8RHAyCiHMXLl4KInUoqNKioKLpjklhpRbaBsM1a/XzVKKDnwci/T3Qwd8THf080Onduei1/jO7eYRFRHQzeIeF6DpqdI04UdoUYI6X1CCnuCnIaHT6Vtt7u7mgS4g3ogI8UXEiH4Vuvig6X4dGw9W/1WQyIEztjg5+Hujo74EO/h6I9Pc0/t27tdV6160Dxo0zVzeJiCTBR0JEFiSKIoqrdS3uxBwvqUFuee1Vg4lSIUd7P3d09PNAR39PdPS/FE78PBHh526Wx0tERPaGj4SILEgmkyHM1x1hvu64KzbIeLxBLyCvohbHS2pQUFmHwFUr0HHODHQM8ESIjwoKPr4hImoTBhYiM1K6yBEb4oPYkEv/SyjrA3TiPkZERDdLLnUBRA6NY0yIiMyCgYWIiIhsHgMLERER2TwGFiIiIrJ5DCxERERk8xhYiIiIyOYxsBAREZHNY2AhIiIim8fAQkRERDaPgYWIiIhsHgMLERER2TwGFiIiIrJ5DCxERERk89oUWJYtW4bIyEioVCokJCRg3759V23b2NiIefPmISYmBiqVCnFxcUhNTb2paxIREZFzMTmwbNiwASkpKZg7dy4OHjyIuLg4JCUloaysrNX2s2fPxvLly7FkyRJkZ2djypQpGDVqFDIyMtp8TSIiInIuMlEURVNOSEhIwIABA7B06VIAgCAIiIiIwHPPPYeZM2de0T4sLAyvvfYapk2bZjw2evRouLu7Y82aNW26Zms0Gg3UajWqq6vh4+NjSpeIiIhIIjf6+9ukOywNDQ1IT09HYmLi5QvI5UhMTMTu3btbPae+vh4qlarFMXd3d+zYsaPN12y+rkajafEiIiIix2RSYKmoqIDBYEBwcHCL48HBwSgpKWn1nKSkJCxcuBAnT56EIAjYtm0btmzZguLi4jZfEwDmz58PtVptfEVERJjSFSIiIrIjFp8ltHjxYnTu3BmxsbFQKpWYPn06kpOTIZff3KeeNWsWqqurja+ioiIzVUxERES2xqTUEBAQAIVCgdLS0hbHS0tLERIS0uo5gYGB2Lp1K7RaLQoKCpCTkwMvLy9ER0e3+ZoA4ObmBh8fnxYvIiIickwmBRalUol+/fohLS3NeEwQBKSlpWHQoEHXPFelUiE8PBx6vR6bN2/GyJEjb/qaRERE5BxcTD0hJSUFkyZNQv/+/TFw4EAsWrQIWq0WycnJAICJEyciPDwc8+fPBwDs3bsXZ8+eRXx8PM6ePYs33ngDgiBgxowZN3xNIiIicm4mB5YxY8agvLwcc+bMQUlJCeLj45GammocNFtYWNhifIpOp8Ps2bORl5cHLy8vDB8+HKtXr4avr+8NX5OIiIicm8nrsNgqrsNCRERkfyyyDgsRERGRFBhYiIiIyOYxsBAREZHNM3nQra1qHorDJfqJiIjsR/Pv7esNqXWYwFJTUwMAXKKfiIjIDtXU1ECtVl/14w4zS0gQBJw7dw7e3t6QyWRSl2MxGo0GERERKCoqcvjZUOyr43Km/rKvjsuZ+mvJvoqiiJqaGoSFhV1z2x6HucMil8vRvn17qcuwGmfajoB9dVzO1F/21XE5U38t1ddr3VlpxkG3REREZPMYWIiIiMjmMbDYGTc3N8ydOxdubm5Sl2Jx7Kvjcqb+sq+Oy5n6awt9dZhBt0REROS4eIeFiIiIbB4DCxEREdk8BhYiIiKyeQwsREREZPMYWCS2bNkyREZGQqVSISEhAfv27btm+02bNiE2NhYqlQq9evXCf//73xYfnzx5MmQyWYvX0KFDLdkFk5jS36NHj2L06NGIjIyETCbDokWLbvqa1mTuvr7xxhtXvLexsbEW7MGNM6WvK1aswG233YZ27dqhXbt2SExMvKK9KIqYM2cOQkND4e7ujsTERJw8edLS3bhh5u6vLX/fmtLXLVu2oH///vD19YWnpyfi4+OxevXqFm1s+b01d19t+X0F2v6zc/369ZDJZHjooYdaHLf4eyuSZNavXy8qlUrx888/F48ePSo+9dRToq+vr1haWtpq+507d4oKhUJ87733xOzsbHH27Nmiq6urePjwYWObSZMmiUOHDhWLi4uNr6qqKmt16ZpM7e++ffvEl19+WVy3bp0YEhIi/utf/7rpa1qLJfo6d+5csUePHi3e2/Lycgv35PpM7eujjz4qLlu2TMzIyBCPHTsmTp48WVSr1eKZM2eMbRYsWCCq1Wpx69atYlZWljhixAgxKipKvHjxorW6dVWW6K+tft+a2tdff/1V3LJli5idnS2eOnVKXLRokahQKMTU1FRjG1t9by3RV1t9X0Wx7T878/PzxfDwcPG2224TR44c2eJjln5vGVgkNHDgQHHatGnGvxsMBjEsLEycP39+q+0feeQR8f77729xLCEhQXzmmWeMf580adIVX0S2wtT+/lnHjh1b/SV+M9e0JEv0de7cuWJcXJwZqzSPm30P9Hq96O3tLa5atUoURVEUBEEMCQkR//nPfxrbXLhwQXRzcxPXrVtn3uLbwNz9FUXb/b41x/dXnz59xNmzZ4uiaNvvrbn7Koq2+76KYtv6q9frxcGDB4ufffbZFX2zxnvLR0ISaWhoQHp6OhITE43H5HI5EhMTsXv37lbP2b17d4v2AJCUlHRF++3btyMoKAhdu3bF1KlTUVlZaf4OmKgt/ZXimuZgybpOnjyJsLAwREdHY/z48SgsLLzZcm+KOfpaV1eHxsZG+Pn5AQDy8/NRUlLS4ppqtRoJCQmSvq+AZfrbzNa+b2+2r6IoIi0tDcePH8ftt98OwHbfW0v0tZmtva9A2/s7b948BAUF4YknnrjiY9Z4bx1m80N7U1FRAYPBgODg4BbHg4ODkZOT0+o5JSUlrbYvKSkx/n3o0KH4y1/+gqioKOTm5uLVV1/FsGHDsHv3bigUCvN35Aa1pb9SXNMcLFVXQkICVq5cia5du6K4uBhvvvkmbrvtNhw5cgTe3t43W3abmKOvr7zyCsLCwow/6Jq/nq/3tS4FS/QXsM3v27b2tbq6GuHh4aivr4dCocBHH32Ee++9F4DtvreW6Ctgm+8r0Lb+7tixA//+97+RmZnZ6set8d4ysDiYsWPHGv/cq1cv9O7dGzExMdi+fTvuueceCSujmzVs2DDjn3v37o2EhAR07NgRGzdubPV/PPZgwYIFWL9+PbZv3w6VSiV1ORZ3tf460vett7c3MjMzUVtbi7S0NKSkpCA6Ohp33nmn1KWZ3fX66ijva01NDSZMmIAVK1YgICBAsjoYWCQSEBAAhUKB0tLSFsdLS0sREhLS6jkhISEmtQeA6OhoBAQE4NSpU5J+g7Slv1Jc0xysVZevry+6dOmCU6dOme2aprqZvr7//vtYsGABfvnlF/Tu3dt4vPm80tJShIaGtrhmfHy8+YpvA0v0tzW28H3b1r7K5XJ06tQJABAfH49jx45h/vz5uPPOO232vbVEX1tjC+8rYHp/c3Nzcfr0aTz44IPGY4IgAABcXFxw/Phxq7y3HMMiEaVSiX79+iEtLc14TBAEpKWlYdCgQa2eM2jQoBbtAWDbtm1XbQ8AZ86cQWVlZYsvICm0pb9SXNMcrFVXbW0tcnNzJX1v29rX9957D2+99RZSU1PRv3//Fh+LiopCSEhIi2tqNBrs3btX0vcVsEx/W2ML37fm+joWBAH19fUAbPe9tURfW2ML7ytgen9jY2Nx+PBhZGZmGl8jRozAXXfdhczMTERERFjnvTXL0F1qk/Xr14tubm7iypUrxezsbPHpp58WfX19xZKSElEURXHChAnizJkzje137twpuri4iO+//7547Ngxce7cuS2mNdfU1Igvv/yyuHv3bjE/P1/85ZdfxL59+4qdO3cWdTqdJH38M1P7W19fL2ZkZIgZGRliaGio+PLLL4sZGRniyZMnb/iaUrFEX1966SVx+/btYn5+vrhz504xMTFRDAgIEMvKyqzevz8zta8LFiwQlUql+NVXX7WY7llTU9Oija+vr/jNN9+Ihw4dEkeOHGkTU19F0fz9teXvW1P7+s4774g///yzmJubK2ZnZ4vvv/++6OLiIq5YscLYxlbfW3P31ZbfV1E0vb//X2szoCz93jKwSGzJkiVihw4dRKVSKQ4cOFDcs2eP8WN33HGHOGnSpBbtN27cKHbp0kVUKpVijx49xB9++MH4sbq6OvG+++4TAwMDRVdXV7Fjx47iU089Jfkv7z8zpb/5+fkigCted9xxxw1fU0rm7uuYMWPE0NBQUalUiuHh4eKYMWPEU6dOWbFHV2dKXzt27NhqX+fOnWtsIwiC+Prrr4vBwcGim5ubeM8994jHjx+3Yo+uzZz9tfXvW1P6+tprr4mdOnUSVSqV2K5dO3HQoEHi+vXrW1zPlt9bc/bV1t9XUTT998+ftRZYLP3eykRRFM1zr4aIiIjIMjiGhYiIiGweAwsRERHZPAYWIiIisnkMLERERGTzGFiIiIjI5jGwEBERkc1jYCEiIiKbx8BCRERENo+BhYiIiGweAwsRERHZPAYWIiIisnkMLERERGTz/g8wdQiGwPbthAAAAABJRU5ErkJggg=="
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 5
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "As we can see, it is good to choose $\\lambda_1$ around 0.15."
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:09:57.714109Z",
     "start_time": "2024-08-26T16:09:57.636907Z"
    }
   },
   "cell_type": "code",
   "source": [
    "val_avg = np.mean(val_err_l1, axis=0)\n",
    "l1_est = l1_list[np.argmin(val_avg)]\n",
    "print(l1_est)"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.13999999999999999\n"
     ]
    }
   ],
   "execution_count": 6
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "We use the one-standard error rule, which gives 0.16. "
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:09:57.824723Z",
     "start_time": "2024-08-26T16:09:57.747560Z"
    }
   },
   "cell_type": "code",
   "source": "parameter_tuning.get_lambda_one_se_1d(val_err_l1.squeeze(), l1_list)",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.16"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "execution_count": 7
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "Also note that, if the data is large (especially node size), we may wish to divide the search to two parts.\n",
    "First, we can search a wider range of $\\lambda_1$ with large step size.\n",
    "After we roughly know the range of good $\\lambda_1$ value, we can use a smaller range and finer step size."
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Step 2: search $\\lambda_2$ using cross validation\n",
    "\n",
    "We now search $\\lambda_2$ with $\\lambda_1$ fixed to be 0.15. The code is almost the same as the previous one.\n",
    "Usally, we do not need to search for too large $\\lambda_2$, and here an upper bound of 0.16 is sufficient."
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:09:57.918591Z",
     "start_time": "2024-08-26T16:09:57.857635Z"
    }
   },
   "cell_type": "code",
   "source": "l2_list = np.arange(0.01, 0.16, 0.01)",
   "outputs": [],
   "execution_count": 8
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:10:11.275054Z",
     "start_time": "2024-08-26T16:09:57.950710Z"
    }
   },
   "cell_type": "code",
   "source": [
    "val_err_l2 = parameter_tuning.cv_2d(\n",
    "    dat1,\n",
    "    dat2,\n",
    "    n_cv=5,\n",
    "    ratio_val=0.5, \n",
    "    lambda1_lst=[0.15],\n",
    "    lambda2_lst=l2_list,\n",
    "    dep_mat=dep_mat,\n",
    "    cores=cores,\n",
    ")"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.15 0.01 s0\n",
      "0 0.15 0.02 s0\n",
      "0 0.15 0.03 s0\n",
      "0 0.15 0.04 s0\n",
      "0 0.15 0.05 s0\n",
      "0 0.15 0.060000000000000005 s0\n",
      "0 0.15 0.06999999999999999 s0\n",
      "0 0.15 0.08 s0\n",
      "0 0.15 0.09 s0\n",
      "0 0.15 0.09999999999999999 s0\n",
      "0 0.15 0.11 s0\n",
      "0 0.15 0.12 s0\n",
      "0 0.15 0.13 s0\n",
      "0 0.15 0.14 s0\n",
      "0 0.15 0.15000000000000002 s0\n",
      "1 0.15 0.01 s0\n",
      "1 0.15 0.02 s0\n",
      "1 0.15 0.03 s0\n",
      "1 0.15 0.04 s0\n",
      "1 0.15 0.05 s0\n",
      "1 0.15 0.060000000000000005 s0\n",
      "1 0.15 0.06999999999999999 s0\n",
      "1 0.15 0.08 s0\n",
      "1 0.15 0.09 s0\n",
      "1 0.15 0.09999999999999999 s0\n",
      "1 0.15 0.11 s0\n",
      "1 0.15 0.12 s0\n",
      "1 0.15 0.13 s0\n",
      "1 0.15 0.14 s0\n",
      "1 0.15 0.15000000000000002 s0\n",
      "2 0.15 0.01 s0\n",
      "2 0.15 0.02 s0\n",
      "2 0.15 0.03 s0\n",
      "2 0.15 0.04 s0\n",
      "2 0.15 0.05 s0\n",
      "2 0.15 0.060000000000000005 s0\n",
      "2 0.15 0.06999999999999999 s0\n",
      "2 0.15 0.08 s0\n",
      "2 0.15 0.09 s0\n",
      "2 0.15 0.09999999999999999 s0\n",
      "2 0.15 0.11 s0\n",
      "2 0.15 0.12 s0\n",
      "2 0.15 0.13 s0\n",
      "2 0.15 0.14 s0\n",
      "2 0.15 0.15000000000000002 s0\n",
      "3 0.15 0.01 s0\n",
      "3 0.15 0.02 s0\n",
      "3 0.15 0.03 s0\n",
      "3 0.15 0.04 s0\n",
      "3 0.15 0.05 s0\n",
      "3 0.15 0.060000000000000005 s0\n",
      "3 0.15 0.06999999999999999 s0\n",
      "3 0.15 0.08 s0\n",
      "3 0.15 0.09 s0\n",
      "3 0.15 0.09999999999999999 s0\n",
      "3 0.15 0.11 s0\n",
      "3 0.15 0.12 s0\n",
      "3 0.15 0.13 s0\n",
      "3 0.15 0.14 s0\n",
      "3 0.15 0.15000000000000002 s0\n",
      "4 0.15 0.01 s0\n",
      "4 0.15 0.02 s0\n",
      "4 0.15 0.03 s0\n",
      "4 0.15 0.04 s0\n",
      "4 0.15 0.05 s0\n",
      "4 0.15 0.060000000000000005 s0\n",
      "4 0.15 0.06999999999999999 s0\n",
      "4 0.15 0.08 s0\n",
      "4 0.15 0.09 s0\n",
      "4 0.15 0.09999999999999999 s0\n",
      "4 0.15 0.11 s0\n",
      "4 0.15 0.12 s0\n",
      "4 0.15 0.13 s0\n",
      "4 0.15 0.14 s0\n",
      "4 0.15 0.15000000000000002 s0\n"
     ]
    }
   ],
   "execution_count": 9
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Again, we can plot the error curves along $\\lambda_2$. And we can choose values from 0.04 to 0.06, considering the uncertainty."
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:10:11.546329Z",
     "start_time": "2024-08-26T16:10:11.324920Z"
    }
   },
   "cell_type": "code",
   "source": "parameter_tuning.plot_error_1d(val_err_l2.squeeze(), lambda_lst=l2_list)",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGdCAYAAADqsoKGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNIklEQVR4nO3deVxU5eIG8GdmYBh2ZBEcFnHfBQUhzbLu5YZp7qlZKVGZdtUWStMk7Xq7l26LP0256jWT0nIphUwLr5G5K7K5L6jIvivDJjDMnN8fGMUVlcEZzjDzfD+f+YPjWZ53Ynmac857JIIgCCAiIiJq56RiByAiIiLSB5YaIiIiMgksNURERGQSWGqIiIjIJLDUEBERkUlgqSEiIiKTwFJDREREJoGlhoiIiEyChdgB2opWq0VeXh7s7e0hkUjEjkNEREQtIAgCKioqoFQqIZXe+7MYsyk1eXl58Pb2FjsGERERtUJ2dja8vLzuuY7ZlBp7e3sADW+Kg4ODyGmIiIioJcrLy+Ht7d34d/xezKbU/HbKycHBgaWGiIionWnJpSOtulA4Ojoavr6+UCgUCA4ORmJi4l3XVavVWLZsGbp16waFQgE/Pz/Ex8c3WefgwYMYM2YMlEolJBIJ4uLi7thPZWUl5s6dCy8vL1hbW6Nv375Yu3Zta+ITERGRCdK51Gzbtg0RERFYunQpUlJS4Ofnh9DQUBQVFTW7fmRkJNatW4dVq1bh/PnzmD17NiZMmIDU1NTGdaqqquDn54fo6Oi7HjciIgLx8fHYvHkzLly4gDfeeANz587Frl27dB0CERERmSCJIAiCLhsEBwdjyJAhWL16NYCGu4q8vb0xb948LFy48I71lUolFi9ejDlz5jQumzRpEqytrbF58+Y7A0kkiI2Nxfjx45ss79+/P6ZOnYr33nuvcVlAQACefPJJfPDBB/fNXV5eDkdHR6hUKp5+IiIiaid0+fut0yc1dXV1SE5ORkhIyO87kEoREhKCY8eONbtNbW0tFApFk2XW1tY4fPiwLofGsGHDsGvXLuTm5kIQBOzfvx+XL1/GE088cdfjlpeXN3kRERGR6dKp1JSUlECj0cDd3b3Jcnd3dxQUFDS7TWhoKJYvX4709HRotVrs27cPO3fuRH5+vk5BV61ahb59+8LLywtyuRwjR45EdHQ0Hn300WbXj4qKgqOjY+OLt3MTERGZNoPPKLxy5Ur06NEDvXv3hlwux9y5cxEeHn7fCXT+16pVq3D8+HHs2rULycnJ+PTTTzFnzhz8/PPPza6/aNEiqFSqxld2drY+hkNERERGSqdbul1dXSGTyVBYWNhkeWFhITw8PJrdxs3NDXFxcaipqUFpaSmUSiUWLlyIrl27tvi4t27dwrvvvovY2FiMHj0aADBw4ECkpaXhk08+aXI67DdWVlawsrLSYXRERETUnun0cYlcLkdAQAASEhIal2m1WiQkJGDo0KH33FahUMDT0xP19fXYsWMHxo0b1+LjqtVqqNXqOz7dkclk0Gq1ugyBiIiITJTOk+9FREQgLCwMgYGBCAoKwooVK1BVVYXw8HAAwIwZM+Dp6YmoqCgAwIkTJ5Cbmwt/f3/k5ubi/fffh1arxYIFCxr3WVlZiStXrjR+nZGRgbS0NDg7O8PHxwcODg4YMWIE5s+fD2tra3Tu3BkHDhzAV199heXLlz/oe0BEREQmQOdSM3XqVBQXF2PJkiUoKCiAv78/4uPjGy8ezsrKavKJSk1NDSIjI3Ht2jXY2dlh1KhR2LRpE5ycnBrXSUpKwuOPP974dUREBAAgLCwMMTExAICtW7di0aJFeO6553Djxg107twZ//jHPzB79uzWjJuIiIhMjM7z1LRXnKeGiIio/THYPDVERERExoqlhoiIiEwCSw0RERGZBJYaIiIiMgksNURERGQSWGqIiIjIJLDUEBERkUlgqSEiIiKTwFJDREREJoGlhoiIiEwCSw0RERHpx5Ytoh6epYaIiIj0g6WGiIiI6MGx1BAREZFJYKkhIiIik2AhdgAiIiJq33JuVuP7tDzcdH8IkSLmYKkhIiIinaluqfHTmXzsTM1FYsYNAIClc3/Mra6Dk41clEwsNURERNQidfVaHLhcjNjUHPx8oQh19VoAgEQCDO3qggkHv4OVxWjR8rHUEBERmYotW4Bp0/S6S0EQkJpdhtiUXOw+nYeb1erGf+vlbo8Jgz0xzl+JTo7WwA//BOQyvR5fFyw1REREpkKPpSaztAqxqbmIS83F9dLqxuVu9lYY56fEhMGe6NvJARKJRC/H0weWGiIiIgIA3Kyqw+4z+YhLzUVy5s3G5daWMozs74EJgzzxcHdXyKTGU2T+iKWGiIjIjNXWa/DLhSLEpuZi/6UiqDUCAEAqAR7u7oqJgz3xRF8P2FoZf2Uw/oRERESkV1qtgKTMm4hNzcWe03kor6lv/Le+nRwwcbAnxvop0dFBIWJK3bHUEBERmYmrxZWIS81FbGoucm7ealzeyVGBcf6emDDIE7087Ft/AD1fpKwrlhoiIiITVlpZix9O5SE2NRenclSNy+2sLPDk7etkgru66Oc6GZYaIiIi0qcatQb7zhciNjUXBy4XQ6NtuE5GJpXg0R6umDDYC3/p4w5rEW+/NgSWGiIiIhOg1Qo4bqNE7Len8NPZAlTW/n6dzEAvR0wY5Ikxfkq42lmJmNKwWGqIiIjasczSKnyXnIMdyTnI8x0DJOcAADydrDFhkCfGD/JE9452IqdsGyw1RERE7Ux1XT1+OlOA7UnZOHH7uUsAYK+pxVNDu2PCIC8Edu4AqZHOJ2MoLDVERET6YoDHFPxGEASkZJXh26Rs7D6d33h6SSIBHunhhskBXvjLkjlQfBxrkOO3Byw1RERE+mKAUlNUXoOdqbn4NikbV4urGpf7ONtgSqAXJg72gtLJumGhoNHrsdsblhoiIiIjo9Zo8cvFInyblI39l36/e8naUoZRAzphcqAXgnydze700v2w1BARERmJSwUV+DYpG7GpuSitqmtcHtC5AyYHeGH0wE6wV1iKmNC4sdQQERGJSHVLjR9O5eHbpOwmk+O52Vth4mBPTA7wNpu7lx4USw0REVEb02oFHLtWiu1J2Yg/W4Daei0AwEIqwZ/7dMSUQG+M6OkGC5lU5KTtC0sNERFRG8m+UY3vknPwXXIOcst+f/ZST3c7TAn0xvhBniY9OZ6hsdQQEREZUI1ag73nGuaUOXKltHG5vcIC4/yVmBzgjYFejpBIeNHvg2KpISIi0jNBEHA6R4XtSdnYdSoPFTW/P7JgeHdXTA70Qmg/DygsTevZS2JjqSEiItKTEpkCcYeuYXtSNi4XVjYu9+pgjacDvDBpsBe8nW0MF0Dkp2SLjaWGiIjoAQiCgCNXSrHp+HUk9Hwe9XsuAACsLKR4sr8HpgR646GuLm0zpwxLDREREemqRq3BrrQ8fHEkAxcLKhoWSmTw83bClEAvPDVQCUdrzinTllhqiIiIdFBcUYvNxzPx9YlMlFQ2TJBnI5dhcoAXntv8MXp++KXICc0XSw0REVELXCwox4ZDGfg+LQ91moZ5ZZSOCoQN88UzQT4Nn8psuClySvPGUkNERHQXWq2AXy8XYcPhjCa3Y/t7O+Gl4V0wsr8HLDlBntFgqSEiIvof1XX12JGSi41HMnDt9pOxpRLgyf6d8OLwLgjo3EHkhNQclhoiIqLbClQ1+PLYdXxzIguqW2oAgL2VBZ4J8kbYMF94dTDg7dj0wFhqiIjINGzZ0upbmk/nlGHD4QzsOZ2Peq0AAPBxtkH4w76YHOgNOyv+uWwP+F+JiIhMg46lRqMVsO98ATYczsDJ679f4BvUxRkvDe+CkD7ukLXF3DKkNyw1RERkVipq1NielIOYoxnIvtHwUEkLqQRj/JR4aXgX9Pd0FDkhtRZLDRERmYXsG9WIOXod209mo6K24VlMTjaWeC7YBzOG+sLdQSFyQnpQLDVERGSyBEFAcuZNbDicgb3nCnD7chl0c7PFi8O7YOIgL1jL+VBJU8FSQ0REJket0eLHM/n44nAGTuWoGpc/0sMVLw7vghE93NrmWUzUplo1Y1B0dDR8fX2hUCgQHByMxMTEu66rVquxbNkydOvWDQqFAn5+foiPj2+yzsGDBzFmzBgolUpIJBLExcU1u68LFy5g7NixcHR0hK2tLYYMGYKsrKzWDIGIiEyQqlqNNb9exaMf7cfrW9NwKkcFuYUUUwO9sfeNR7HppWA83quj4QqNmT9QUmw6f1Kzbds2REREYO3atQgODsaKFSsQGhqKS5cuoWPHjnesHxkZic2bN2P9+vXo3bs39u7diwkTJuDo0aMYNGgQAKCqqgp+fn548cUXMXHixGaPe/XqVQwfPhwvvfQS/va3v8HBwQHnzp2DQsFzoERE5i6v7BbWeAzHd1EJuKXWAABc7eSY/pAvnnvIB652Vm0ThKVGVBJBEARdNggODsaQIUOwevVqAIBWq4W3tzfmzZuHhQsX3rG+UqnE4sWLMWfOnMZlkyZNgrW1NTZv3nxnIIkEsbGxGD9+fJPlzzzzDCwtLbFp0yZd4jYqLy+Ho6MjVCoVHBwcWrUPIiIyPnvPFWD+t6dQXtNw8W9vD3u8NLwLxvorYWXB62XaO13+fut0+qmurg7JyckICQn5fQdSKUJCQnDs2LFmt6mtrb3j0xRra2scPny4xcfVarXYs2cPevbsidDQUHTs2BHBwcF3PU3123HLy8ubvIiIyHTU1mvw/q5zmLUpGeU19fC7VYRvXg7GT68/gsmB3iw0ZkinUlNSUgKNRgN3d/cmy93d3VFQUNDsNqGhoVi+fDnS09Oh1Wqxb98+7Ny5E/n5+S0+blFRESorK/Hhhx9i5MiR+O9//4sJEyZg4sSJOHDgQLPbREVFwdHRsfHl7e3d8oESEZFRyyqtxtNrjiHm6HUAwMxHuuDbjO8xrLsrJBJeAGyuDP5o0ZUrV6JHjx7o3bs35HI55s6di/DwcEilLT+0VtvwiPdx48bhzTffhL+/PxYuXIinnnoKa9eubXabRYsWQaVSNb6ys7P1Mh4iIhLXT2fyMfqzQziTq4KTjSU+nxGIxaP7Qg6t2NFIZDqVGldXV8hkMhQWFjZZXlhYCA8Pj2a3cXNzQ1xcHKqqqpCZmYmLFy/Czs4OXbt21em4FhYW6Nu3b5Plffr0uevdT1ZWVnBwcGjyIiKi9qtGrcGS78/i1a9TUFFbj4DOHbDntUcQ0tf9/huTWdCp1MjlcgQEBCAhIaFxmVarRUJCAoYOHXrPbRUKBTw9PVFfX48dO3Zg3LhxOh13yJAhuHTpUpPlly9fRufOnXUZAhERtUPXS6owac1RfHUsEwAwa0RXbH3lIXg6WYucjIyJzrd0R0REICwsDIGBgQgKCsKKFStQVVWF8PBwAMCMGTPg6emJqKgoAMCJEyeQm5sLf39/5Obm4v3334dWq8WCBQsa91lZWYkrV640fp2RkYG0tDQ4OzvDx8cHADB//nxMnToVjz76KB5//HHEx8fjhx9+wK+//vog4yciIiP3w6k8LNp5BpW19ehgY4nlU/zxeO87pxAh0rnUTJ06FcXFxViyZAkKCgrg7++P+Pj4xouHs7KymlwvU1NTg8jISFy7dg12dnYYNWoUNm3aBCcnp8Z1kpKS8Pjjjzd+HRERAQAICwtDTEwMAGDChAlYu3YtoqKi8Nprr6FXr17YsWMHhg8f3ppxExGRkatRa7Bs93l8c6LhMoMhvh3w2bRB6OTIT2eoeTrPU9NecZ4aIqL242pxJeZ8nYKLBRUAgL8+1g0Rf+kJC9k9rpoYOxbYtauNElJb0eXvN5/9RERERiUuNRfvxp5BdZ0GLrZyLJ/qjxE93e6/IWfzNXssNUREZBRq1A2T6W092TAFR3AXZ3w2bRDcHVr4OByWGrNn8HlqiIjITGzZ0upNrxRVYtzqI9h6MhsSCfDan7rj65eDW15oiMBSQ0RE+tLKUrMjOQdjVh3GpcIKuNpZYdOLwYh4ote9r58hagZPPxERkSiq6+qx9Ptz+DY5BwAwrJsLVjzjj472/HSGWoelhoiI2tzlwgrM+ToF6UWVkEqA1//cE3P/1B0yKZ/bRK3HUkNERG1GEAR8m5yDJd+fRY1aCzd7K3z2zCAM7eYidjQyASw1RETUJqpq6/Fe3FnsTM0FADzSwxXLp/jDzd5K5GRkKlhqiIjI4C4WlGPO1ym4WlwFqQSI+EtP/PWx7pDydBPpEUsNEREZjCAI2HYyG0t3nUNtvRbuDg2nm4K78nQT6R9LDRERGURlbT0Wx57B92l5AIARPd2wfIofXOx4uokMg6WGiIj07nxeOeZ+k4JrJVWQSSV4+4lemPVoV55uIoNiqSEiIr0RBAFfn8jCst3nUVevRSdHBVZNG4RAX2exo5EZYKkhIiK9qJBaYuGWVOw5nQ8A+FPvjvh0sh862MpFTkbmgqWGiIge2NlcFeZ0nYTM0/mwkEqwYGQvvDycp5uobbHUEBHRA8koqcK09cdRIXeEp5M1Vj07CIN9Oogdi8wQnxZGREStVllbj1e+SkJFTT0GVxdgz2vDWWhINCw1RETUKoIgYP63p5BeVImO9lZYm/1fONnw+hkSD0sNERG1ypoDV/HT2QJYyiRY83wAOmpuiR2JzBxLDRER6ezA5WJ8vPcSAOD9sf0Q0JmnnEh8LDVERKSTrNJqvLYlFYIAPDPEG88G+YgdiQgASw0REemguq4er2xKguqWGv7eTvjbuH6QSHjbNhkHlhoiImoRQRCwcMcZXCyogKudHGueHwwrC5nYsYgasdQQEVGLbDicgV2n8mAhlSD62cHo5GgtdiSiJlhqiIjovo5eKcE/f7wAAHjvqb4I7uoiciKiO7HUEBHRPeXcrMbcLanQCsDEwZ6YMbSz2JGImsVSQ0REd1Wj1mD25mTcqKpDf08H/HPCAF4YTEaLpYaIiJolCALejT2Ds7nlcLaVY+3zAVBY8sJgMl4sNURE1KyvjmViZ0oupBJg9bRB8OpgI3YkontiqSEiojskZtzA33efBwC8O6oPhnV3vf9G06YZOBXRvbHUEBFRE/mqW/jr18mo1woY46fES8O7tGxDlhoSGUsNERE1qq3X4NXNKSiprENvD3v8axIvDKb2g6WGiIgavb/rHNKyy+BobYn/TA+EjdxC7EhELcZSQ0REAIBvTmRhS2I2JBLgs2mD4OPCC4OpfWGpISIiJGfexNJdZwEA80N7YURPN5ETEemOpYaIyMwVldfg1c3JUGsEPNnfA6+O6CZ2JKJWYakhIjJjdfVa/PXrFBRV1KJHRzt8PNmPFwZTu8VSQ0Rkxv6++zySMm/CXmGB/8wIhJ0VLwym9oulhojITG1Pysam45kAgBVT/dHF1VbkREQPhqWGiMgMnc4pQ2Rcw4XBb4b0xJ/7uIuciOjBsdQQEZmZkspazN6UjLp6LUL6uGPen7qLHYlIL1hqiIjMiFqjxZyvU5CnqkFXV1ssn+oHqZQXBpNpYKkhIjIVW7bcd5WoHy/iRMYN2Mpl+M+MADgoLNsgGFHbYKkhIjIV9yk1cam5+OJIBgDg0yn+6N7Rvi1SEbUZlhoiIjNwNleFhTtPAwDmPt4dI/t7iJyISP9YaoiITNzNqjrM3pyMGrUWI3q64c2/9BQ7EpFBsNQQEZmweo0W87akIufmLfg42+CzZwZBxguDyUSx1BARmbCP/3sJh6+UwNqy4cJgRxteGEymi6WGiMhE7T6dh3UHrgEAPp48EL09HERORGRYLDVERCboYkE55n/bcGHwrEe74qmBSpETERkeSw0RkYlRVasxa1Mybqk1eLi7C+aH9hI7ElGbYKkhIjIhGq2A17elIrO0Gp5O1lg1bTAsZPxVT+aB3+lERCZkxc+X8eulYlhZSLFuegCcbeViRyJqM60qNdHR0fD19YVCoUBwcDASExPvuq5arcayZcvQrVs3KBQK+Pn5IT4+vsk6Bw8exJgxY6BUKiGRSBAXF3fP48+ePRsSiQQrVqxoTXwiIpO0194Xq365AgD4cNIA9Pd0FDkRUdvSudRs27YNERERWLp0KVJSUuDn54fQ0FAUFRU1u35kZCTWrVuHVatW4fz585g9ezYmTJiA1NTUxnWqqqrg5+eH6Ojo+x4/NjYWx48fh1LJi96IiH5zpagCbykfBwCEP+yLCYO8RE5E1PYkgiAIumwQHByMIUOGYPXq1QAArVYLb29vzJs3DwsXLrxjfaVSicWLF2POnDmNyyZNmgRra2ts3rz5zkASCWJjYzF+/Pg7/i03NxfBwcHYu3cvRo8ejTfeeANvvPFGi3KXl5fD0dERKpUKDg68rZGITIeqWo0Ja47gWnEVgrs4Y/PLwbDkdTRkInT5+63Td31dXR2Sk5MREhLy+w6kUoSEhODYsWPNblNbWwuFQtFkmbW1NQ4fPqzLoaHVajF9+nTMnz8f/fr1u+/6tbW1KC8vb/IiIjI1NWoNZm5KwrXiKnRSVyL6ucEsNGS2dPrOLykpgUajgbu7e5Pl7u7uKCgoaHab0NBQLF++HOnp6dBqtdi3bx927tyJ/Px8nYL+61//goWFBV577bUWrR8VFQVHR8fGl7e3t07HIyIydlqtgLe+PYXEjBuwt7LAF1k/wdXOSuxYRKIxeJ1fuXIlevTogd69e0Mul2Pu3LkIDw+HVNryQycnJ2PlypWIiYmBRNKyZ5YsWrQIKpWq8ZWdnd3aIRARGaV//HgBe07nw1ImwbrpAehTe0PsSESi0qnUuLq6QiaTobCwsMnywsJCeHg0/xh7Nzc3xMXFoaqqCpmZmbh48SLs7OzQtWvXFh/30KFDKCoqgo+PDywsLGBhYYHMzEy89dZb8PX1bXYbKysrODg4NHkREZmKzw9dw4bDGQCATyb7YVh3V5ETEYlPp1Ijl8sREBCAhISExmVarRYJCQkYOnToPbdVKBTw9PREfX09duzYgXHjxrX4uNOnT8fp06eRlpbW+FIqlZg/fz727t2ryxCIiNq93afz8MGeCwCAhU/2xjh/T5ETERkHC103iIiIQFhYGAIDAxEUFIQVK1agqqoK4eHhAIAZM2bA09MTUVFRAIATJ04gNzcX/v7+yM3Nxfvvvw+tVosFCxY07rOyshJXrlxp/DojIwNpaWlwdnaGj48PXFxc4OLi0iSHpaUlPDw80KsXp/8mIvNx4lopIradAgCEDe2MWY+2/FNvIlOnc6mZOnUqiouLsWTJEhQUFMDf3x/x8fGNFw9nZWU1uV6mpqYGkZGRuHbtGuzs7DBq1Chs2rQJTk5OjeskJSXh8ccfb/w6IiICABAWFoaYmJhWDo2IyLRcLqzAzK+SUKfRIrSfO5aM6dfi6wyJzIHO89S0V5ynhojaswJVDSb++wjyVDUI6NwBX78cDIWlrOlKY8cCu3aJE5DIQAw2Tw0REbW98ho1XtiYiDxVDbq62eLzGYF3FhoiYqkhIjJmdfVavLo5GRcLKuBqZ4Uvw4PQgQ+pJGoWSw0RkZESBAELvjuFI1dKYSuXISZ8CLydbcSORWS0WGqIiIzUR3svIS4tDxZSCf79fACfuk10Hyw1RERGaNOx61jz61UAQNTEARjR003kRETGj6WGiMjI7D1XgCW7zgEA3vpLT0wO5LPriFqCpYaIyIgkZ97Aa1tSIQjAtCBvzP1T95ZvPG2a4YIRtQMsNURERuJqcSVe+jIJtfVa/Ll3R/x9XH/dJtdjqSEzx1JDRGQEiipq8MLGRJRVq+Hn7YRVzw6ChYy/ool0wZ8YIiKRVdXW46WYJGTfuIXOLjbYEBYIG7nOT7EhMnssNUREIlJrtPjr1yk4k6uCi60cX4YHwdXOSuxYRO0SSw0RkUgEQcC7O8/gwOViKCyl2PDCEPi62oodi6jdYqkhItKXLVt0Wn3Fz+n4NjkHUgkQ/exg+Hs7GSYXkZlgqSEi0hcdSs3WxCysTEgHAHwwfgD+3MfdUKmIzAZLDRFRG9t/sQiL484CAOb9qTueDfYRORGRaWCpISJqQ6eyy/DXr1Og0Qp4OsALEX/pKXYkIpPBUkNE1EYyS6vwYsxJ3FJr8GhPN0RNHKDb5HpEdE8sNUREbaC0shZhXySitKoO/ZQO+Pdzg2HJyfWI9Io/UUREBnarToOXvkzC9dJqeDpZY+MLQ2Bnxcn1iPSNpYaIyIDqNVrM25KKtOwyONlY4ssXg9DRQSF2LCKTxFJDRGQggiBg6a5z+PlCIawspPh8RiC6d7QTOxaRyWKpISIykH//ehVfn8iCRAKsfGYQAn2dxY5EZNJYaoiIDGBHcg4+3nsJAPD+mH4Y2d9D5EREpo+lhohIzw6lF+OdHacBALNGdEXYMF9xAxGZCZYaIiI9OpenwqubU1CvFTDOX4l3QnuLHYnIbLDUEBHpSY6lHV7YeBKVtfUY2tUFHz09EFIpJ9cjaissNUREelBWXYcwn1EorqhFbw97rJsRACsLmdixiMwKSw0R0QOqUWsw86skXLXqgE6OCmwMHwIHhaXYsYjMDksNEdEDEAQB7+48g5PXb8JeU4svXwxCJ0drsWMRmSWWGiKiB7DhcAZ2puZCJpVgbc4+9HS3FzsSkdliqSEiaqVD6cX4548XAADvje6Dh6tyRU5EZN5YaoiIWuF6SRXmfpMKrQBMCfTiXDRERoClhohIR5W19Zj5VRJUt9QY5OOEv4/vD4mEt24TiY2lhohIB1qtgDe3pSG9qBLuDlZY9zxv3SYyFiw1REQ6WJGQjn3nCyG3kGLd9EB0dFCIHYmIbmOpISJqoZ/O5OOzhHQAQNSEAfD3dhI3EBE1wVJDRNQCF/LL8da3pwAALw3vgkkBXiInIqL/xVJDRHQfN6rqMPOrJFTXaTC8uysWPcmHVBIZI5YaIqJ7UGu0mPN1CnJu3oKPsw1WPzsIFjL+6iQyRvzJJCK6h3/suYBj10phK5fh87BAONnIxY5ERHfBUkNEdBfbk7IRc/Q6AGD5VH8+AoHIyLHUEBE1IyXrJiJjzwIA3gzpidB+HiInIqL7YakhIvofheU1mL0pGXUaLUb288C8P3UXOxIRtQBLDRHRH9SoNXhlUzKKKmrRy90en07xg1TKRyAQtQcsNUREtwmCgMWxZ3EquwxONpZYPyMQtlYWYsciohZiqSEiuu2LI9exIyUHMqkE0c8Oho+LjW47mDbNMMGIqEVYaoiIABxOL8E/9pwHAESO7oOHu7vqvhOWGiJRsdQQkdnLLK3CnG9SoBWAyQFeeGGYr9iRiKgVWGqIyKxV1tZj5ldJUN1SY5CPEz6Y0B8SCS8MJmqPWGqIyGxptQIitqXhcmEl3B2ssO75AFhZyMSORUStxFJDRKZhyxadN1mZkI7/ni+EXCbF2ucD0NFBYYBgRNRWWGqIyDToWGriz+ZjZUI6AOCfEwdgkE8HQ6QiojbUqlITHR0NX19fKBQKBAcHIzEx8a7rqtVqLFu2DN26dYNCoYCfnx/i4+ObrHPw4EGMGTMGSqUSEokEcXFxd+zjnXfewYABA2BrawulUokZM2YgLy+vNfGJyMxdLChHxPZTAIAXH+6CpwO8RE5ERPqgc6nZtm0bIiIisHTpUqSkpMDPzw+hoaEoKipqdv3IyEisW7cOq1atwvnz5zF79mxMmDABqampjetUVVXBz88P0dHRze6juroaKSkpeO+995CSkoKdO3fi0qVLGDt2rK7xicjM3ayqw8yvklBdp8HD3V3w7qjeYkciIj2RCIIg6LJBcHAwhgwZgtWrVwMAtFotvL29MW/ePCxcuPCO9ZVKJRYvXow5c+Y0Lps0aRKsra2xefPmOwNJJIiNjcX48ePvmePkyZMICgpCZmYmfHx87pu7vLwcjo6OUKlUcHBwuO/6RNTOjB0L7Np1z1XqNVrM+CIRR6+WwsfZBt/PeRgdbOVtFJCIWkOXv986fVJTV1eH5ORkhISE/L4DqRQhISE4duxYs9vU1tZCoWh68Z21tTUOHz6sy6HvoFKpIJFI4OTkdNfjlpeXN3kRkXn7x48XcPRqKWzkMqyfEchCQ2RidCo1JSUl0Gg0cHd3b7Lc3d0dBQUFzW4TGhqK5cuXIz09HVqtFvv27cPOnTuRn5/f6tA1NTV45513MG3atLu2tqioKDg6Oja+vL29W308Imr/tidlY+OR6wCA5VP80cvDXtxARKR3Br/7aeXKlejRowd69+4NuVyOuXPnIjw8HFJp6w6tVqsxZcoUCIKANWvW3HW9RYsWQaVSNb6ys7NbOwQiaudSsm4iMvYsAOCNkB4Y2d9D5EREZAg6NQtXV1fIZDIUFhY2WV5YWAgPj+Z/Sbi5uSEuLg5VVVXIzMzExYsXYWdnh65du+oc9rdCk5mZiX379t3z3JqVlRUcHByavIjI/BSW12D2pmTUabQI7eeO1/7UQ+xIRGQgOpUauVyOgIAAJCQkNC7TarVISEjA0KFD77mtQqGAp6cn6uvrsWPHDowbN06noL8VmvT0dPz8889wcXHRaXsiMj81ag1e2ZSMoopa9HK3x6dT/CGV8hEIRKbKQtcNIiIiEBYWhsDAQAQFBWHFihWoqqpCeHg4AGDGjBnw9PREVFQUAODEiRPIzc2Fv78/cnNz8f7770Or1WLBggWN+6ysrMSVK1cav87IyEBaWhqcnZ3h4+MDtVqNp59+GikpKdi9ezc0Gk3jNTzOzs6Qy3mxHxE1JQgCFseexansMjjZWGL9jEDYWen8K4+I2hGdf8KnTp2K4uJiLFmyBAUFBfD390d8fHzjxcNZWVlNrpepqalBZGQkrl27Bjs7O4waNQqbNm1qctdSUlISHn/88cavIyIiAABhYWGIiYlBbm4udt2+VdPf379Jnv379+Oxxx7TdRhEZOK+OHIdO1JyIJNKEP3sYPi42IgdiYgMTOd5atorzlNDZOL+ME/N4fQShG1MhEYrYMlTffHi8C4ihyOi1jLYPDVERMYus7QKc75JgUYr4OkAL4Q/7Ct2JCJqIyw1RGQyKmvrMfOrJKhuqeHv7YQPxveHRMILg4nMBUsNEZkELYC3tqfhcmElOtpbYd30ACgsZWLHIqI2xFJDRCbhM9cA7D1XCLlMirXTA+DuoLj/RkRkUlhqiKjd23e+ECs6BgIAPpjQH4N9OoiciIjEwFJDRO3ajao6LNp5GgDwwjBfTAnkc96IzBVLDRG1a0t3nUNJZR161tzAolG9xY5DRCJiqSGidiv+bD5+OJUHmVSCT/L2w8qCFwYTmTOWGiJql25W1SEyruHJ27Me7YqBNSUiJyIisbHUEFG79P4PDaedenS0w+shfPI2EbHUEFE7tPdcAb5Py4NUAnwy2Y+nnYgIAEsNEbUzZdV1WBx7+7TTiG7w83YSNxARGQ2WGiJqV/72w3mUVNaie0c7vP5nnnYiot+x1BBRu7HvfCFiU3MhlQAfPz2Qj0EgoiZYaoioXSirrsO7sWcAADMf7YpBnDWYiP4HSw0RtQvLfjiP4opadHOzxZshPcWOQ0RGiKWGiIzez+cLsfO3006T/XjaiYiaxVJDREZNVa1uPO308iNd+bBKIrorlhoiMmrLdp9HUUUturrZIuIvPO1ERHfHUkNERuuXi4XYkZIDiQT4+GmediKie2OpISKjpLqlxqKdt087De+CgM487URE98ZSQ0RG6YPd51FYXouurrZ464leYschonaApYaIjM7+i0X4NrnhtNNHnGSPiFqIpYaIjMofTzu9+HAXBPo6i5yIiNoLlhoiMir/2HMeBeU18HWxwdu6nHaaNs1woYioXWCpISKj8eulImxP+u20kx+s5TqcdmKpITJ7LDVEZBTKa34/7fTCMF8EdeFpJyLSDUsNERmFf+65gHxVDTq72GB+KO92IiLdsdQQkegOXi7G1pPZjZPs2cgtxI5ERO0QSw0Riaq8Ro2FO04DAMKG8rQTEbUeSw0RiSrqxwvIU9XAx9kGC0bytBMRtR5Ljb5s2SJ2AqJ25+DlYmxJzAbQMMkeTzsR0YNgqdEXlhoinVT84W6nsKGd8VBXF5ETEVF7x1JDRKKI+ukicstuwdvZGgtG9hY7DhGZAJYaImpzh9NL8M2JLADAvyYNhK0VTzsR0YNjqSGiNlVZW493bt/tNP2hzhjWzVXkRERkKlhqiKhNffjTBeSW3YJXB2ssfJKnnYhIf1hqiKjNHL1Sgs3HG047fcTTTkSkZyw1RNQmKmvrMf+7htNOzz/kg2HdedqJiPSLpYaI2sS/bt/t5OlkjYVP9hE7DhGZIJYaIjK4o1dLsOl4JoCGSfbseNqJiAyApYaI9OMuE1BW/eFup2eDffAwTzsRkYGw1BCRftyl1HwUfxHZNxpOOy3i3U5EZEAsNURkMMevleLLYw2nnT6cNAD2CkuRExGRKWOpISKDqK6rx4LbdztNC/LGIz3cRE5ERKaOpYaIDOKj+EvIulENpaMC747i3U5EZHgsNUSkd8evlSLm6HUAwIeTBvK0ExG1CZYaItKr6rrf73Z6Zog3Hu3J005E1DZYaohIrz7eewmZpdXo5KjAu6N52omI2g5LDRHpTWLGjcbTTlETB8CBp52IqA2x1DwgQRCw8UgG/u4+VOwoRKK6JbHAgu9OQRCAKYFeeKxXR7EjEZGZYal5QGdyVfjbD+exwWUgfjyTL3YcItF80nEIrpdWw8NBgcWj+4odh4jMEEvNAxro5YRZI7oCAN757jSySqtFTkTU9k5ev4EvnAcAAKImDYCjNU87EVHba1WpiY6Ohq+vLxQKBYKDg5GYmHjXddVqNZYtW4Zu3bpBoVDAz88P8fHxTdY5ePAgxowZA6VSCYlEgri4uDv2IwgClixZgk6dOsHa2hohISFIT09vTXy9e/uJXgioLkBFbT3mbklBbb1G7EhEbaZAVYPXtqRCkEgwOcALj/O0ExGJROdSs23bNkRERGDp0qVISUmBn58fQkNDUVRU1Oz6kZGRWLduHVatWoXz589j9uzZmDBhAlJTUxvXqaqqgp+fH6Kjo+963I8++gifffYZ1q5dixMnTsDW1hahoaGoqanRdQh6ZymT4rOcBDhaW+J0jgof/nRR7EhEbaKyth7hMSeRr6pB19qbiHyKp52ISDwSQRAEXTYIDg7GkCFDsHr1agCAVquFt7c35s2bh4ULF96xvlKpxOLFizFnzpzGZZMmTYK1tTU2b958ZyCJBLGxsRg/fnzjMkEQoFQq8dZbb+Htt98GAKhUKri7uyMmJgbPPPPMfXOXl5fD0dERKpUKDg4Ougy5ZcaOxc8frsfLXyUBANZND0BoPw/9H4fISKg1Wrz0ZRIOXi6Gq50csakx8N7xtdixiMjE6PL3W6dPaurq6pCcnIyQkJDfdyCVIiQkBMeOHWt2m9raWigUiibLrK2tcfjw4RYfNyMjAwUFBU2O6+joiODg4Hset7y8vMnL0EL6uuPl4V0AAPO/PYXsG7y+hkyTIAiIjD2Lg5eLYW0pwxcvDIG3ukLsWERk5nQqNSUlJdBoNHB3d2+y3N3dHQUFBc1uExoaiuXLlyM9PR1arRb79u3Dzp07kZ/f8juFftu3LseNioqCo6Nj48vb27vFx3sQC0b2hp+3E8pr6jFvSyrq6rVtclyithS9/wq2JWVDKgFWTRuEgV5OYkciIjL83U8rV65Ejx490Lt3b8jlcsydOxfh4eGQSg176EWLFkGlUjW+srOzDXq838gtpFg9bRAcFBZIyy7Dx3t5fQ2ZltjUHHzy38sAgL+N7YeQvu732YKIqG3o1CxcXV0hk8lQWFjYZHlhYSE8PJq/fsTNzQ1xcXGoqqpCZmYmLl68CDs7O3Tt2rXFx/1t37oc18rKCg4ODk1ebcXb2QYfT/YDAKw/lIGEC4X32YKofTh6pQQLvmt4rtOsR7ti+lBfcQMREf2BTqVGLpcjICAACQkJjcu0Wi0SEhIwdOi9Z9RVKBTw9PREfX09duzYgXHjxrX4uF26dIGHh0eT45aXl+PEiRP3Pa5YQvt54IVhvgCAt749hbyyW+IGInpAlwoqMGtzMtQaAaMHdsI7I3uLHYmIqAmdzwFFRERg/fr1+PLLL3HhwgW8+uqrqKqqQnh4OABgxowZWLRoUeP6J06cwM6dO3Ht2jUcOnQII0eOhFarxYIFCxrXqaysRFpaGtLS0gA0XBiclpaGrKwsAA13RL3xxhv44IMPsGvXLpw5cwYzZsyAUqlscpeUsVk0qjcGeDqirFqN17akQq3h9TXUPhWW1yB8YyIqauoxxLcDPp3sB6lUInYsIqImLHTdYOrUqSguLsaSJUtQUFAAf39/xMfHN17Em5WV1eR6mZqaGkRGRuLatWuws7PDqFGjsGnTJjg5OTWuk5SUhMcff7zx64iICABAWFgYYmJiAAALFixAVVUVXnnlFZSVlWH48OGIj4+/484qY2JlIcPqZwfhqc8OIynzJpbvu8z/u6V2p7K2Hi/GnESeqgZd3WyxfkYgFJYysWMREd1B53lq2qu2mKcGu3Y1+097TudjzjcpAICY8CF80B+1G/W356I5cHsump2vPgwfF5vmV77HzwARUWsZbJ4aap3RAzvh+Yd8AAAR20+hQCX+LMhE9yMIAiLjzuLA5WIoLKXYEDbk7oWGiMgIsNS0kcjRfdG3kwNuVNXhta2pqOf1NWTk/v3rVWw9+dtcNIPh5+0kdiQiontiqWkjCksZop8bDFu5DIkZN7AywTgexknUnLjUXHy89xIA4P2x/fAXzkVDRO0AS00b6uJqi39OHAAAWL3/Cg6nl4iciOhOR6+WYP53pwAArzzaFTNaOhfNtGmGC0VE1AIsNW1snL8npgV5QxCAN7aloqiC19eQ8bhcWIFZm27PRTOgExbqcrceSw0RiYylRgRLx/RDbw97lFTW4Y2tadBozeIGNDJyReU1CN94EhU19Qjs3AGfTuFcNETUvrDUiEBhKcPqZwfDRi7D0aulWP3LFbEjkZmrqq1HeMxJ5JbdQhdXzkVDRO0TS41Iune0wwfj+wMAViZcxrGrpSInInNVr9Fi7jcpOJdXDhdbOWLCh6CDrVzsWEREOmOp0ZdWXE8wcbAXJgd4QSsAr29NRUllrQGCEd2dIAh47/tz2H+pYS6az8MC0dnFVuxYREStwlKjL628SPJv4/qhR0c7FFXU4s1tadDy+hpqQ2sOXMWWxCxIJMBnzwzCIJ8OYkciImo1lhqR2cgtEP3cYCgspTiUXoI1B66KHYnMxPdpufgo/vZcNGP64Yl+HiInIiJ6MCw1RqCnuz2WjW24vubT/15CYsYNkRORqTt+rRTzvz0NAHh5eBeEDfMVNxARkR6w1BiJyYFemDDIE1oBeG1LKm5U1YkdiUzUlaIKvPJVEuo0WjzZ3wPvjuojdiQiIr1gqTESEokEH4zvj65utigor0HEdh2vr9myxXDhqH1owfdAUUUNwr44ifKaegR07oD/m+rPuWiIyGSw1BgRWysLRD87GFYWUvx6qRjrD11r+cYsNXSf74Gq2nq8yLloiMiEsdQYmT6dHLB0TD8AwEd7LyE586bIicgU1Gu0mLclFWdzf5+Lxplz0RCRiWGpMULTgrwxxk8JjVbAa1tSUVbN62uo9QRBwNJd5/DLxSJYWUixnnPREJGJYqkxQhKJBP+c0B++LjbILbuFt789DUHg/DXUOmsPXMPXJxrmoln5zCAM5lw0RGSiWGqMlL3CEqufHQy5TIqfLxRiw+EMsSNRO7TrVB7+FX8RALDkqb4Y2Z9z0RCR6WKpMWL9PR0R+VTD7bb/ir+ItOwycQNRu3LiWine3n4KAPDS8C4If7iLyImIiAyLpcbITX+oM0YN8IBaI2DuNylQ3VKLHYnagStFFZj5h7loFnMuGiIyAyw1Rk4ikeDDSQPh7WyNnJu3sOC7U7y+hu6pqKIGL2xsmItmsI8T56IhIrPBUtMOOCgsEf3sYFjKJNh7rhBfHr0udiQyUtV19XgpJgk5N2/B18WGc9EQkVlhqWknBno5YdGTDacQ/vnjRZzJUYmciIxNPSSY900qzuSq4GwrR0x4EFzsrMSORUTUZlhq2pHwh33xRF931Gm0mLslBRU1vL6GGgiCgPc9HkbCb3PRzAiEryvnoiEi88JS045IJBJ8/LQfPJ2skVlajYU7z/D6GgLQMBfNZud+t+ei8UdAZ85FQ0Tmh6WmnXG0scTqZwfBQirBntP5+PpEltiRSESCIOCTvZca56KJHN0XI/t3EjkVEZE4WGraoUE+HfDOyN4AgGW7z+NcHq+vMUdqjRbv7DiN1fuvAADeLDqJl4ZzLhoiMl8sNe3Uy490wZ97d0RdvRZzv0lFpdRS7EjUhqrr6vHKV0nYnpQDqQSImjgAr5ekiB2LiEhULDXtlEQiwSeT/dDJUYGMkios7vQIr68xE6WVtZi2/gT2XyqGwlKK/0wPxLQgH7FjERGJjqWmHetgK8eqaYMgk0rwvWMP/N++y2JHIgPLKq3G02uP4VR2GTrYWOKbmQ8hpK+72LGIiIwCS007F+jrjCVP9QUAfPbLFaxKSBc5ERnK2VwVJq45ioySKng6WeO7V4fxidtERH/AUmMCwob5YlHhcQDAp/suY+2BqyInIn07eLkYU9cdQ0llLfp0ckDsX4ehm5ud2LGIiIwKS42JmFV6Cm8/0RMA8OFPF7HhcIbIiUhfYlNz8GLMSVTVafBwdxdsn/UQOjooxI5FRGR0LMQOQPoz9089UKcR8FlCOv6++zwsZRLMGOordixqJUEQsO7gNXz4U8McNGP9lPhksh/kFvx/ESKi5rDUmJg3Q3pArdFiza9XseT7c7CUSXlnTDuk1QpYtvs8Ym4/vHTmI12w6Mk+fNo2EdE9sNSYGIlEggWhvaCu1+Lzwxl4N/YMLKQSTA70FjsatVCNWoO3tp/CnjP5AIDI0X3w8iNdRU5FRGT8WGpMkEQiweLRfVCvFRBz9DoW7DgNS5kU4wd5ih2N7kN1S41XvkrCiYwbsJRJ8OkUf4z1U4odi4ioXeDJeRMlkUiwdExfPBfsA0EAIranYc/pfLFjmbYtWx5o8wJVDaauO4YTGTdgZ2WBL8ODWGiIiHTAUmPCJBIJ/j6uP6YEekErAK9tTUX82QKxY5muByg1V4oqMPHfR3CxoAId7a2wfdZQDOvuqsdwRESmj6XGxEmlEkRNHIiJgzyh0QqYtyUFCRcKxY5Ff3Dy+g1MWnMMeaoadHWzxY5Xh6Gv0kH3HU2bpv9wRETtCEuNGZBJJfh4sh/G+Cmh1gh4dXMKfr1UJHYsAhB/tgDPf34CqltqDPJxwo7Zw+DtbNO6nbHUEJGZY6kxEzKpBMun+OHJ/h6o02jxyqZkHE4vETuWWdt0PBN//ToZtfVahPTpiG9efggdbOVixyIiardYasyIpUyKlc8MQkgfd9TVa/HyVydx/Fqp2LHMjiAI+GTvJbwXdxZaAZgW5I21zwfAWi4TOxoRUbvGUmNm5BZSRD83CI/1ckONWosXY04i6foNsWOZDbVGi3d2nMbq/VcAAG+E9MA/JwyAhYw/ikRED4q/Sc2QlYUMa58PwCM9XFFdp8ELG08iNeum2LFMXnVdPV75Kgnbk3IglQBREwfgjZCekEg4SzARkT6w1JgphaUM/5keiIe6OqOyth4zvkjEmRyV2LFMVmllLaatP4H9l4qhsJTiP9MD+fgKIiI9Y6kxY9ZyGTaEDcEQ3w6oqKnH8xtO4Hxeeet29oATz5myrNJqPL32GE5ll6GDjSW+mfkQQvq6ix2LiMjksNSYilbezmtrZYGN4UEY5OME1S01nt9wApcKKnTfEUtNs87kqDBxzRFklFTB08ka3706DIN9Oogdi4jIJLHUmIoHmKPEzsoCX74YhIFejrhRVYfnPj+OK0WVegxnng5eLsYz/zmGkso69OnkgNi/DkM3NzuxYxERmSyWGgIAOCgssenFYPTt5ICSyjo8u/44MkqqxI7VbsWm5uDFmJOoqtPg4e4u2D7rIXR0UIgdi4jIpLHUUCNHG0tsfjkYvT3sUVRRi2n/OY7MUhYbXQiCgLUHruLNbadQrxUw1k+JjS8EwV5hKXY0IiKT16pSEx0dDV9fXygUCgQHByMxMfGu66rVaixbtgzdunWDQqGAn58f4uPjdd5nQUEBpk+fDg8PD9ja2mLw4MHYsWNHa+LTPTjbyrH55WB072iHgvIaPLv+BHJuVosdq13QAvjbD+fx4U8XAQAzH+mCFVP9Ibfg/zsQEbUFnX/bbtu2DREREVi6dClSUlLg5+eH0NBQFBU1/yyhyMhIrFu3DqtWrcL58+cxe/ZsTJgwAampqTrtc8aMGbh06RJ27dqFM2fOYOLEiZgyZUqT/ZB+uNpZ4ZuXg9HV1Ra5Zbcwbf1x5KtuiR3LqNWoNZjnGYKYo9cBAJGj+2Dx6L6QSjkHDRFRmxF0FBQUJMyZM6fxa41GIyiVSiEqKqrZ9Tt16iSsXr26ybKJEycKzz33nE77tLW1Fb766qsm+3F2dhbWr1/fotwqlUoAIKhUqhatT4KQX3ZLePSjX4TO7+wWHvt4v1CgunX3lceMabtgRiYt66bw9JojQud3dgvd390jfJ+WK3YkIiKTocvfb50+qamrq0NycjJCQkIal0mlUoSEhODYsWPNblNbWwuFoukFktbW1jh8+LBO+xw2bBi2bduGGzduQKvVYuvWraipqcFjjz121+OWl5c3eZFuPBwV+GbmQ/DqYI2Mkio8u/44iitqxY5lFARBwOH0Ejz3+XGMiz6Ck9dvwk5Thy/DgzDWTyl2PCIis6RTqSkpKYFGo4G7e9OJw9zd3VFQUNDsNqGhoVi+fDnS09Oh1Wqxb98+7Ny5E/n5+Trtc/v27VCr1XBxcYGVlRVmzZqF2NhYdO/evdnjRkVFwdHRsfHl7e2ty1DpNk8na2yZ+RA6OSpwtbgKz31+HKWVRlps2mCuHK1WwE9n8jEu+gie33ACR66UQiaVYOJgT/yQsQPDursaPAMRETXP4Fcwrly5Ej169EDv3r0hl8sxd+5chIeHQyrV7dDvvfceysrK8PPPPyMpKQkRERGYMmUKzpw50+z6ixYtgkqlanxlZ2frYzhmydvZBltmPoSO9la4XFiJ5zckoqy6TuxYdzJgqamr12J7UjZC/u8AXv06BadzVFBYSvHCMF8cmP8Ylk/xR5c6fhpIRCQmC11WdnV1hUwmQ2FhYZPlhYWF8PDwaHYbNzc3xMXFoaamBqWlpVAqlVi4cCG6du3a4n1evXoVq1evxtmzZ9GvXz8AgJ+fHw4dOoTo6GisXbv2juNaWVnByspKl+HRPfi62uKbmQ/hmf8cx4X8ckzfkIjNLwfD0dq0b1Wuqq3H1pPZ+PzQNeSragAADgoLhA3zxQvDfOFix+8xIiJjodPHJXK5HAEBAUhISGhcptVqkZCQgKFDh95zW4VCAU9PT9TX12PHjh0YN25ci/dZXd1wS/H/frojk8mg1Wp1GQI9gO4d7fDNzGA428pxJleFsC8SUVGjFjuWQdysqsOKny/j4X/9gr/vPo98VQ062lvh3VG9cWThn/DWE71YaIiIjIxOn9QAQEREBMLCwhAYGIigoCCsWLECVVVVCA8PB9Bw67WnpyeioqIAACdOnEBubi78/f2Rm5uL999/H1qtFgsWLGjxPnv37o3u3btj1qxZ+OSTT+Di4oK4uDjs27cPu3fv1sf7QC3U090em18KxrOfH0dadhnCN57Ely8GwVbsYHqSr7qFzw9lYEtiFqrrNAAAXxcbzBrRDRMGeUJhKRM5IRER3Y3OpWbq1KkoLi7GkiVLUFBQAH9/f8THxzde6JuVldXkE5WamhpERkbi2rVrsLOzw6hRo7Bp0yY4OTm1eJ+Wlpb48ccfsXDhQowZMwaVlZXo3r07vvzyS4waNeoB3wLSVV+lQ0OxWX8cSZk3ER5zEjESC9iIHewBXC2uxLoDVxGbmgu1RgAA9O3kgL8+3g1P9u8EGeebISIyehJBEASxQ7SF8vJyODo6QqVSwcHBQew4JuFUdhme//wEKmrr8XBlDhYtega9POxhKRNpBt2xY4Fdu3Ta5HROGdb8ehXx5wrw209CcBdnvPpYN4zo6QaJRIcy04rjExHRveny91vnT2qIfuPn7YSYF4dgxoZEHLHzwlOrDkNhKUV/pSP8vJ3g5+0Efy8neDtb61YODEwQBBy9Woo1v17F4SsljctD+rjj1ce6IaBzBxHTERFRa7HU0AMJ6OyMTS8H4/8+2oY01y6oqKlHUuZNJGXebFzH2VYOP6+mRaeDrbzNs2q1Av57vhBrfr2CUzkqAIBMKsE4PyVmjeiGXh72bZ6JiIj0h6WGHthgnw7YlPUjtJ99j4zSKpzKLsOp7DKk5ahwIa8cN6rqsP9SMfZfKm7cprOLDfy8nOB/u+j0UzoY7CLcunotvk/LxdoDV3G1uOGp41YWUjwzxBsvP9IV3s7t+WogIiL6DUsN6Y1UKkE3Nzt0c7PDxMFeAIDaeg0u5Fc0lJzbZedaSRUyS6uRWVqNXafyAAAWUgn6dHKAn7djY9np5mb3QA+ErK6rx9bEhjlm8m7PMWOvsEDYUF+88LAvXHlLNhGRSWGpIYOyspDB37uhpITdXqaqVuN0bhnSsspwKqeh7JRU1uFMrgpnclXYjCwAgJ2VBQb+dtrKywmDfJzg7qC4+8FuK6uuw5dHMxFzNAM3qxvm0XGzt8LLw7vg2WAf2CtMe8JAIiJzxVJDbc7RxhKP9HDDIz3cADRcuJunqmlScs7kqFBZW4+jV0tx9Gpp47YeDoqGT3NuF6UBno6NJaXAwgaf7z6Pb/4wx0xnFxvMerQbJg7mHDNERKaOpYZEJ5FI4OlkDU8na4we2AkAUK/RIr2osuH6nJwypGaV4XJhBQrKa1BwrgZ7zxXe3hbo7mYHH2cbHOzxLNSHMwAAfTo54NXHumFUfw9YiHWLORERtSmWGjJKFjIp+nRyQJ9ODngmyAdAwzUyZ3PLG6/PScsuQ27ZLaQXVSK9qBKQyBB0e46Zx3SdY4aIiNo9lhpqN2zkFgjq4oygLs6Ny4oranE6pwzpRZUY8u8PEfDhehETEhGRmFhqqF1zs7fCn/u448993IFPC++/gSFNmybu8YmIzBwvNiDSF5YaIiJRsdQQERGRSWCpISIiIpPAUkNEREQmgaWGiIiITAJLDREREZkElhoiIiIyCSw1pB+8nZmIiETGUkP6wVJDREQiY6khIiIik8BSQ0RERCaBpYaIiIhMAksNmQ5e10NEZNZYash0sNQQEZk1lhoiIiIyCSw1REREZBJYaoiIiMgksNQQERGRSWCpISIiIpPAUkNEREQmgaWGiIiITAJLDREREZkElhoiIiIyCSw1REREZBJYaoiIiMgksNQQERGRSWCpISIiIpNgIXaAtiIIAgCgvLxc5CRERETUUr/93f7t7/i9mE2pqaioAAB4e3uLnISIiIh0VVFRAUdHx3uuIxFaUn1MgFarRV5eHuzt7SGRSMSOo1fl5eXw9vZGdnY2HBwcxI4jCnN/Dzh+8x4/wPfA3McPmO57IAgCKioqoFQqIZXe+6oZs/mkRiqVwsvLS+wYBuXg4GBS38itYe7vAcdv3uMH+B6Y+/gB03wP7vcJzW94oTARERGZBJYaIiIiMgksNSbAysoKS5cuhZWVldhRRGPu7wHHb97jB/gemPv4Ab4HgBldKExERESmjZ/UEBERkUlgqSEiIiKTwFJDREREJoGlhoiIiEwCS42Rio6Ohq+vLxQKBYKDg5GYmHjP9b/99lv07t0bCoUCAwYMwI8//tj4b2q1Gu+88w4GDBgAW1tbKJVKzJgxA3l5eYYeRqvpc/z/a/bs2ZBIJFixYoWeU+uXId6DCxcuYOzYsXB0dIStrS2GDBmCrKwsQw3hgeh7/JWVlZg7dy68vLxgbW2Nvn37Yu3atYYcwgPRZfznzp3DpEmT4Ovre8/vbV3fU7Hp+z2IiorCkCFDYG9vj44dO2L8+PG4dOmSAUfwYAzxPfCbDz/8EBKJBG+88YZ+Q4tNIKOzdetWQS6XC1988YVw7tw5YebMmYKTk5NQWFjY7PpHjhwRZDKZ8NFHHwnnz58XIiMjBUtLS+HMmTOCIAhCWVmZEBISImzbtk24ePGicOzYMSEoKEgICAhoy2G1mL7H/0c7d+4U/Pz8BKVSKfzf//2fgUfSeoZ4D65cuSI4OzsL8+fPF1JSUoQrV64I33///V33KSZDjH/mzJlCt27dhP379wsZGRnCunXrBJlMJnz//fdtNawW03X8iYmJwttvvy1s2bJF8PDwaPZ7W9d9is0Q70FoaKiwceNG4ezZs0JaWpowatQowcfHR6isrDTwaHRniPH/cV1fX19h4MCBwuuvv26YAYiEpcYIBQUFCXPmzGn8WqPRCEqlUoiKimp2/SlTpgijR49usiw4OFiYNWvWXY+RmJgoABAyMzP1E1qPDDX+nJwcwdPTUzh79qzQuXNnoy41hngPpk6dKjz//POGCaxnhhh/v379hGXLljVZZ/DgwcLixYv1mFw/dB3/H93te/tB9ikGQ7wH/6uoqEgAIBw4cOBBohqEocZfUVEh9OjRQ9i3b58wYsQIkys1PP1kZOrq6pCcnIyQkJDGZVKpFCEhITh27Fiz2xw7dqzJ+gAQGhp61/UBQKVSQSKRwMnJSS+59cVQ49dqtZg+fTrmz5+Pfv36GSa8nhjiPdBqtdizZw969uyJ0NBQdOzYEcHBwYiLizPYOFrLUN8Dw4YNw65du5CbmwtBELB//35cvnwZTzzxhGEG0kqtGb8Y+zSktsqrUqkAAM7Oznrbpz4Ycvxz5szB6NGj7/h5MRUsNUampKQEGo0G7u7uTZa7u7ujoKCg2W0KCgp0Wr+mpgbvvPMOpk2bZnQPPTPU+P/1r3/BwsICr732mv5D65kh3oOioiJUVlbiww8/xMiRI/Hf//4XEyZMwMSJE3HgwAHDDKSVDPU9sGrVKvTt2xdeXl6Qy+UYOXIkoqOj8eijj+p/EA+gNeMXY5+G1BZ5tVot3njjDTz88MPo37+/XvapL4Ya/9atW5GSkoKoqKgHjWi0zOYp3dRArVZjypQpEAQBa9asETtOm0hOTsbKlSuRkpICiUQidhxRaLVaAMC4cePw5ptvAgD8/f1x9OhRrF27FiNGjBAzXptYtWoVjh8/jl27dqFz5844ePAg5syZA6VSabL/10p3N2fOHJw9exaHDx8WO0qbyM7Oxuuvv459+/ZBoVCIHcdgWGqMjKurK2QyGQoLC5ssLywshIeHR7PbeHh4tGj93wpNZmYmfvnlF6P7lAYwzPgPHTqEoqIi+Pj4NP67RqPBW2+9hRUrVuD69ev6HcQDMsR74OrqCgsLC/Tt27fJOn369DG6X+qGGP+tW7fw7rvvIjY2FqNHjwYADBw4EGlpafjkk0+MqtS0Zvxi7NOQDJ137ty52L17Nw4ePAgvL68H3p++GWL8ycnJKCoqwuDBgxuXaTQaHDx4EKtXr0ZtbS1kMtkD5TYGPP1kZORyOQICApCQkNC4TKvVIiEhAUOHDm12m6FDhzZZHwD27dvXZP3fCk16ejp+/vlnuLi4GGYAD8gQ458+fTpOnz6NtLS0xpdSqcT8+fOxd+9eww2mlQzxHsjlcgwZMuSO21cvX76Mzp0763kED8YQ41er1VCr1ZBKm/7Kk8lkjZ9iGYvWjF+MfRqSofIKgoC5c+ciNjYWv/zyC7p06aKPuHpniPH/+c9/xpkzZ5r8HgwMDMRzzz2HtLQ0kyg0AHhLtzHaunWrYGVlJcTExAjnz58XXnnlFcHJyUkoKCgQBEEQpk+fLixcuLBx/SNHjggWFhbCJ598Ily4cEFYunRpk9tZ6+rqhLFjxwpeXl5CWlqakJ+f3/iqra0VZYz3ou/xN8fY734yxHuwc+dOwdLSUvjPf/4jpKenC6tWrRJkMplw6NChNh/f/Rhi/CNGjBD69esn7N+/X7h27ZqwceNGQaFQCP/+97/bfHz3o+v4a2trhdTUVCE1NVXo1KmT8PbbbwupqalCenp6i/dpbAzxHrz66quCo6Oj8Ouvvzb5PVhdXd3m47sfQ4z/f5ni3U8sNUZq1apVgo+PjyCXy4WgoCDh+PHjjf82YsQIISwsrMn627dvF3r27CnI5XKhX79+wp49exr/LSMjQwDQ7Gv//v1tNCLd6HP8zTH2UiMIhnkPNmzYIHTv3l1QKBSCn5+fEBcXZ+hhtJq+x5+fny+88MILglKpFBQKhdCrVy/h008/FbRabVsMR2e6jP9uP+MjRoxo8T6Nkb7fg7v9Hty4cWPbDUoHhvge+CNTLDUSQRCENvpQiIiIiMhgeE0NERERmQSWGiIiIjIJLDVERERkElhqiIiIyCSw1BAREZFJYKkhIiIik8BSQ0RERCaBpYaIiIhMAksNERERmQSWGiIiIjIJLDVERERkElhqiIiIyCT8P/aGEEIWqGlLAAAAAElFTkSuQmCC"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 10
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:10:11.655407Z",
     "start_time": "2024-08-26T16:10:11.579241Z"
    }
   },
   "cell_type": "code",
   "source": [
    "val_avg = np.mean(val_err_l2, axis=0)\n",
    "l2_est = l2_list[np.argmin(val_avg)]\n",
    "print(l2_est)"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.04\n"
     ]
    }
   ],
   "execution_count": 11
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "We can also use the one standard error rule."
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:10:11.781466Z",
     "start_time": "2024-08-26T16:10:11.692294Z"
    }
   },
   "cell_type": "code",
   "source": "parameter_tuning.get_lambda_one_se_1d(val_err_l2.squeeze(), l2_list)",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.05"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "execution_count": 12
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Alternative: use grid search for two hyperparameters\n",
    "\n",
    "As discussed above, we can use 2D grid search for different combinations of $\\lambda_1$ and $\\lambda_2$.\n",
    "This will usually take long time, and is not recommended unless the network (node number and sample size) is small."
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:10:11.890324Z",
     "start_time": "2024-08-26T16:10:11.814212Z"
    }
   },
   "cell_type": "code",
   "source": [
    "RUN_ME = False\n",
    "\n",
    "if RUN_ME:\n",
    "    l1_list = np.arange(0.02, 0.401, 0.02)\n",
    "    l2_list = np.arange(0.01, 0.16, 0.01)\n",
    "    val_err_2d = parameter_tuning.cv_2d(\n",
    "        dat1,\n",
    "        dat2,\n",
    "        n_cv=5,\n",
    "        ratio_val=0.5, \n",
    "        lambda1_lst=l1_list,\n",
    "        lambda2_lst=l2_list,\n",
    "        dep_mat=dep_mat,\n",
    "        cores=cores,\n",
    "    )\n",
    "    parameter_tuning.plot_error_2d(val_err_2d)"
   ],
   "outputs": [],
   "execution_count": 13
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Alternative: use statistical approaches\n",
    "\n",
    "We can also estimate $\\lambda_1$ and $\\lambda_2$ based on some statistics of the data.\n",
    "More details can be found in the DDN3.0 paper. Usually, these estimates are not accurate enough, and we suggest the above cross validation protocol."
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-26T16:10:12.001842Z",
     "start_time": "2024-08-26T16:10:11.924026Z"
    }
   },
   "cell_type": "code",
   "source": [
    "n = int((len(dat1) + len(dat2))/2)\n",
    "p = dat1.shape[1]\n",
    "alpha1 = 0.05\n",
    "alpha2 = 0.01\n",
    "\n",
    "# using the method in the M-B algorithm to estimate lambda1\n",
    "l1_est = parameter_tuning.get_lambda1_mb(alpha1, n, p, mthd=0)\n",
    "\n",
    "# using Bai's method in the UAI paper to estimate lambda2\n",
    "l2_est = parameter_tuning.get_lambda2_bai(dat1, dat2, alpha=alpha2)\n",
    "\n",
    "print(l1_est, l2_est)"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.2711602723929704 0.08082354445935315\n"
     ]
    }
   ],
   "execution_count": 14
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ai",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}