From 065450d8c54ecac5024bda24fea669ef73109ce4 Mon Sep 17 00:00:00 2001 From: SU Xiaochen <164285740+DataboyUsen@users.noreply.github.com> Date: Mon, 27 Apr 2026 20:44:51 +0800 Subject: [PATCH] Adaptive lasso (#10) - feat(adaptive-lasso): add omega parameter for adaptive lasso weights * Change `rho` from scalar scalar to array * Add `omega` parameter in `plqERM_ElasticNet` to apply adaptive lasso - perf(core): optimize adaptive LASSO conditional logic in C++ layer * Reduce redundant computations when omega is not provided - feat(validation): add parameter validation for omega * Validate length matches number of features * Ensure all omega values are positive * Add informative error messages - test(ci): add Adaptive ElasticNet CI tests * Test adaptive weighting with synthetic data - docs(examples): add Adaptive ElasticNet example * Compare adaptive vs standard ElasticNet --- doc/source/example.rst | 1 + doc/source/examples/Adaptive_ElasticNet.ipynb | 201 ++++++++++++++++++ rehline/_base.py | 32 +-- rehline/_class.py | 27 ++- rehline/_internal.pyi | 2 +- src/rehline.cpp | 8 +- src/rehline.h | 63 +++--- tests/test_elastic_net.py | 120 ++++++++++- 8 files changed, 402 insertions(+), 52 deletions(-) create mode 100644 doc/source/examples/Adaptive_ElasticNet.ipynb diff --git a/doc/source/example.rst b/doc/source/example.rst index 7a3dcea..7e52813 100644 --- a/doc/source/example.rst +++ b/doc/source/example.rst @@ -32,3 +32,4 @@ Example Gallery examples/GridSearchCV_reg_losses.ipynb examples/ElasticNet.ipynb + examples/Adaptive_ElasticNet.ipynb diff --git a/doc/source/examples/Adaptive_ElasticNet.ipynb b/doc/source/examples/Adaptive_ElasticNet.ipynb new file mode 100644 index 0000000..4570d14 --- /dev/null +++ b/doc/source/examples/Adaptive_ElasticNet.ipynb @@ -0,0 +1,201 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1b4bedb4", + "metadata": {}, + "source": [ + "## Adaptive ElasticNet\n", + "\n", + "[![Slides](https://img.shields.io/badge/🦌-ReHLine-blueviolet)](https://rehline-python.readthedocs.io/en/latest/)\n", + "\n", + "Adaptive ElasticNet solves the following optimization problem:\n", + "$$\n", + "\\min_{\\beta \\in \\mathbb{R}^d} \\; C \\sum_{i=1}^{n} \\text{PLQ}(y_i, \\mathbf{x}_i^T \\beta) + \\ell_1\\text{ratio} \\sum_{j=1}^{d} \\omega_j |\\beta_j| + \\frac{1}{2}(1 - \\ell_1\\text{ratio})\\|\\beta\\|_2^2, \\quad \\text{s.t.} \\quad \\mathbf{A}\\beta + \\mathbf{b} \\geq \\mathbf{0},\n", + "$$\n", + "\n", + "where\n", + "\n", + "- $\\text{PLQ}(\\cdot)$ is a piecewise linear-quadratic loss function (e.g., hinge, quantile, Huber),\n", + "- $\\mathbf{x}_i \\in \\mathbb{R}^d$ is a feature vector,\n", + "- $y_i$ is the response variable (class label or continuous value),\n", + "- $C > 0$ is the regularization strength (larger $C$ = less regularization),\n", + "- $\\ell_1\\text{ratio} \\in [0, 1)$ is the mixing parameter: $\\ell_1\\text{ratio} = 0$ gives Ridge,\n", + "- $\\omega_j > 0$ is a data-dependent $\\ell_1$-penalty weight applied to $\\beta_j$,\n", + "- $\\mathbf{A}\\beta + \\mathbf{b} \\geq \\mathbf{0}$ represents optional linear constraints on $\\beta$.\n", + "\n", + "In this example, we use squared loss and let $\\omega_j = {|\\beta^{\\text{OLS}}_j|^{-1}}$" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "bd45500e", + "metadata": {}, + "outputs": [], + "source": [ + "## packages\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from sklearn.datasets import make_regression\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.preprocessing import StandardScaler\n", + "from sklearn.preprocessing import StandardScaler\n", + "from sklearn.linear_model import LinearRegression\n", + "\n", + "from rehline import plqERM_ElasticNet" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "2142dda8", + "metadata": {}, + "outputs": [], + "source": [ + "## data simulation\n", + "n, d = 1000, 12\n", + "seed = 42\n", + "C = 0.001\n", + "l1_ratio = 0.5\n", + "\n", + "X, y, beta_true = make_regression(\n", + " n_samples=n,\n", + " n_features=d,\n", + " noise=0.1,\n", + " random_state=seed,\n", + " n_informative=4,\n", + " coef=True\n", + ")\n", + "X_train, X_test, y_train, y_test = train_test_split(\n", + " X,\n", + " y,\n", + " test_size=0.2,\n", + " random_state=seed\n", + ")\n", + "scaler = StandardScaler()\n", + "X_train_standered = scaler.fit_transform(X_train)\n", + "X_test_standered = scaler.transform(X_test)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "6fab08c9", + "metadata": {}, + "outputs": [], + "source": [ + "## fit the OLS solution\n", + "ols = LinearRegression(fit_intercept=False)\n", + "ols.fit(X_train, y_train)\n", + "beta_ols = ols.coef_\n", + "omega = 1 / np.abs(beta_ols)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "aa324935", + "metadata": {}, + "outputs": [], + "source": [ + "## apply the adaptive LASSO weights\n", + "clf_ENAL = plqERM_ElasticNet(\n", + " loss={\"name\":\"MSE\"},\n", + " C=C,\n", + " l1_ratio=l1_ratio,\n", + " omega=omega,\n", + " max_iter=30000,\n", + " shrink=seed,\n", + " verbose=1\n", + ")\n", + "clf_ENAL.fit(X_train, y_train)\n", + "beta_ENAL = clf_ENAL.coef_.flatten()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "600529d7", + "metadata": {}, + "outputs": [], + "source": [ + "## fit an elastic net model for comparison\n", + "clf_EN = plqERM_ElasticNet(\n", + " loss={\"name\":\"MSE\"},\n", + " C=C,\n", + " l1_ratio=l1_ratio,\n", + " max_iter=30000,\n", + " shrink=seed,\n", + " verbose=1\n", + ")\n", + "clf_EN.fit(X_train, y_train)\n", + "beta_EN = clf_EN.coef_.flatten()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "f257c0aa", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAAIhCAYAAABwnkrAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAsQRJREFUeJzs3Xl4U1XiPvD3JmmbdKGlQFcKLaXIYpEibqBSF0AEBRVQwBEUUAeXQWcENxSVAWG+8sNdB9kcNxQVHRlBXEAUERCBDgqU0pZSulC6b2mTe35/dBopbSEpTW5uzvt5njzc3qTJOXlPbjk5556rCCEEiIiIiIiIyMGgdQGIiIiIiIi8DTtKREREREREp2FHiYiIiIiI6DTsKBEREREREZ2GHSUiIiIiIqLTsKNERERERER0GnaUiIiIiIiITsOOEhERERER0WnYUSIiIiIiIjoNO0pEbfTSSy9BURScf/75Lv3eqlWroCgKsrKy3FOw/1mwYAHWrVvXbP/mzZuhKAo2b97s1tc/XVZWFhRFafU2b948x2OnTp2K+Ph4t5XFE+/NqfX94IMPmt0/b948KIqCoqIil59727ZtmDdvHkpLS8+5nO1h6tSpZ8y2UeN7smrVKreU40zvS2pqKlJTU9vldRrr269fP9jt9mb3K4qC+++/v03P3VrbbIsjR47g/vvvR69evWCxWBAYGIh+/frhySefRG5ubru8Rmu++eYbDBo0CEFBQVAUxVGnNWvWoF+/frBYLFAUBXv27HF8Flzl7uMEABw/fhzz5s3Dnj17nHp84zGktZsrbb+6uhrz5s1r8Xjkqb8jLTnTe9LWLIm8lUnrAhDp1YoVKwAA+/fvx88//4xLLrlE4xI1tWDBAowbNw5jx45tsn/gwIH46aef0LdvX03K9cADD2DSpEnN9nft2tVjZfD0e/PEE0/glltugZ+fX7s837Zt2/DMM89g6tSpCAsLa5fnPFcWiwXffvutpmU40/vy2muvtfvr/fbbb1i1ahWmTZvWbs/ZWtt01RdffIHbbrsNnTt3xv3334+UlBQoioK0tDSsWLEC69evx6+//to+hT6NEAITJkxAr1698PnnnyMoKAjnnXceTpw4gT/96U+47rrr8NprryEgIAC9evXC9OnTcd1117n8OnPnzsVf/vIXN9TgD8ePH8czzzyD+Ph4DBgwwOnfW7BgAa666qpm+xMTE51+jurqajzzzDMA0KyTP2rUKPz000+Ijo52+vnay5nek7ZmSeSt2FEiaoNdu3Zh7969GDVqFNavX4/ly5d7XUepNR06dMCll16q2et369ZN09c/E3e8NyNHjsSXX36JN954Aw888EC7Prc3MRgMXpsrgHbv/AYFBWHgwIF4+umnMWnSJFgslnZ9/nORmZmJ2267Db169cJ3332H0NBQx31XX301HnzwQXz66adue/3jx4+juLgYN910E6655hrH/h9//BH19fW4/fbbMXToUMf+wMDANn1R4kqnw9OSkpLc+nno0qULunTp4rbnb6uuXbt69EsvInfj1DuiNli+fDkA4Pnnn8fgwYPxwQcfoLq6utnjtm/fjiFDhsBsNiMmJgaPPfYY6uvrmz1uzZo1GD58OKKjo2GxWNCnTx88+uijqKqqavK4qVOnIjg4GPv378c111yDoKAgdOnSBffff3+T11cUBVVVVVi9erVjykfjN5KnTy9bunQpFEXB4cOHm5Vrzpw58Pf3bzI97Ouvv8Y111yDDh06IDAwEEOGDME333zj8nvoqldffRVXXnklIiIiEBQUhOTkZCxevLjZ+/nrr79i9OjRiIiIQEBAAGJiYjBq1CgcO3YMgGvvTaOff/4ZN9xwAzp16gSz2YzExETMmjXLqXJfffXVGDFiBJ577jlUVFSc9fFne3/nzZuHRx55BACQkJDgqENr0wVdyfds7507HD58GHfeeSeSkpIQGBiI2NhY3HDDDUhLS2vyOFVVMX/+fJx33nmwWCwICwtD//798eKLLwI4+/vS0tQ7q9WKZ599Fn369IHZbEanTp1w1VVXYdu2bU6VfdGiRcjNzXWU4UzKy8vxt7/9DQkJCfD390dsbCxmzZrV5DN+prbpiiVLlqCqqgqvvfZak07Sqa9z8803N9m3YsUKXHDBBTCbzQgPD8dNN92E33//vdnv7tq1CzfeeCPCw8NhNpuRkpKCDz/80HH/vHnzHP9RnjNnDhRFQXx8PKZOnYrLL78cAHDrrbc2qVtr07Xee+89XHbZZQgODkZwcDAGDBjgOPYCLU+9E0Lgtddew4ABA2CxWNCxY0eMGzcOR44cafK41NRUnH/++di5cyeuuOIKBAYGokePHnj++eehqiqAhuPBRRddBAC48847W5wmfC6+/fZbpKamolOnTrBYLOjWrRtuueUWVFdXIysry9EReuaZZxyvPXXqVAAtT71rrNNPP/2EwYMHw2KxID4+HitXrgQArF+/HgMHDkRgYCCSk5OxYcOGJuVx5rN4tvekpSxVVcXixYvRu3dvBAQEICIiAnfccUez44ozmRB5GjtKRC6qqanB+++/j4suugjnn38+7rrrLlRUVOCjjz5q8rjffvsN11xzDUpLS7Fq1Sq88cYb+PXXXzF//vxmz5meno7rr78ey5cvx4YNGzBr1ix8+OGHuOGGG5o9tr6+Htdffz2uueYarFu3Dvfffz/efPNN3HrrrY7H/PTTT7BYLLj++uvx008/4aeffmp16tHtt98Of3//ZnPn7XY73nnnHdxwww3o3LkzAOCdd97B8OHD0aFDB6xevRoffvghwsPDMWLECKc7S6qqwmazNbudTUZGBiZNmoR//etf+OKLLzBt2jT84x//wD333ON4TFVVFYYNG4aCggK8+uqr2LRpE5YuXYpu3bo5OimuvDcAsHHjRlxxxRU4evQolixZgi+//BJPPvkkCgoKnKov0PAf6qKiIvzjH/844+OceX+nT5/uGJn65JNPHHUYOHBgi8/pbL7OvHdn01KuZ/sPzvHjx9GpUyc8//zz2LBhA1599VWYTCZccsklOHjwoONxixcvxrx58zBx4kSsX78ea9aswbRp0xznI7n6vthsNowcORLPPfccRo8ejU8//RSrVq3C4MGDcfToUafqe9lll+Gmm27CokWLUFxc3OrjqqurMXToUKxevRoPPvggvvzyS8yZMwerVq3CjTfeCCEEgLO3zfj4eKfOyfnqq68QGRnp9IjGwoULMW3aNPTr1w+ffPIJXnzxRezbtw+XXXYZ0tPTHY/77rvvMGTIEJSWluKNN97AZ599hgEDBuDWW291tK/p06fjk08+AdAwzfann37Cp59+irlz5+LVV18F0DAt7Wyfu6eeegqTJ09GTEwMVq1ahU8//RRTpkxBdnb2Getyzz33YNasWbj22muxbt06vPbaa9i/fz8GDx7c7DObn5+PyZMn4/bbb8fnn3+OkSNH4rHHHsM777wDoGEqbmMn48knn3RkMn369LO+p2c7zmVlZWHUqFHw9/fHihUrsGHDBjz//PMICgpCXV0doqOjHR2ZadOmOV577ty5Z3zd/Px83HnnnZg+fTo+++wzJCcn46677sKzzz6Lxx57DLNnz8bHH3+M4OBgjB07FsePH3f8rjOfxba8J3/+858xZ84cDBs2DJ9//jmee+45bNiwAYMHD252jubZMiHyOEFELnn77bcFAPHGG28IIYSoqKgQwcHB4oorrmjyuFtvvVVYLBaRn5/v2Gez2UTv3r0FAJGZmdni86uqKurr68WWLVsEALF3717HfVOmTBEAxIsvvtjkd/7+978LAOKHH35w7AsKChJTpkxp9vzfffedACC+++47x76bb75ZdO3aVdjtdse+//znPwKA+Pe//y2EEKKqqkqEh4eLG264ocnz2e12ccEFF4iLL764xfo0yszMFABavW3durVJPbt3797qc9ntdlFfXy/efvttYTQaRXFxsRBCiF27dgkAYt26dWcsiyvvTWJiokhMTBQ1NTVnfM7W6vuPf/xDCCHE5MmTRVBQkMjLyxNCCPH0008LAOLEiRNCCNfe33/84x9nbEOncyZfZ9+7ljS2y5Zu11xzTbP3ZOXKla0+l81mE3V1dSIpKUk89NBDjv2jR48WAwYMOGM5zvS+DB06VAwdOtTxc+PneNmyZU7Xs9GUKVNEUFCQEEKIAwcOCKPRKP7617867gcg7rvvPsfPCxcuFAaDQezcubPJ86xdu1YAEP/5z38c+1prm0L80RbPxmw2i0svvdSpupSUlAiLxSKuv/76JvuPHj0qAgICxKRJkxz7evfuLVJSUkR9fX2Tx44ePVpER0c72tfpbb9R4+fro48+arK/8bPQ6MiRI8JoNIrJkyefseynHyd++uknAUC88MILTR6Xk5MjLBaLmD17tmPf0KFDBQDx888/N3ls3759xYgRIxw/79y586xttqU6tnbLyckRQvyR/Z49e1p9rhMnTggA4umnn25238qVK5u19cY67dq1y7Hv5MmTwmg0CovFInJzcx379+zZIwCIl156qdXXb+2zeKb35PQsf//9dwFAzJw5s8njfv75ZwFAPP74483Kf7ZMiDyJI0pELlq+fDksFgtuu+02AEBwcDDGjx+PrVu3Nvv29ZprrkFkZKRjn9FobDLy0+jIkSOYNGkSoqKiYDQa4efn55jD39L0l8mTJzf5uXFxhO+++65Ndbrzzjtx7NgxfP311459K1euRFRUFEaOHAmg4UT54uJiTJkypdmIwXXXXYedO3c2myrYkr/85S/YuXNns9vZTpT+9ddfceONN6JTp06O9+iOO+6A3W7HoUOHAAA9e/ZEx44dMWfOHLzxxhv47bff2vR+NDp06BAyMjIwbdo0mM3mc3qu+fPno76+3nFy9una6/1tiTP5nut7Z7FYWsz1bIso2Gw2LFiwAH379oW/vz9MJhP8/f2Rnp7epO1ffPHF2Lt3L2bOnImNGzeivLzcpfKd7ssvv4TZbMZdd911Ts9z3nnnYdq0aXjllVdaHYn64osvcP7552PAgAFNsh0xYoRLqywePny4xSmU5+Knn35CTU2NY0pXo7i4OFx99dWOkczDhw/jwIEDjmPPqfW4/vrrkZeX12QE8Fxs2rQJdrsd9913n0u/98UXX0BRFNx+++1NyhcVFYULLrig2fscFRWFiy++uMm+/v37n3XUyhmLFi1q8fPQ+PdgwIAB8Pf3x913343Vq1c3mxrYVtHR0bjwwgsdP4eHhyMiIgIDBgxATEyMY3+fPn0AoEldnf0suqLxb9Lp7eviiy9Gnz59ms1EcGcmRG3BjhKRCw4fPozvv/8eo0aNghACpaWlKC0txbhx4wD8sRIeAJw8eRJRUVHNnuP0fZWVlbjiiivw888/Y/78+di8eTN27tzpmL5SU1PT5PEmkwmdOnVq8TlPnjzZpnqNHDkS0dHRjikVJSUl+Pzzz3HHHXfAaDQCgGPayrhx4+Dn59fktmjRIgghzjj9qFHXrl0xaNCgZrfg4OBWf+fo0aO44oorHOeDbN26FTt37nRM5Wl8j0JDQ7FlyxYMGDAAjz/+OPr164eYmBg8/fTTLZ4bdjYnTpxwlPlcxcfHY+bMmXjrrbeadKgbtdf72xJn8j3X985gMLSYa69evc74ew8//DDmzp2LsWPH4t///jd+/vln7Ny5ExdccEGTtv/YY4/h//7v/7B9+3aMHDkSnTp1wjXXXINdu3a16T05ceIEYmJiYDCc+5/BefPmwWg0tjotqqCgAPv27WuWa0hICIQQbVoi/ky6deuGzMxMpx7beMxoafW0mJgYx/2N7fNvf/tbs3rMnDkTANqtHm393BUUFEAIgcjIyGZl3L59e7PynX4cBYCAgIBmx9y26NGjR4ufh8aVLxMTE/H1118jIiIC9913HxITE5GYmOjU+W5nEh4e3myfv79/s/3+/v4AgNraWsc+Zz+LrnC2fTVyZyZEbcFV74hcsGLFCgghsHbtWqxdu7bZ/atXr8b8+fNhNBrRqVMn5OfnN3vM6fu+/fZbHD9+HJs3b26yElRr18ix2Ww4efJkkz8ojc/Z0h8ZZxiNRvzpT3/CSy+9hNLSUrz33nuwWq248847HY9pPE/p5ZdfbvXch1NHz9rTunXrUFVVhU8++QTdu3d37G/pOh7Jycn44IMPIITAvn37sGrVKjz77LOwWCx49NFHXXrdxpOp22sxgyeffBIrVqxwdERO5c7315l8gfZ975z1zjvv4I477sCCBQua7C8qKmqyxLfJZMLDDz+Mhx9+GKWlpfj666/x+OOPY8SIEcjJyUFgYKBLr9ulSxf88MMPUFX1nDtL0dHRmDVrFp5//nn89a9/bXZ/586dYbFYmnyRcvr97WnEiBF4+eWXsX379rOep9R4zMjLy2t23/Hjxx1la/z3sccea7YQRKPzzjvvXIrtcOrnLi4uzunf69y5MxRFwdatWxEQENDs/pb2aemKK67AFVdcAbvdjl27duHll1/GrFmzEBkZ6Zix4EnOfhZdcWr7Or3je2r7IvJWHFEicpLdbsfq1auRmJiI7777rtntr3/9K/Ly8vDll18CAK666ip88803TU4gttvtWLNmTZPnbVwh6PQ/4m+++WarZXn33Xeb/Pzee+8BaHqtDVe/hbvzzjtRW1uL999/H6tWrcJll12G3r17O+4fMmQIwsLC8Ntvv7X4TemgQYMc31K2t5beIyEEli1bdsbfueCCC/D//t//Q1hYGHbv3u24z9n3plevXkhMTMSKFStgtVrPoQYNOnXqhDlz5mDt2rXYsWNHk/tceX8b34f2zPdUZ3rv2puiKM3a/vr16894QdSwsDCMGzcO9913H4qLix0rf7nyvowcORK1tbXtdvHbOXPmIDw8vMUO5ejRo5GRkYFOnTq1mOupCzS0x7fnDz30EIKCgjBz5kyUlZU1u18I4Vge/LLLLoPFYml2svyxY8fw7bffOpb3Pu+885CUlIS9e/e22j5DQkLOqdyNhg8fDqPRiNdff92l3xs9ejSEEMjNzW2xfMnJyS6XpS2fNVcZjUZccskljhHyxs+bJ177VM5+Fl0p19VXXw0AzdrXzp078fvvvzdZPp7IG3FEichJX375JY4fP45Fixa1uGTv+eefj1deeQXLly/H6NGj8eSTT+Lzzz/H1VdfjaeeegqBgYF49dVXm51nMnjwYHTs2BH33nsvnn76afj5+eHdd9/F3r17WyyHv78/XnjhBVRWVuKiiy7Ctm3bMH/+fIwcOdKx/C7QMDqwefNm/Pvf/0Z0dDRCQkLO+I1v7969cdlll2HhwoXIycnBP//5zyb3BwcH4+WXX8aUKVNQXFyMcePGISIiAidOnMDevXtx4sQJp/5jc/ToUWzfvr3Z/i5durR6XZRhw4bB398fEydOxOzZs1FbW4vXX38dJSUlTR73xRdf4LXXXsPYsWPRo0cPCCHwySefoLS0FMOGDWvTe/Pqq6/ihhtuwKWXXoqHHnoI3bp1w9GjR7Fx48ZmHVZnzJo1C6+++qqjQ93Ilfe38T98L774IqZMmQI/Pz+cd955Z/yP6tnydfa9a42qqi3mCgApKSmtfps/evRorFq1Cr1790b//v3xyy+/4B//+Eezb59vuOEGnH/++Rg0aBC6dOmC7OxsLF26FN27d0dSUpLL78vEiROxcuVK3HvvvTh48CCuuuoqqKqKn3/+GX369HH5G/0OHTrgiSeewEMPPdTsvlmzZuHjjz/GlVdeiYceegj9+/eHqqo4evQovvrqK/z1r391XIftTG2zZ8+eAHDW85QSEhLwwQcf4NZbb8WAAQMcF5wFGlbjbBwZv+mmmxAWFoa5c+fi8ccfxx133IGJEyfi5MmTeOaZZ2A2m/H00087nvfNN9/EyJEjMWLECEydOhWxsbEoLi7G77//jt27dzdb+bOt4uPj8fjjj+O5555DTU0NJk6ciNDQUPz2228oKipq9Ty/IUOG4O6778add96JXbt24corr0RQUBDy8vLwww8/IDk5GX/+859dKktiYiIsFgveffdd9OnTB8HBwYiJiWlyvk9L0tPTW/w8NF5n6I033sC3336LUaNGoVu3bqitrXWMOF577bUAgJCQEHTv3h2fffYZrrnmGoSHh6Nz585OrXzYFs5+Fl15T8477zzcfffdePnll2EwGDBy5EhkZWVh7ty5iIuLa/HzQuRVtFhBgkiPxo4dK/z9/UVhYWGrj7ntttuEyWRyrHT3448/iksvvVQEBASIqKgo8cgjj4h//vOfzVYr2rZtm7jssstEYGCg6NKli5g+fbrYvXt3s5WFGlfb2rdvn0hNTRUWi0WEh4eLP//5z6KysrJJWfbs2SOGDBkiAgMDBQDHil8trezWqLFsFotFlJWVtVjHLVu2iFGjRonw8HDh5+cnYmNjxahRo5qtZHW6s616d+oKVy2tevfvf/9bXHDBBcJsNovY2FjxyCOPiC+//LJJXQ4cOCAmTpwoEhMThcViEaGhoeLiiy8Wq1atOqf35qeffhIjR44UoaGhIiAgQCQmJjZZBepM9T195S8h/nifccqqd42cfX8fe+wxERMTIwwGQ6t5tva6LeXr7HvXkjOtegdApKenN3lPTm3TJSUlYtq0aSIiIkIEBgaKyy+/XGzdurXZKnUvvPCCGDx4sOjcubPw9/cX3bp1E9OmTRNZWVlOvS+nP58QQtTU1IinnnpKJCUlCX9/f9GpUydx9dVXi23btp21vo2r3p3KarWKhISEZqveCSFEZWWlePLJJ8V5550n/P39RWhoqEhOThYPPfRQk5UxW2ubQgjRvXv3M64GebqMjAwxc+ZM0bNnTxEQECAsFovo27evePjhh5utDPjWW2+J/v37O8o2ZswYsX///mbPuXfvXjFhwgQREREh/Pz8RFRUlLj66qsdq4AKce6r3jV6++23xUUXXSTMZrMIDg4WKSkpzY6HLb0fK1asEJdccokICgoSFotFJCYmijvuuKPJanBDhw4V/fr1a/a7LT3n+++/L3r37i38/PxaXYXu9Dq2dnviiSeEEA3HlJtuukl0795dBAQEiE6dOomhQ4eKzz//vMnzff311yIlJUUEBAQIAI4VEVtb9a6lOnXv3l2MGjWq2f7T26mzn8UzvSctZWm328WiRYtEr169hJ+fn+jcubO4/fbbHSsAnq38Z1sFlcidFCH+dwEHIvJ6U6dOxdq1a1FZWal1UYiIiIh8Gs9RIiIiIiIiOg07SkRERERERKfh1DsiIiIiIqLTcESJiIiIiIjoNOwoERERERERnYYdJSIiIiIiotP4/AVnVVXF8ePHERISAkVRtC4OERERERFpRAiBiooKxMTEwGA485iRz3eUjh8/jri4OK2LQUREREREXiInJwddu3Y942M07yjl5uZizpw5+PLLL1FTU4NevXph+fLluPDCCwE09PqeeeYZ/POf/0RJSQkuueQSvPrqq+jXr59Tzx8SEgKg4c3o0KGD2+pBRERERETerby8HHFxcY4+wplo2lEqKSnBkCFDcNVVV+HLL79EREQEMjIyEBYW5njM4sWLsWTJEqxatQq9evXC/PnzMWzYMBw8eNCpCjZOt+vQoYPmHSW73Y6MjAwkJibCaDRqWhbyPOYvN+ZPbANyY/5yY/7ex5lTcjS9jtKjjz6KH3/8EVu3bm3xfiEEYmJiMGvWLMyZMwcAYLVaERkZiUWLFuGee+4562uUl5cjNDQUZWVlmneUhBCoq6uDv78/z5eSEPOXG/MntgG5MX+5MX/v4UrfQNNV7z7//HMMGjQI48ePR0REBFJSUrBs2TLH/ZmZmcjPz8fw4cMd+wICAjB06FBs27atxee0Wq0oLy9vcgMaFnVo/NeZ7cb+Y2vbdrv9jNtCiGbbiqLA39/f8TqN+xu3T93f0razZfdknVqqB+vUcnmFEPDz84OiKD5TJ1/MyV11AhqOX42P8YU6+WJO7qyTEAIBAQHnXD9vqpMv5uSuOimKAj8/PzTyhTr5Yk7uqpOqqggICHA8vy/USc85OUvTjtKRI0fw+uuvIykpCRs3bsS9996LBx98EG+//TYAID8/HwAQGRnZ5PciIyMd951u4cKFCA0NddwaF3IoLCwEABQVFaGoqAgAUFBQgOLiYgBAXl4eSktLATScN1VWVgag4dymiooKAEBWVhaqqqoANHTiamtrAQAZGRmoq6sDAKSnp8Nms0FVVaSnp0NVVdhsNqSnp8Nut2P//v04fPgwAKC2thaZmZkAgKqqKmRlZQEAKioqkJOTAwAoKytDbm4uAKC0tBR5eXkAgOLiYhQUFGheJwCoq6tDRkYG63SWOhUUFODXX3+F3W73mTr5Yk7uqlNNTQ0OHDiAgwcP+kydfDEnd9YpJycHBw4cQHFxsc/UyRdzcled7HY7fv75Z0f9fKFOvpiTu+p05MgRHDhwABUVFT5TJz3n5CxNp975+/tj0KBBTUaHHnzwQezcuRM//fQTtm3bhiFDhuD48eOIjo52PGbGjBnIycnBhg0bmj2n1WqF1Wp1/Nx4wlZJSQnCwsIcvVODwXDGbUVRHN/8t7Rtt9thMBha3QYaerqnb9fX18NgMMBkMjl62kaj0dHDNRgMrW47W3ZP1qmx7K1ts05/bNvtdthsNvj7+zu+zdB7nXwxJ3fV6dR/G+en671OvpiTO+vU+A1p42v4Qp18MSd31Qlo+I+mn5+f42+CJ+ukKArq6uocz8+cPFunxtGMxuO/L9TJW3MCGvoYp37OTq1TeXk5wsLCnJp6p+liDtHR0ejbt2+TfX369MHHH38MAIiKigLQMLJ0akepsLCw2ShTo4CAAMfQ5qka36TGf891+9QT8ZzdbvyAND7Pqf9hagz6TNvtVfb2rFNL9WCdWi+XyWTyuTq1x7YMdWr8Y3HqgVzvdXJ2m3WCo2PU+PrurB9z8s46CSFgMpnO+Pl3V53q6+uRl5eH6upqkHYaO6zkfoGBgYiOjoa/v79j36mfFWdp2lEaMmQIDh482GTfoUOH0L17dwBAQkICoqKisGnTJqSkpABo+DZmy5YtWLRokcfLe65UtWHYLykpiSueSIj5y435E9uA3LTKX1VVZGZmwmg0IiYmhosJaEQIAavVioCAAL7/biREw6IZJ06cQGZmJpKSkpp8ceEqTTtKDz30EAYPHowFCxZgwoQJ2LFjB/75z3/in//8J4CGHt+sWbOwYMECJCUlISkpCQsWLEBgYCAmTZqkZdHbxGAwnHNgpF/MX27Mn9gG5KZV/nV1dVBVFXFxcQgMDPToa9MfhBAwm80AXBvRINdZLBb4+fkhOzsbdXV1jve9LTTtKF100UX49NNP8dhjj+HZZ59FQkICli5dismTJzseM3v2bNTU1GDmzJmOC85+9dVXTl1DyRudOueT5MP85cb8iW1Ablrmz3anPU6985z2au+aLubgCd50HSW73c5pFxJj/nJj/sQ2IDet8m9cnSwhIeGcvlmncyOEQG1tLcxmMztLHnCmdu9K30DTESXZGI1G9O7dW+tikEaYv9yYP7ENyI35y01RFFgsFq2LQS5iR8mDGk8w44mUcmL+cmP+xDYgN73nX1BQgM2bN6OiogIhISFITU1tdQViau7UZdr1mL+sOGHVg1RVRXZ2tmNteJKLzPkXFBRgzZo1eOutt7BmzRrHheZkInP+1IBtQG56zT8tLQ0TJ05E165dcdttt2HGjBm47bbb0LVrV0ycOBFpaWnt/pqNnYnWblOnTm331/SExougkn7wHCUicpu0tDQsWLAAa9eubXIlbJPJhHHjxuHxxx9HcnKyhiUkIvJt53KO0saNGzF27FjYbLYmx/BGJpMJJpMJ69atw4gRI9qryMjPz3dsr1mzBk899VSTy8lYLBaEhoY6fq6vr4efn1+7vT7pX3udo8QRJQ8SQqCmpgY+3jelVsiW/8aNG3HxxRc36yQBgM1mw9q1a3HxxRdj48aNGpXQs2TLn5pjG5Cb3vJPS0vD2LFjYbVaW+wkAQ3HcqvVirFjx7bryFJUVJTjFhoaCkVRHD/X1tYiLCwMH374IVJTU2E2m/HOO+9g3rx5GDBgQJPnWbp0KeLj45vsW7lyJfr06QOz2YzevXvjtddea7dyn4kQAqqq6iZ/asCOkgepqorc3FzdDbtT+5Apfy3/wHormfKnlrENyE1v+S9YsAA2m+2s/7EXQsBms2HhwoUeKlmDOXPm4MEHH8Tvv//u9GjWsmXL8MQTT+Dvf/87fv/9dyxYsABz587F6tWr3VzaBpx6pz/sKHmQ0WhEz549uSyspGTK39v/wGpBpvypZWwDctNT/gUFBS3OBmiNzWbDRx99hMLCQjeX7A+zZs3CzTffjISEBMTExDj1O8899xxeeOEFx+/dfPPNeOihh/Dmm2+6ubQN511xaXD9YUfJg4QQqKys5LCrpGTJXw9/YLUgS/7UOrYBuekp/82bNzt9DG9ks9mwefNm9xSoBYMGDXLp8SdOnEBOTg6mTZuG4OBgx23+/PnIyMhwUyn/IISA3W7XRf70B3aUPEgIgcLCQn5IJCVL/nr4A6sFWfKn1rENyE1P+VdUVLTp98rLy9u5JK0LCgpq8rPBYGj23tbX1zu2G6c8Llu2DHv27HHc/vvf/2L79u3uLzDg8t9G0h6vo+RBBoMBPXr00LoYpBFZ8tfDH1gtyJI/tY5tQG56yj8kJKRNv6fl6sJdunRBfn6+41pFALBnzx7H/ZGRkYiNjcWRI0cwefJkj5dPURQEBAR4/HXp3LCj5EFCCMeF2jhHVT6y5K/HP7CeIEv+1Dq2AbnpKf/U1FSYTCaXRkBMJhNSU1PdV6izSE1NxYkTJ7B48WKMGzcOGzZswJdfftnkb8u8efPw4IMPokOHDhg5ciSsVit27dqFkpISPPzww24tX+OqdwaDwevzpz9w6p0HCSFQUlKii2F3an+y5N/4B9YVWv+B9QRZ8qfWsQ3ITU/5R0ZGYty4cU4fy00mE8aPH4+IiAg3l6x1ffr0wWuvvYZXX30VF1xwAXbs2IG//e1vTR4zffp0vPXWW1i1ahWSk5MxdOhQrFq1CgkJCR4pI6fe6Q8vOEtE7W7ixIlOL+jQ+Af2vffe80DJiIjk0tYLzqalpeHiiy+G1Wo9Y+eucUrZjh07eAFx8hq84KwOCSFQWlqqi2+TqP3JlP/jjz8Ok8l01ukFiqLAZDLhscce81DJtCNT/tQytgG56S3/5ORkrFu3DgEBAa2OLJlMJgQEBGDdunXsJJ1F4+Uw9JI/NWBHyYMa5yfzQyInmfLnH9jmZMqfWsY2IDc95j9ixAjs2LED48ePb3Ysb5wNsGPHDqcv+Co7vVxsmP7AqXdE5DZpaWlYuHAhPvrooybT8Br/wD722GNSdJKIiLTS1ql3pyssLMTmzZtRXl6ODh06IDU1VdNzkojOpL2m3rGj5EGqqqK0tBRhYWEwGDiYJxuZ8+cfWLnzpwZsA3LTKv/26ijRuWm84KzRaOSqdx7QXh0lLg/uYTU1NQgLC9O6GKQRWfOPiIjAhAkTtC6G5mTNn/7ANiA35i83VVVhNBq1Lga5gB0lDzIYDIiNjdW6GKQR5i835k9sA3Jj/nJTFAX+/v5aF4NcxLF/D1JVFUVFRTyZT1LMX27Mn9gG5Mb85SaEQH19va4W8yB2lDyuvr5e6yKQhpi/3Jg/sQ3IjfnLjZ0k/eHUOw8yGAyIjo7WuhikEeYvN+ZPbANyY/5y49Q7feKIkgepqorCwkIOu0uK+cuN+RPbgNz0nn9pKbB1K/DVVw3/lpZqXaKGzse6devc+hrx8fFYunTpOT8Pp97pEztKRERERNSirCxg8WJgzBjg3nuBBx5o+HfMmIb9WVnued2pU6dCUZRmt+uuu84tr7dq1aoWVyTcuXMn7r777jY/b3x8PBRFwfbt25vsnzVrFlJTU51+nqysLCiKgj179rS5LOQ6Tr3zIIPBIN21Y+gPzF9uzJ/YBuSmx/x37wbmzAGOHgU6dADi4gCTCbDZgJMngeXLgU2bgEWLgIED2//1r7vuOqxcubLJvoCAgPZ/oTPo0qXLOT+H2WzGo48+ii1btrRDiciTOKLkQaqqIi8vT7fD7nRumL/cmD+xDchNb/lnZTV0ko4dAxITgagowM8PUJSGf6OiGvYfO9bwOHeMLAUEBCAqKqrJrWPHjq0+fs6cOejVqxcCAwPRo0cPzJ07t8kCGnv37sVVV12FkJAQdOjQARdeeCF27dqFzZs3484770RZWZlj5GrevHkAmk+9Ky0txd13343IyEiYzWacf/75+OKLL85Yj3vuuQfbt2/HZ599dsapdytXrkSfPn1gNpvRu3dvvPbaa477EhISAAApKSlQFMWl0ShqO44oeZifn5/WRSANMX+5MX9iG5CbnvL/8MOGkaTERKC1a6QajUBCApCRAXz0EfDII54t4+lCQkKwatUqxMTEIC0tDTNmzEBISAhmz54NAJg8eTJSUlLw+uuvw2g0Ys+ePfDz88PgwYOxdOlSPPXUUzh48CAAIDg4uNnzq6qKkSNHoqKiAu+88w4SExPx22+/nfUisvHx8bjnnnswd+5cjB49usXHL1u2DE8//TReeeUVpKSk4Ndff8WMGTMQFBSEKVOmYMeOHbj44ovx9ddfo1+/flwYwkPYUfIgg8GAzp07a10M0gjzlxvzJ7YBuekp/9JSYP36hul2Z+kDwGhseNz69cDddwOhoe1Xji+++KJZh2XOnDmYO3dui49/8sknHdvx8fH461//ijVr1jg6SkePHsUjjzyC3r17AwCSkpIcjw8NDYWiKIiKimq1PF9//TV27NiB33//Hb169QIA9OjRw6m6zJ07F6tWrcJ7772HP/3pT83uf+655/DCCy/g5ptvBtAwgvTbb7/hzTffxJQpUxxTADt16nTGMlL7YkfJgxqH3aOjo2EwcNajbJi/3Jg/sQ3ITU/5p6UBRUUN5yQ5o1MnICen4fcuv7z9ynHVVVfh9ddfb7IvPDy81cevXbsWS5cuxeHDh1FZWQmbzYYOHTo47n/44Ycxffp0/Otf/8K1116L8ePHIzEx0eny7NmzB127dnV0klzRuXNnPPTQQ3jqqadw6623NrnvxIkTyMnJwbRp0zBjxgzHfpvNhtD27HmSy7z7k+qDLBaL1kUgDTF/uTF/YhuQm17yr6lpWLDB5OTX6SYTYLcD1dXtW46goCD07Nmzya21jtL27dtx2223YeTIkfjiiy/w66+/4oknnkBdXZ3jMfPmzcP+/fsxatQofPvtt+jbty8+/fRTp8tzrvk99NBDqKmpaXLuEQDHeWvLli3Dnj17HLf//ve/zVbLI8/iiJIHGQyGM34TQr6N+cuN+RPbgNz0lL/F8sfqds6cVmWzNUzBCwx0f9la8+OPP6J79+544oknHPuys7ObPa5Xr17o1asXHnroIUycOBErV67ETTfdBH9/f9jt9jO+Rv/+/XHs2DEcOnTI5VElRVEQFhaGuXPnYt68ebjhhhsc90VGRiI2NhZHjhzB5MmTW/z9xnOSzlZGal8cUfIgVVWRk5OjmxVvqH0xf7kxf2IbkJue8k9OBjp3blgC3BknTwJdujT8XnuyWq3Iz89vcisqKmrxsT179sTRo0fxwQcfICMjAy+99FKT0aKamhrcf//92Lx5M7Kzs/Hjjz9i586d6NOnD4CGc5oqKyvxzTffoKioCNUtDI8NHToUV155JW655RZs2rQJmZmZ+PLLL7Fhw4az1kUIgbq6OsyYMQOhoaF4//33m9w/b948LFy4EC+++CIOHTqEtLQ0rFy5EkuWLAEAREREwGKxYMOGDSgoKEBZWZnT7yO1HTtKHqQoCkJCQqAoitZFIQ0wf7kxf2IbkJue8g8LA0aNAsrLG6bUnYnd3vC4UaPadyEHANiwYQOio6Ob3C5v5SSoMWPG4KGHHsL999+PAQMGYNu2bU0WfTAajTh58iTuuOMO9OrVCxMmTMDIkSPxzDPPAAAGDx6Me++9F7feeiu6dOmCxYsXt/g6H3/8MS666CJMnDgRffv2xezZs50e5TEYDPDz88Nzzz2H2traJvdNnz4db731FlatWoXk5GQMHToUq1atciwLbjKZ8NJLL+HNN99ETEwMxowZ49Rr0rlRxJkWdPcB5eXlCA0NRVlZWZMT+oiIiIh8XW1tLTIzM5GQkACz2ez072VlATNmNFwnKSGh5dXv7HYgMxPo2hVYtgyIj2+3YhOdkzO1e1f6BhxR8iBVVZGdna2LYXdqf8xfbsyf2Abkprf84+OBRYsaOkEZGUB+PlBfDwjR8G9+fsP+rl0bHsdO0pkJIWC1Ws94wVnyPlzMwYMURUHHjh11MexO7Y/5y435E9uA3PSY/8CBDSNFH33UcJ2knJyGUSSjseGcpAkTgPHj2UlylsnZZQTJazAxD1IUhdP/JMb85cb8iW1AbnrNPz4eeOSRhovJpqU1LAEeGNiwcAMv8eM8RVFgPNvVe8nrcOqdB6mqiiNHjuhm2J3aF/OXG/MntgG56T3/0NCGi8kOH97wLztJruHUO31iR8mDFEVBRESErobdqf0wf7kxf2IbkBvzJ0690x8m5kGKoiA4OFjrYpBGmL/cmD+xDciN+cuNU+/0iSNKHmS323H48GFeVVlSzF9uzJ/YBuTG/OUmhEBtbS2n3ukMO0oeZDAYEBsbC4OBb7uMmL/cmD+xDciN+ZO/v7/WRSAXceqdBymKAovFonUxSCPMX27Mn9gG5Kb7/EtLG5a9q6kBLJaGZe/CwrQulW4oisLz03SIX2t4kN1ux6FDhzjsLinmLzfmT2wDctNt/llZwOLFwJgxwL33Ag880PDvmDEN+7OytC5hE/PmzcOAAQM88lpTp07F2LFjnXqsO6fexcfHY+nSpe3+vKdKTU3FrFmz3Poa3ogdJQ8yGAzo3r07h90lxfzlxvyJbUBuusx/925gxgxg+fKGCyjFxQFJSQ3/Vlc37J8xo+FxbrJt2zYYjUZcd911bnuNs8nKyoKiKNizZ0+T/S+++CJWrVrl9PO0ZerdvHnzHKNRp9569+7t8nM5Y/PmzVAUBaWlpU32f/LJJ3juuefa/LypqalQFAUffPBBk/1Lly5FvItXLFYUBevWrWtzWVyho0+r/imKgoCAAA69Sor5y435E9uA3HSXf1YWMGcOcOwYkJgIREUBfn6AojT8GxXVsP/YsYbHuWlkacWKFXjggQfwww8/4OjRo255jbYKDQ1FmJPTDxVFgcFgaFP+/fr1Q15eXpPbDz/84PLznIvw8HCEhISc03OYzWY8+eSTqK+vb6dSuR87Sh5kt9tx4MAB/Q27U7tg/nJj/sQ2IDfd5f/hh8DRo0BCAtDastZGY8P9R48CH33U7kWoqqrChx9+iD//+c8YPXp0i6M3zz//PCIjIxESEoJp06ahtra2yf07d+7EsGHD0LlzZ4SGhmLo0KHYfdoImKIoeP311zFy5EhYLBYkJCTgo1Pqk5CQAABISUmBoihITU0F0HTq3ZtvvonY2NhmFxS+8cYbMWXKFAghUFNTg88//xwXXnghzGYzevTogWeeeQY2m+2M74PJZEJUVFSTW+fOnVt9/JIlS5CcnIygoCDExcVh5syZqKysdNyfnZ2NG264AR07dkRQUBD69euH//znP8jKysJVV10FAOjYsSMURcHUqVMBNJ96Z7VaMXv2bMTFxSEgIABJSUlYvnz5GesxceJElJWVYdmyZWd83L///e9W36PG0aebbroJiqK4PBrlKnaUPMhgMCAxMVFfw+7Ubpi/3Jg/sQ3ITVf5l5YC69cDHTq03klqZDQ2PG79eqCsrF2LsWbNGpx33nk477zzcPvtt2PlypVNzvH58MMP8fTTT+Pvf/87du3ahejoaLz22mtNnqOiogJTpkzB1q1bsX37diQlJeH6669HRUVFk8fNnTsXt9xyC/bu3Yvbb78dEydOxO+//w4A2LFjBwDg66+/Rl5eHj755JNmZR0/fjyKiorw3XffOfaVlJRg48aNmDx5MoCGaW1/+tOf8OCDD+K3337Dm2++iVWrVuHvf/97+7xh/2MwGPDSSy/hv//9L1avXo1vv/0Ws2fPdtx/3333wWq14vvvv0daWhoWLVqE4OBgxMXF4eOPPwYAHDx4EHl5eXjxxRdbfI077rgDH3zwAV566SX8/vvveOONN856nbAOHTrg8ccfx7PPPouqqqoWH7Nx40bcfvvtrb5HO3fuBACsXLkSeXl5jp/dRQefVt+iiwMkuQ3zlxvzJ7YBuekm/7Q0oKgI6NTJucd36gScONHwe+1o+fLluP322wEA1113HSorK/HNN9847l+6dCnuuusuTJ8+Heeddx7mz5+Pvn37NnmOq6++Grfffjv69OmDPn364M0330R1dTW2bNnS5HHjx4/H9OnT0atXLzz33HMYNGgQXn75ZQBAly5d/lfNToiKikJ4eHizsoaHh+O6667De++959j30UcfITw8HNdccw0AYOHChZgzZw6mTJmCHj16YNiwYXjuuefw5ptvnvF9SEtLQ3BwcJPb9OnTW338rFmzcNVVVyEhIQFXX301nnvuOXz44YeO+48ePYohQ4YgOTkZPXr0wOjRo3HllVfCaDQ66hYREYGoqCiEhoY2e/5Dhw7hww8/xIoVK3DTTTehR48euOaaa3DrrbeesR4AMHPmTJjNZixZsqTF+//+97/j0UcfbfU9aswiLCwMUVFRjp/dRSefWN+gqirS09ObDcuSHJi/3Jg/sQ3ITVf519QANhtgcvIqMiYTYLc3LPDQTg4ePIgdO3bgtttu+99LmHDrrbdixYoVjsf8/vvvuOyyy5r83uk/FxYW4t5770WvXr0QGhqK0NBQVFZWNjvfqaXnaRxRctbkyZPx8ccfw2q1AgDeffdd3HbbbTD+b1Tul19+wXPPPdekwzNjxgzk5eWh+gzv3XnnnYc9e/Y0uZ1pFOq7777DsGHDEBsbi5CQENxxxx04efKkYxTnwQcfxPz58zFkyBA8/fTT2Ldvn0v13LNnD4xGI4YOHerS7wFAQEAAnn32WfzjH/9AUVFRs/t/+eUXPPvssy6/R+7C6yh5kMFgQFJSkn6+UaJ2xfzlxvyJbUBuusrfYmno/NhsDQs3nI3N1jAFLzCw3YqwfPly2Gw2xMbGOvYJIeDn54eSkhJ07NjRqeeZOnUqTpw4gaVLl6J79+4ICAjAZZddhrq6urP+rqsLL9xwww1QVRXr16/HRRddhK1btzYZOVFVFfPmzcMtt9zS7HfNZnOrz+vv74+ePXs6VYbs7Gxcf/31uPfee/Hcc88hPDwcP/zwA6ZNm+ZYRGH69OkYMWIE1q9fj6+++goLFy7ECy+8gAceeMCp1zjX64Hdfvvt+L//+z/Mnz+/2TlGqqrimWeewc0339zs9870HrmLDj6tvkUX3ySR2zB/uTF/YhuQm27yT04GOncGTp507vEnTwJdujT8Xjuw2Wx4++238cILLzQZRdm7dy+6d++Od999FwDQp08fbN++vcnvnv7z1q1b8eCDD+L6669Hv379EBAQ0OJIRkvP07gEd+Oy3mdbiMNiseDmm2/Gu+++i/fffx+9evXChRde6Lh/4MCBOHjwIHr27Nns1l4d6F27dsFms+GFF17ApZdeil69euH48ePNHhcXF4d7770Xn3zyCf761786Flhwpq7JyclQVbXZ9EVnGQwGLFy4EK+//jqyTlst0Zn3yM/Pz2OLonBEyYNUVUVGRgaSkpIcw7AkD+YvN+ZPbANy01X+YWHAqFEN10nq0uXMCzrY7UB5OTBhAtDC+Sxt8cUXX6CkpATTpk1rdo7MuHHjsHz5ctx///34y1/+gilTpmDQoEG4/PLL8e6772L//v3o0aOH4/E9e/bEv/71LwwaNAjl5eV45JFHWhwR+eijj5o8z44dOxyruEVERMBisWDDhg3o2rUrzGZzi+fuAA3T72644Qbs37/fcX5Vozlz5uCWW25Bt27dMH78eBgMBuzbtw9paWmYP39+q++HzWZDfn5+k32KoiAyMrLZYxMTE2Gz2fDyyy/jhhtuwI8//og33nijyWNmzZqFkSNHolevXigpKcG3336LPn36AAC6d+8ORVHwxRdf4Prrr4fFYmm2SEN8fDymTJmCu+66Cy+99BIuuOACZGdno7CwEBMmTGi1HqcaNWoULrnkErz55ptN6vHUU09h9OjRiIuLa/U9io+PxzfffIMhQ4YgICDA6dHFtuCIkgcZjUb07t3b+w+Q5BbMX27Mn9gG5Ka7/CdMALp1AzIzGzpDLbHbG+7v1g0YP77dXnr58uW49tprW+yM3HLLLdizZw92796NW2+9FU899RTmzJmDCy+8ENnZ2fjzn//c5PErVqxASUkJUlJSHCvORURENHveZ555Bh988AH69++P1atX491333UsDGEymfDSSy/hzTffRExMDMaMGdNq2a+++mqEh4fj4MGDmDRpkmO/oii48cYb8cUXX2DTpk246KKLcOmll2LJkiXo3r37Gd+P/fv3Izo6usmttd8ZMGAAlixZgkWLFuH888/Hu+++i4ULFzZ5jN1ux3333Yc+ffrguuuuw3nnnedYLTA2NhbPPPMMHn30UURGRuL+++9v8XVef/11jBs3DjNnzkTv3r0xY8aMVleya82iRYuaLec+YsSIs75HL7zwAjZt2oS4uDikpKS49JquUsSp6yz6oPLycoSGhqKsrAwdOnTQtCxCCNTV1cHf318/F5yjdsP85cb8iW1AblrlX1tbi8zMTCQkJLh+jsfu3Q0Xkz16tGEJ8E6d/jh36eTJhpGkbt2ARYuAgQPdUwEPUBQFn376qeOaSO4ghIAQAoqi8PPvAWdq9670DTii5EGqqiI7O1s/c5SpXTF/uTF/YhuQmy7zHzgQWLYMmD4dCAoCcnKAw4cb/g0Kati/bJmuO0me5MwCEuRdeI6SBxmNRvTq1UvrYpBGmL/cmD+xDchNt/nHxwOPPALcfXfDdZKqqxtWt0tObrdzkmSgKIomq7bRuWFHyYOEEKitrYXZbOawq4SYv9yYP7ENyE33+YeGApdfrnUp3MITZ6Fw6p0+ceqdB6mqitzcXH0Nu1O7Yf5yY/7ENiA35k+ceqc/HFHyIKPR6PQFw8j3MH+5MX9iG5Cb1vn7+NpdXo9T7zyrvdq7piNK8+bNcwxBNt6ioqIc9wshMG/ePMTExMBisSA1NRX79+/XsMTnRgiByspKHqwkxfzlxvyJbUBuWuXv5+cHAKiurvbo61JTQgjY7XZ+/j2ksb03tv+20nxEqV+/fvj6668dP596fYHFixdjyZIlWLVqFXr16oX58+dj2LBhOHjwIEJCQrQo7jkRQqCwsBDx8fGcnyoh5i835k9sA3LTKn+j0YiwsDAUFhYCAAIDA9n+NCCEQH19Pfz8/Pj+u5EQAtXV1SgsLERYWNg5X7dM846SyWRqMorUSAiBpUuX4oknnsDNN98MAFi9ejUiIyPx3nvv4Z577vF0Uc+ZwWBocrVokgvzlxvzJ7YBuWmZf+P/sxo7S0S+LiwsrMX+has07yilp6cjJiYGAQEBuOSSS7BgwQL06NEDmZmZyM/Px/Dhwx2PDQgIwNChQ7Ft27ZWO0pWqxVWq9Xxc3l5OQA4Tp5s/NdgMJxxu3EqYGvbdrsdBoOh1e3G1zp9u7y8HEFBQTCZTBBCQFVVGI1Gx2ooBoOh1W1ny+7JOjWWvbVt1umPbbvdjoqKCoSGhjqG3vVeJ1/MyV11UhQFlZWVCAwMdHzDpfc6+WJO7qyTqqqoqqpCUFAQFEXxiTr5Yk7uqhMAxwUuG/8meLJO0dHR6Ny5M+rr65mTBnWy2+2ora1FYGCg4/X1XidvzclgMDhG7lqqkyvTHzU9R+mSSy7B22+/jY0bN2LZsmXIz8/H4MGDcfLkSeTn5wMAIiMjm/xOZGSk476WLFy4EKGhoY5bXFwcgD++RSkqKkJRUREAoKCgAMXFxQCAvLw8lJaWAgByc3NRVlYGAMjJyUFFRQUAICsrC1VVVQCAzMxM1NbWAgAyMjIcK5mkp6fDZrNBVVWkp6dDVVXYbDakp6dDCIETJ04gIyMDwB9XDQaAqqoqZGVlAQAqKiqQk5MDoOGgmpubCwAoLS1FXl4eAKC4uBgFBQWa1wloWMWFdTp7nU6cOIGcnBwIIXymTr6Yk7vqZLVaUVJS4lN18sWc3FmnY8eOoaSkBCUlJT5TJ1/MyV11EkLg8OHDmtapoqICJSUlMJvNqK6uRllZGcxmMyorK1FZWQmz2YyysjJUV1fDbDajpKTEsaT5yZMnYbVaYTabceLECdTX18NsNiM/Px92ux1ms9nxfpnNZuTm5sJgMMBsNiMnJwcmkwn+/v7IycmBv78/TCYTcnJyYDabYTAYkJub61jsIC8vD2azGXa7Hfn5+TCbzaivr8eJEydgNpthtVpx8uRJmM1m1NbW6qZOVVVVPlcnb8wpKyvrrJ8nZynCi84qq6qqQmJiImbPno1LL70UQ4YMwfHjxxEdHe14zIwZM5CTk4MNGza0+BwtjSjFxcWhpKQEYWFh7I2zTqwT68Q6sU6sE+vEOrFOrJOkdSovL0dYWJhjhPdMvKqjBADDhg1Dz5498cgjjyAxMRG7d+9GSkqK4/4xY8YgLCwMq1evdur5ysvLERoa6tSb4W5CCJSVlSE0NBSKwhP5ZMP85cb8iW1Absxfbszfe7jSN/CqC85arVb8/vvviI6ORkJCAqKiorBp0ybH/XV1ddiyZQsGDx6sYSnbTgjhGH4n+TB/uTF/YhuQG/OXG/PXJ00Xc/jb3/6GG264Ad26dUNhYSHmz5+P8vJyTJkyBYqiYNasWViwYAGSkpKQlJSEBQsWIDAwEJMmTdKy2G1mMBgc50yRfJi/3Jg/sQ3IjfnLjfnrk6YdpWPHjmHixIkoKipCly5dcOmll2L79u3o3r07AGD27NmoqanBzJkzUVJSgksuuQRfffWVLq+hBDTMjSwtLUVYWJhjviTJg/nLjfkT24DcmL/cmL8+ed05Su3Nm85RUlUVeXl5iI6O5odEQsxfbsyf2Abkxvzlxvy9hyt9A3aUiIiIiIhICrpdzMHXqaqKoqIix5KHJBfmLzfmT2wDcmP+cmP++sSOkofV19drXQTSEPOXG/MntgG5MX+5MX/94dQ7IiIiIiKSAqfeeSlVVVFYWMhhV0kxf7kxf2IbkBvzlxvz1yd2lIiIiIiIiE7DqXdERERERCQFTr3zUo1r6HPYVU7MX27Mn9gG5Mb85cb89YkdJQ/z8/PTugikIeYvN+ZPbANyY/5yY/76w6l3REREREQkBU6981KqqiI3N5fDrpJi/nJj/sQ2IDfmLzfmr0/sKHmYxWLRugikIeYvN+ZPbANyY/5yY/76w6l3REREREQkBU6981KqqiInJ4fDrpJi/nJj/sQ2IDfmLzfmr0/sKHmQoigICQmBoihaF4U0wPzlxvyJbUBuzF9uzF+fTFoXQCaKoiAsLEzrYpBGmL/cmD+xDciN+cuN+esTR5Q8SFVVZGdnc9hVUsxfbsyf2Abkxvzlxvz1iR0lD1IUBR07duSwq6SYv9yYP7ENyI35y4356xOn3nmQoihceU9izF9uzJ/YBuTG/OXG/PWJI0oepKoqjhw5wmFXSTF/uTF/YhuQG/OXG/PXJ3aUPEhRFERERHDYVVLMX27Mn9gG5Mb85cb89YlT7zxIURQEBwdrXQzSCPOXG/MntgG5MX+5MX994oiSB9ntdhw+fBh2u13ropAGmL/cmD+xDciN+cuN+esTO0oeZDAYEBsbC4OBb7uMmL/cmD+xDciN+cuN+esTp955kKIosFgsWheDNML85cb8iW1Absxfbsxfn9it9SC73Y5Dhw5x2FVSzF9uzJ/YBuTG/OXG/PVJEUIIrQvhTuXl5QgNDUVZWZnm69cLIVBXVwd/f3+ueiIh5i835k9sA3Jj/nJj/t7Dlb4Bp955kKIoCAgI0LoYpBHmLzfmT2wDcmP+cmP++sSpdx5kt9tx4MABDrtKSur8S0uBrVuBr75q+Le0VOsSeZzU+RMAtgHZMX+5MX994tQ7DxJCwGazwWQycdhVQlLmn5UFfPghsH49UFQE2GyAyQR07gyMGgVMmADEx2tdSo+QMn9qgm1Absxfbszfe3DqnRfjspBykyr/3buBOXOAo0dhC+qA8uA42GGCETZ0qDgJ0/LlwKZNwKJFwMCBWpfWI6TKn5ooKCjA5s2bUVFRgZCQEKSmpiIyMlLrYpGH8RggN+avP0zMg1RVRXp6OlRV1boopAGp8s/KAubMQV3mMRxWErEjOwr7fvND2n4F+37zw47sKBxWElGXeayhM5WVpXWJ3U6q/MkhLS0NEydORNeuXXH77bfjm2++we23346uXbti4sSJSEtL07qI5CE8BsiN+esTp955kBACqqrCYDBw2FVCUuW/eDFqX12OPVWJqK41ws8E+PkDBgVQBVBfB9TbgECzHQOCM2C+bzrwyCNal9qtpMqfAAAbN27E2LFjYbPZYLPZAAAmk6nJtslkwrp16zBixAgti0oewGOA3Ji/93Clb8ARJQ/jNwlykyL/0lLUfrweWSUdUFNrRFAgEBDQ0EkCGv4NCACCAoGaWiOyijugZu16oKxM23J7gBT5E4CGkaSxY8fCarU6OkYA4Ofn59i22WywWq0YO3YsR5YkwWOA3Ji//rCj5EGqqiIjI4MfFElJk39aGiqzilBg74TAQKC1L84UBQgMBArsnVCVdQLw8f8oSpM/AQAWLFgAm82GUydtmEwmjBkzBibTH6cHN57gvXDhQi2KSR7EY4DcmL8+ceodEbWryk++QvHtDyDbLwkB5rNPL7DWCnSvP4xO776EoJuGe6CERO5VUFCArl27NhlJOhuTyYTc3FxERES4sWRERMSpd15KCAGr1Qof75tSK2TJP+O4BXV2E8x+zv0n0exnQ53diIy8QDeXTFuy5E/A5s2bW+wkKYqC0NDQFs9PsNls2Lx5swdKR1rhMUBuzF+f2FHyIFVVkZ2dzWFXScmS/8mYZBQbOyPUdtKpx4faTqLY2AVF0cluLpm2ZMmfgIqKihb3G41GXHvttTAajS3eX15e7s5ikcZ4DJAb89cndpQ8yGg0olevXq3+kSTfJkv+fl3CsDVkFALry6GIM1+BXBF2WOrL8X3IKPh3CfVQCbUhS/4EhISEtLjfZrPh448/bnVKHqeH+zYeA+TG/PWJHSUPEkKgpqaGw66SkiX/5GRgR/wE5Bq7Iao6s9XOkiLsiKrOxHFjN+xKGI9k3x5QkiZ/AlJTU5ss2NBIURR07ty5xal3JpMJqampHigdaYXHALkxf31iR8mDVFVFbm4uh10lJUv+YWHAhbfE4/mOi3DC3BXR1RkIs+bDqNYDQsCo1iPMmo/o6gycMHfF8x0X4cJb4hHq2wNK0uRPQGRkJMaNG9ess2Q0GjF48OBm3yibTCaMHz+eCzn4OB4D5Mb89Ymr3hFRu8vKAmbMANQjWZjk9xEuPrEeoXUnYBB2qIoRZf5dsKPLKLxXPx6GHvFYtgyIj9e61ETtJy0tDRdffPFZT95WFAUBAQHYsWMHkn19WJWIyAu40jdgR8mDhBCoqqpCUFAQr8osIdny370bmDMHOHoUiA4sw4X+aQhENaoRiF/qkpFXHYpu3YBFi4CBA7UurfvJlj8BGzduxNixY2Gz2WCz2aAoCiIjI1FQUAAhBEwmE0wmE9atW4cRI0ZoXVxyMx4D5Mb8vQeXB/dSQggUFhZyfqqkZMt/4EBg2TJg+nQAoaHYUHk5Piobjg2VlwOhoZg+veF+GTpJgHz5EzBixAjs2LED48ePh8lkgtFoREpKCoxGo2O63Y4dO9hJkgSPAXJj/vrEESUicruyMiAtDaiuBgIDGxZ88PVzkohOVVhYiM2bN6O8vBwdOnRAamoqz0kiItIAp96dwps6SkIIVFRUICQkhMOuEmL+cmP+xDYgN+YvN+bvPTj1zksJIVBSUsJhV0kxf7kxf2IbkBvzlxvz1yeOKBERERERkRQ4ouSlhBAoLS3ltwmSYv5yY/7ENiA35i835q9P7Ch5UOP8VH5I5MT85cb8iW1Absxfbsxfnzj1joiIiIiIpMCpd15KVVUUFxdDVVWti0IaYP5yY/7ENiA35i835q9P7Ch5WE1NjdZFIA0xf7kxf2IbkBvzlxvz1x9OvSMiIiIiIilw6p2XUlUVRUVFHHaVFPOXG/MntgG5MX+5MX99YkfJw+rr67UuAmmI+cuN+RPbgNyYv9yYv/5w6h0REREREUmBU++8lKqqKCws5LCrpJi/3Jg/sQ3IjfnLjfnrEztKREREREREp+HUOyIiIiIikgKn3nkpVVWRl5fHYVdJMX+5MX9iG5Ab85cb89cnr+koLVy4EIqiYNasWY59QgjMmzcPMTExsFgsSE1Nxf79+7UrZDvw8/PTugikIeYvN+ZPbANyY/5yY/764xUdpZ07d+Kf//wn+vfv32T/4sWLsWTJErzyyivYuXMnoqKiMGzYMFRUVGhU0nNjMBjQuXNnGAxe8baThzF/uTF/YhuQG/OXG/PXJ83TqqysxOTJk7Fs2TJ07NjRsV8IgaVLl+KJJ57AzTffjPPPPx+rV69GdXU13nvvPQ1L3HaqqiI3N5fDrpJi/nJj/sQ2IDfmLzfmr0+ad5Tuu+8+jBo1Ctdee22T/ZmZmcjPz8fw4cMd+wICAjB06FBs27at1eezWq0oLy9vcgPgaJiqqjq13bjGRWvbdrv9jNtCiGbbAGA2mx3bp+4XQjhev7VtZ8vu6TqdaZt1arptNpt9rk6+mJO76mSxWHyuTr6YkzvrZLFYfK5OvpiTu+oUEBDgc3XyxZzcVSeLxeJzddJrTs7StKP0wQcfYPfu3Vi4cGGz+/Lz8wEAkZGRTfZHRkY67mvJwoULERoa6rjFxcUBAAoLCwEARUVFKCoqAgAUFBSguLgYAJCXl4fS0lIAQG5uLsrKygAAOTk5jql+WVlZqKqqAtDQkautrQUAZGRkoK6uDgCQnp4Om80GVVWRnp4OVVVhs9mQnp4Og8GA4OBgZGZmAgBqa2sd21VVVcjKygIAVFRUICcnBwBQVlaG3NxcAEBpaSny8vIAAMXFxSgoKNC8TgBQV1eHjIwM1uksdSouLobdbofBYPCZOvliTu6qk81mQ3h4uGPbF+rkizm5s055eXkIDw9HeXm5z9TJF3NyV50MBgPKy8tRU1PjM3XyxZzcVafs7GyEh4ejpqbGZ+qk55ycpdny4Dk5ORg0aBC++uorXHDBBQCA1NRUDBgwAEuXLsW2bdswZMgQHD9+HNHR0Y7fmzFjBnJycrBhw4YWn9dqtcJqtTp+Li8vR1xcHEpKShAWFubonRoMhjNuK4oCRVFa3W78D29r20BDT/fUbUVRcOzYMURHR8PPz8/R0zYajY4ersFgaHXb2bJ7sk6NZW9tm3X6Y9tms+H48ePo2rWro33qvU6+mJO76gQAx48fR1RUFEwmk0/UyRdzcmed7HY78vLyEB0dDYPB4BN18sWc3FUnIQRycnLQtWtXGI1Gn6iTL+bkrjrZbDbk5+cjJiYGiqL4RJ30mlN5eTnCwsKcWh5cs47SunXrcNNNN8FoNDr22e12R+M5ePAgevbsid27dyMlJcXxmDFjxiAsLAyrV6926nW86TpKQgiUlZUhNDQUiqJoWhbyPOYvN+ZPbANyY/5yY/7eQxfXUbrmmmuQlpaGPXv2OG6DBg3C5MmTsWfPHvTo0QNRUVHYtGmT43fq6uqwZcsWDB48WKtinxNFURAWFsYPiKSYv9yYP7ENyI35y43565NmHaWQkBCcf/75TW5BQUHo1KkTzj//fChKwzWVFixYgE8//RT//e9/MXXqVAQGBmLSpElaFfucqKqK7Oxsx3AiyYX5y435E9uA3Ji/3Ji/Ppm0LsCZzJ49GzU1NZg5cyZKSkpwySWX4KuvvkJISIjWRWsTRVHQsWNHfpsgKeYvN+ZPbANyY/5yY/76pNk5Sp7iTecoERERERGRdnRxjpKMVFXFkSNHOOwqKeYvN+ZPbANyY/5yY/76xI6SBymKgoiICA67Sor5y435E9uA3Ji/3Ji/Pnn1OUq+RlEUBAcHa10M0gjzlxvzJ7YBuTF/uTF/feKIkgfZ7XYcPnzYcfFJkgvzl1dBQQHWrFmDd999F2vWrHFckZzkwmOA3Ji/3Ji/PnFEyYMMBgNiY2MdVwgmuTB/+aSlpWHBggVYu3Yt7HY7OnXqhJMnT8JoNGLcuHF4/PHHkZycrHUxyUN4DJAb85cb89cnrnpHROQGGzduxNixY2Gz2WCz2ZrdbzKZYDKZsG7dOowYMUKDEhIREcmHq955KbvdjkOHDnHYVVLMXx5paWkYO3YsrFaro5NkMplwyy23wGRqGMi32WywWq0YO3Ys0tLStCwueQiPAXJj/nJj/vrEjpIHGQwGdO/encOukmL+8liwYAFsNhtOHbC32+34+uuvm/yRFELAZrNh4cKFWhSTPIzHALkxf7kxf33i1DsionZUUFCArl27tjjdrjUmkwm5ubmIiIhwY8mIiIiIU++8lN1ux4EDBzjsKinmL4fNmze3ek7Sbbfd5ph6dyqbzYbNmzd7oHSkJR4D5Mb85cb89YkdJQ8yGAxITEzksKukmL8cKioqWtxvs9nw2WeftTrSVF5e7s5ikRfgMUBuzF9uzF+fmJaH8QMiN+bv+0JCQlq9r76+vtX7ODVYDjwGyI35y4356w8T8yBVVZGeng5VVbUuCmmA+cshNTW1xel1JpMJ48aNa/W+1NRUD5SOtMRjgNyYv9yYvz5xMQcPEkJAVVUYDAYoiqJpWcjzmL88Jk6ciLVr1zabZmcymVrcN378eLz33nueLCJpgMcAuTF/uTF/78HFHLwYv0mQG/OXw+OPPw6TydTsj6Gfn1+TnxVFgclkwmOPPebJ4pGGeAyQG/OXG/PXH3aUPEhVVWRkZPCDIinmL4/k5GSsW7cOAQEBjql2JpMJY8aMafJzQEAA1q1bh+TkZC2LSx7CY4DcmL/cmL8+ceodEZGbpKWlYeHChfjoo4+aTLlrnG732GOPsZNERETkQa70DdhR8iAhBOrq6uDv78/5qRJi/vIqLCzE5s2bUVFRgZCQEKSmpvLishLiMUBuzF9uzN978BwlL6WqKrKzsznsKinmL6+IiAjccsstuOKKK3DLLbewkyQpHgPkxvzlxvz1qU0jShkZGVi5ciUyMjLw4osvIiIiAhs2bEBcXBz69evnjnK2mTeNKBERERERkXbcOqK0ZcsWJCcn4+eff8Ynn3yCyspKAMC+ffvw9NNPt63EkhBCoKamBj4+25FawfzlxvyJbUBuzF9uzF+fXO4oPfroo5g/fz42bdoEf39/x/6rrroKP/30U7sWzteoqorc3FwOu0qK+cuN+RPbgNyYv9yYvz65PPUuODgYaWlpSEhIQEhICPbu3YsePXogKysLvXv3Rm1trbvK2iacekdERERERICbp96FhYUhLy+v2f5ff/0VsbGxrj6dVIQQqKys5LCrpJi/3Jg/sQ3IjfnLjfnrk8sdpUmTJmHOnDnIz8+HoihQVRU//vgj/va3v+GOO+5wRxl9hhAChYWF/JBIivnLjfkT24DcmL/cmL8+uTz1rr6+HlOnTsUHH3wAIQRMJhPsdjsmTZqEVatWwWg0uqusbcKpd0REREREBHjogrMZGRn49ddfoaoqUlJSkJSU1KbCups3dZSEEI4LTvJiY/Jh/nJj/sQ2IDfmLzfm7z1c6RuY2voiiYmJSExMbOuvS0kIgZKSEgQHB/NDIiHmLzfmT2wDcmP+cmP++uTyiNJdd911xvtXrFhxTgVqb940okRERERERNpx64hSSUlJk5/r6+vx3//+F6Wlpbj66qtdfTqpCCFQVlaG0NBQfpsgIeYvN+ZPbANyY/5yY/765HJH6dNPP222T1VVzJw5Ez169GiXQvmqxvmpHTp04IdEQsxfbsyf2Abkxvzlxvz1qc2LOZzu4MGDSE1NbfEaS1ri1DsiIiIiIgLcfMHZ1mRkZMBms7XX0/kkVVVRXFwMVVW1LgppgPnLjfkT24DcmL/cmL8+uTz17uGHH27ysxACeXl5WL9+PaZMmdJuBfNVNTU1CAsL07oYpBHmLzfmT2wDcmP+cmP++uPy1Lurrrqqyc8GgwFdunTB1VdfjbvuugsmU5tXHHcLTr0jIiLNlZYCaWlATQ1gsQDJyQD/w0RE5HFuXfXuu+++a3PBZNc47BoeHg6Dod1mPZJOMH+5MX9JZWUBH34IrF8PtbgYxV27IvzYMRjCw4FRo4AJE4D4eK1LSR7AY4DcmL8+MSkPq6+v17oIpCHmLzfmL5ndu4EZM4Dly4HqaiAuDvUJCUBcXMPPy5c33L97t9YlJQ/hMUBuzF9/nJp6l5KS4vRShru97IDPqXdERORxWVkNnaBjx4CEBNSrRpRXAKodMBiBDiGAn8EOZGYCXbsCy5ZxZImIyAPaferd2LFj26Nc0lNVFUVFRejcuTOHXSXE/OXG/CXz4YfA0aOojknE8SwjCgsAq01Bfd8Y+P33OAJMAhGRRsREJSDwaAbw0UfAI49oXWpyIx4D5Mb89cmpjtLTTz/t7nIQERH5htJSYP16VJk6YN9/jaipAfxMgNkCmMyAyQLYaoGjR4ETJ4zo36UDgtavB+6+GwgN1br0RET0P+zSepDBYEBERAS/SZAU85cb85dIWhrq8oqwP78TamuAoEAgIAAwCoGA33Ib/g1o2F9bA+wv6IS64ycaVsUjn8VjgNyYvz65nJbdbsf//d//4eKLL0ZUVBTCw8Ob3Kh1qqoiLy+PFxuTFPOXG/OXSE0NKkpsqKw1ITAQaDzFVxgV1F7QDcLYsENRgMBAoLLWhIpSe8MCD+SzeAyQG/PXJ5c7Ss888wyWLFmCCRMmoKysDA8//DBuvvlmGAwGzJs3zw1F9C1+fn5aF4E0xPzlxvzlUGGzoLzaBLPJhibrIAlAqakHTllCSVEAs8mG8iojKtVAj5eVPIvHALkxf/1xuaP07rvvYtmyZfjb3/4Gk8mEiRMn4q233sJTTz2F7du3u6OMPsNgMPAkPokxf7kxf3mkIRlF6IxOONlkv6IKBBzKg6I2XWy2E07iBLogDcmeLCZ5GI8BcmP++uRyWvn5+UhObjiYBwcHo6ysDAAwevRorF+/vn1L52NUVUVubi6HXSXF/OXG/OVRaQrDd4GjEGQvhyLsjv3CqKDmwgTH1DsAUIQdQfZybA4chQoDF3LwZTwGyI3565PLHaWuXbsiLy8PANCzZ0989dVXAICdO3ciICCgfUvngywWi9ZFIA0xf7kxfzlYLMCmsAkoDOiGqOrMPzpLAjCWVDmm3inCjqjqTBQGdMOmsPEI5Mw7n8djgNyYv/643FG66aab8M033wAA/vKXv2Du3LlISkrCHXfcgbvuuqvdC+hLDAYDwsPDOewqKeYvN+Yvj+RkoC4mHksiF6HI0hXR1RkIs+bDZKuDf0YBTLY6hFnzEV2dgSJLVyyJXIT62Hgkc+adT+MxQG7MX58UIYQ4+8OApUuX4o477mi2st327duxbds29OzZEzfeeKNbCnkuXLn6rrs1DrvGxsbygyIh5i835i+XxYuB5cuBy6KzcGXBR7iocD1C7MUou7QfQrfvR4UxHDsjRuH7yPH4KS8e06fzerO+jscAuTF/7+FK38DpjlLHjh1RU1ODMWPGYNq0aRg2bBiUJsv5eCdv6igJIVBWVobQ0FBdvHfUvpi/3Ji/XLKygBkzgGPHgIQEIEQtQ3zlPoR1qkLpySBkBfdHhSEUmZlA167AsmVAfLzWpSZ34jFAbszfe7ilo2S1WrF27VqsXLkS3333HWJjY3HnnXdi6tSpSEhIaJeCu4M3dZSIiEgeu3cDc+YAR48CHToAnToBJhNgswEnTwLl5UC3bsCiRcDAgVqXlohIDq70DZwe+wsICMDkyZPx9ddfIyMjA3feeSfefvttJCUl4dprr8X7778Pq9V6zoX3ZaqqIjs7myueSIr5y435y2fgwIaRounTgaAg4PhxFR07ZuP4cRVBQQ37ly1jJ0kWPAbIjfnrk9MjSq35+uuvsXLlSqxbtw5msxknT548+y95kDeNKAkhUFFRgZCQEA67Soj5y435y62sDNi3T6CqqgJBQSHo319BKFcDlwqPAXJj/t7Dlb6B6VxfzGAwQFEUCCHYSz4LRVE076yRdpi/3Ji/3EJDgSuuUACwDciKxwC5MX99atOyG9nZ2XjmmWeQkJCA4cOH4/jx41i2bJnj+krUMlVVceTIEXYoJcX85cb8iW1Absxfbsxfn5weUaqtrcXHH3+MFStWYMuWLYiOjsaUKVNw1113oUePHu4so89QFAUREREccpUU85cb8ye2Abkxf7kxf31yuqMUFRWF2tpajB49Gv/+978xYsQIrgPvIkVREBwcrHUxSCPMX27Mn9gG5Mb85cb89cnpns5TTz2FY8eOYe3atRg5ciQ7SW1gt9tx+PBh2O12rYtCGmD+cmP+xDYgN+YvN+avT+e86p2387ZV72pra2E2mzn0KiHmLzfmT2wDcmP+cmP+3sOjq96R8xRFgcVi0boYpBHmLzfmT2wDcmP+cmP++sT5cx5kt9tx6NAhDrtKivnLjfkT24DcmL/cmL8+adpRev3119G/f3906NABHTp0wGWXXYYvv/zScb8QAvPmzUNMTAwsFgtSU1Oxf/9+DUt8bgwGA7p3787zuyTF/OXG/IltQG7MX27MX59cTuvZZ59FdXV1s/01NTV49tlnXXqurl274vnnn8euXbuwa9cuXH311RgzZoyjM7R48WIsWbIEr7zyCnbu3ImoqCgMGzYMFRUVrhbbKyiKgoCAAM5NlRTzlxvzJ7YBuTF/uTF/fXJ5MQej0Yi8vDxEREQ02X/y5ElERESc85BieHg4/vGPf+Cuu+5CTEwMZs2ahTlz5gAArFYrIiMjsWjRItxzzz1OPZ83LeZgt9uRnp6OpKQkGI1GTctCnsf85cb8iW1Absxfbszfe7jSN3B5REkI0WJveO/evQgPD3f16Rzsdjs++OADVFVV4bLLLkNmZiby8/MxfPhwx2MCAgIwdOhQbNu2rdXnsVqtKC8vb3ID4LgSsqqqTm039h9b27bb7WfcFkI02zYYDOjRo4fjsY37G7cbX7+1bWfL7sk6tVQP1qnl8gJAQkICDAaDz9TJF3NyV50URUFiYqKjDr5QJ1/MyZ11AoDExERHGXyhTr6Yk7vqZDAYkJCQ4Pg/lC/UyRdzcledhBBITEyEoig+Uyc95+QspztKHTt2RHh4OBRFQa9evRAeHu64hYaGYtiwYZgwYYLTL9woLS0NwcHBCAgIwL333otPP/0Uffv2RX5+PgAgMjKyyeMjIyMd97Vk4cKFCA0Nddzi4uIAAIWFhQCAoqIiFBUVAQAKCgpQXFwMAMjLy0NpaSkAIDc3F2VlZQCAnJwcx1S/rKwsVFVVAQAyMzNRW1sLAMjIyEBdXR0AID09HTabDaqqIj09HaqqwmazIT09HQBgs9lw5MgRAEBtbS0yMzMBAFVVVcjKygIAVFRUICcnBwBQVlaG3NxcAEBpaSny8vIAAMXFxSgoKPCKOtXV1SEjI4N1cqJOjfXwpTr5Yk7uqpPBYPC5OvliTu6sk8Fg8Lk6+WJO7qrTsWPHfK5OvpiTO+qUlZUFg8HgU3XSc07Ocnrq3erVqyGEwF133YWlS5ciNDTUcZ+/vz/i4+Nx2WWXOf3Cjerq6nD06FGUlpbi448/xltvvYUtW7agtLQUQ4YMwfHjxxEdHe14/IwZM5CTk4MNGza0+HxWqxVWq9Xxc3l5OeLi4lBSUoKwsDBH77TxW/3WthVFcfT6W9puHCFqbRto6Omeug0Ahw4dQmJiIvz9/R09baPR6OjhGgyGVredLbsn69RY9ta2Wac/tuvr63H48GH06tXL8Y2i3uvkizm5q05CCBw+fBg9evSAn5+fT9TJF3NyZ51sNhsyMjKQmJgIo9HoE3XyxZzcVSdVVXHo0CEkJSXBZDL5RJ18MSd31am+vh5HjhxBz549YTAYfKJOes2pvLwcYWFhTk29c/kcpS1btmDw4MGOP/Tt7dprr0ViYiLmzJmDxMRE7N69GykpKY77x4wZg7CwMKxevdqp5/Omc5QaG1VjeCQX5i835k9sA3Jj/nJj/t7DrRecHTp0qONbkcLCQkePr9GVV17p6lM2IYSA1WpFQkICoqKisGnTJkdHqa6uDlu2bMGiRYvO6TW0dGoPneTD/OXG/IltQG7MX27MX39c7iht374dkyZNQnZ2drOToRqHuJz1+OOPY+TIkYiLi0NFRQU++OADbN68GRs2bICiKJg1axYWLFiApKQkJCUlYcGCBQgMDMSkSZNcLbZXUFUVGRkZXPFEUsxfYqWlUPftQ0ZFBZJCQmDs3x8IC9O6VORhPAbIjfnLjfnrk8sdpXvvvReDBg3C+vXrER0dfU7DhwUFBfjTn/6EvLw8hIaGon///tiwYQOGDRsGAJg9ezZqamowc+ZMlJSU4JJLLsFXX32FkJCQNr+mloxGI3r37q11MUgjzF9CWVnAhx8C69fDWFSE3jYbYDIBnTsDo0YBEyYA8fFal5I8hMcAuTF/uTF/fXL5HKWgoCDs3bsXPXv2dFeZ2pW3naNUV1cHf39/zk+VEPOXzO7dwJw5wNGjsAV1QJlfJ9QGBsNcXYnQ+pMwVZUD3boBixYBAwdqXVryAB4D5Mb85cb8vYdbr6N0ySWX4PDhw20unMxUVUV2dnaz87pIDsxfIllZwJw5qMs8hsNKInZkR2HfoQDsi+6NfYcCsCM7CoeVRNRlHmvoTP1vWVXybTwGyI35y43565PLI0qffvopnnzySTzyyCNITk5utvpd//7927WA58qbRpSISBKLF6P21eXYU5WI6loj/EyAnz9gUABVAPV1QL0NCDTbMSA4A+b7pgOPPKJ1qYmIiHyeK30DlztKLa3WoSgKhBAuL+bgCd7UURJCoLa2FmazmcOuEmL+kigtRe2IMcj6vRp5ahQCAwFFAYQCqGFBMJRWQRGAEEB1NRBtyEf3PkGwfPUZcMr16cj38BggN+YvN+bvPdw69S4zM7PZ7ciRI45/qXWqqiI3N5fDrpJi/pJIS0NlVhEK7J0cnSQAgNGAmgsTAGPDYVdRgMBAoMDeCVVZJ4C0NO3KTB7BY4DcmL/cmL8+ubzqXffu3d1RDikYjUbdLIJB7Y/5y6HyRA2qK2ww+Jlw6peGik1F8Nf/bfJYRQEMJhOqK+yoOlGNIA+XlTyLxwC5MX+5MX99atNVr/71r39hyJAhiImJQXZ2NgBg6dKl+Oyzz9q1cL5GCIHKyspm158iOTB/OWQct6DOboLZz9Zkv1AAW5cOEKfNuDD72VBnNyIjL9CDpSQt8BggN+YvN+avTy53lF5//XU8/PDDuP7661FaWuo4JyksLAxLly5t7/L5FCEECgsL+SGRFPOXw8mYZBQbOyPUdrLpHQYDrP1igdPO8wy1nUSxsQuKopM9WErSAo8BcmP+cmP++uRyR+nll1/GsmXL8MQTTzS5svCgQYOQxjn2Z2QwGNCjR48WF8Qg38f85eDXJQxbQ0YhsL4civhjcRvFriJo8+9Q7H/MT1eEHZb6cnwfMgr+XbiQg6/jMUBuzF9uzF+f2rSYQ0pKSrP9AQEBqKqqapdC+SohBMrLy/ltgqSYvxySk4Ed8ROQa+yGqOpMR2dJKArqo8Mg/nfikiLsiKrOxHFjN+xKGI9kDij5PB4D5Mb85cb89cnljlJCQgL27NnTbP+XX36Jvn37tkeZfJYQAiUlJfyQSIr5yyEsDLjwlng833ERTpi7Iro6A2HWfBhhQ32PCBhhQ5g1H9HVGThh7ornOy7ChbfEc2VwCfAYIDfmLzfmr08ur3r3yCOP4L777kNtbS2EENixYwfef/99LFy4EG+99ZY7yugzDAYDVw2UGPOXx4QJwKZNA/HIkWWY1PkjXHxiPSKqshH11RGoihFl/l2wNWoC3qsfD0OPeIwfr3WJyRN4DJAb85cb89cnly84CwDLli3D/PnzkZOTAwCIjY3FvHnzMG3atHYv4LnytgvOlpWVITQ0lBcbkxDzl8vu3cCcOcDRo0B0YBkuDNiHzhFVKCoMwi/W/sirDkW3bsCiRcDAgVqXljyBxwC5MX+5MX/v4UrfoE0dpUZFRUVQVRURERFtfQq386aOUuPFxmJjY3kyn4SYv3yysoCPPgLWrweKi1X065eL/ftjER5uwKhRwPjxQHy81qUkT+ExQG7MX27M33t4rKOkB97UUSIiOZWVAWlpQHU1EBjYsOADz0kiIiLyPFf6Bk6dozRw4EB888036NixI1JSUs44ZLh7927XSisRVVVRWlqKsLAwfpsgIeYvr9BQYPBg5i87HgPkxvzlxvz1yamO0pgxYxAQEAAAGDt2rDvL4/NqamoQFhamdTFII8xfbsyf2Abkxvzlxvz1h1PviIiIiIhICq70DVwe+9u5cyd+/vnnZvt//vln7Nq1y9Wnk4qqqo4FMEg+zF9uzJ/YBuTG/OXG/PXJ5Y7Sfffd51gW/FS5ubm477772qVQvqy+vl7rIpCGmL/cmD+xDciN+cuN+euPy1PvgoODsW/fPvTo0aPJ/szMTPTv3x8VFRXtWsBzxal3REREREQEuHnqXUBAAAoKCprtz8vLg8nk1NoQ0lJVFYWFhRx2lRTzlxvzJ7YBuTF/uTF/fXK5ozRs2DA89thjKCsrc+wrLS3F448/jmHDhrVr4YiIiIiIiLTg8tS73NxcXHnllTh58iRSUlIAAHv27EFkZCQ2bdqEuLg4txS0rTj1joiIiIiIADdPvYuNjcW+ffuwePFi9O3bFxdeeCFefPFFpKWleV0nyduoqoq8vDwOu0qK+cuN+RPbgNyYv9yYvz616aSioKAg3H333e1dFin4+flpXQTSEPOXG/MntgG5MX+5MX/9cWrq3eeff46RI0fCz88Pn3/++Rkfe+ONN7Zb4doDp94RERERERHgWt/AqY6SwWBAfn4+IiIiYDC0PltPURTY7XbXS+xG3tRRahx2jY6OPuP7SL6J+cuN+RPbgNyYv9yYv/dwpW/g1NS7U+dTcm7lubFYLFoXgTTE/OXG/IltQG7MX27MX3+c6tKGh4ejqKgIAHDXXXd53UVl9cJgMCA8PJzfJEiK+cuN+RPbgNyYv9yYvz45lVZdXR3Ky8sBAKtXr0Ztba1bC+WrVFVFTk4OR+UkxfzlxvyJbUBuzF9uzF+fnJp6d9lll2Hs2LG48MILIYTAgw8+2Orw4YoVK9q1gL5EURSEhIRAURSti0IaYP5yY/7ENiA35i835q9PTnWU3nnnHfy///f/kJGRAQAoKyvjqFIbKIqCsLAwrYtBGmH+cmP+xDYgN+YvN+avT06teneqhIQE7Nq1C506dXJXmdqVt616l5OTg7i4OM5RlRDzlxvzJ7YBuTF/uTF/7+FK38DlxRyuuuoq+Pv7n3spJaQoCjp27MhhV0kxf7kxf2IbkBvzlxvz1ycu5uBBiqKgQ4cO/JBIivnLjfkT24DcmL/cmL8+cTEHD1JVFVlZWYiPj+ewq4SYv9yYP7ENyI35y43565PLizkoisLFHNpIURRERETw2wRJMX+5MX9iG5Ab85cb89cnLuZARERERERSaPfFHE6VmZnp6CRxVMk1drsdhw8fht1u17oopAHmLzfmT2wDcmP+cmP++uRyR0lVVTz33HOIjY1FcHAwjhw5AgCYO3culi9f3u4F9CUGgwGxsbGcmyop5i835k9sA3Jj/nJj/vrkclrz58/HqlWrsHjx4ibLhCcnJ+Ott95q18L5GkVRYLFYOD9VUsxfbsyf2Abkxvzlxvz1yeWO0ttvv41//vOfmDx5MoxGo2N///79ceDAgXYtnK+x2+04dOgQh10lxfzlxvyJbUBuzF9uzF+fXO4o5ebmomfPns32q6qK+vr6dimUrzIYDOjevTuHXSXF/OXG/IltQG7MX27MX59cTqtfv37YunVrs/0fffQRUlJS2qVQvkpRFAQEBHDYVVLMX27Mn9gG5Mb85cb89cmp6yid6umnn8af/vQn5ObmQlVVfPLJJzh48CDefvttfPHFF+4oo8+w2+1IT09HUlJSk2mLJAfmLzfmT2wDcmP+cmP++uTydZQAYOPGjViwYAF++eUXqKqKgQMH4qmnnsLw4cPdUcZz4k3XURJCwGazwWQy8RsFCTF/uTF/YhuQG/OXG/P3Hq70DdrUUdITb+soqaoKg8HAD4mEmL/cmD+xDciN+cuN+XsPt15wttEvv/yCd955B++++y5+/fXXtj6NVFRVRXp6OlRV1boopAHmLzfmT2wDcmP+cmP++uTyiFJhYSFuu+02bN68GWFhYRBCoKysDFdddRU++OADdOnSxV1lbROOKJG3YP5yY/7ENiA35i835u893Dqi9MADD6C8vBz79+9HcXExSkpK8N///hfl5eV48MEH21xoWfCbBLkxf7kxf2IbkBvzlxvz1x+XO0obNmzA66+/jj59+jj29e3bF6+++iq+/PLLdi2cr1FVFRkZGfygSIr5y435E9uA3Ji/3Ji/Prm8PLiqqvDz82u238/Pj+GfhdFoRO/evbUuBmmE+cuN+RPbgNyYv9yYvz65PKJ09dVX4y9/+QuOHz/u2Jebm4uHHnoI11xzTbsWztcIIWC1WuHjCw1SK5i/3Jg/sQ3IjfnLjfnrk8sdpVdeeQUVFRWIj49HYmIievbsiYSEBFRUVODll192Rxl9hqqqyM7O5sibpJi/3Jg/sQ3IjfnLjfnrU5uvo7Rp0yYcOHAAQgj07dsX1157bXuXrV1406p3RERERESkHV5w9hTe1FESQqC2thZms5lLQ0qI+cuN+RPbgNyYv9yYv/dwy/Lg3377Lfr27Yvy8vJm95WVlaFfv37YunWr66WViKqqyM3N5bCrpJi/3Jg/sQ3IjfnLjfnrk9MjSjfeeCOuuuoqPPTQQy3e/9JLL+G7777Dp59+2q4FPFfeNKJERERERETaccuI0t69e3Hddde1ev/w4cPxyy+/OF9KCQkhUFlZyRVPJMX85cb8iW1Absxfbsxfn5zuKBUUFLR4/aRGJpMJJ06caJdC+SohBAoLC/khkRTzlxvzJ7YBuTF/uTF/fXK6oxQbG4u0tLRW79+3bx+io6PbpVC+ymAwoEePHjAYXF6VnXwA85cb8ye2Abkxf7kxf31yOq3rr78eTz31FGpra5vdV1NTg6effhqjR4926cUXLlyIiy66CCEhIYiIiMDYsWNx8ODBJo8RQmDevHmIiYmBxWJBamoq9u/f79LreAshBMrLy/ltgqSYv9yYP7ENyI35y43565PTHaUnn3wSxcXF6NWrFxYvXozPPvsMn3/+ORYtWoTzzjsPxcXFeOKJJ1x68S1btuC+++7D9u3bsWnTJthsNgwfPhxVVVWOxyxevBhLlizBK6+8gp07dyIqKgrDhg1DRUWFS6/lDYQQKCkp4YdEUsxfbsyf2Abkxvzlxvz1yaXrKGVnZ+PPf/4zNm7c6AhaURSMGDECr732GuLj48+pMCdOnEBERAS2bNmCK6+8EkIIxMTEYNasWZgzZw4AwGq1IjIyEosWLcI999xz1ufkqndERERERAS4adU7AOjevTv+85//oKioCD///DO2b9+OoqIi/Oc//znnThLQcD0mAAgPDwcAZGZmIj8/H8OHD3c8JiAgAEOHDsW2bdtafA6r1Yry8vImNwCOdetVVXVqu7Ej2Nq23W4/47YQosXtkpIS2Gw2AHDsb9xufP3Wtp0tuyfr1FI9WKeWy2u32x3fJvlKnXwxJ3fVSVVVlJaWwmaz+UydfDEnd9bJbrejtLTU0R58oU6+mJO76iSEQHFxcZO/CXqvky/m5K462Ww2lJaWNnl9vddJzzk5q01nlHXs2BEXXXQRLr74YnTs2LEtT9GMEAIPP/wwLr/8cpx//vkAgPz8fABAZGRkk8dGRkY67jvdwoULERoa6rjFxcUBAAoLCwEARUVFKCoqAtCwkl9xcTEAIC8vD6WlpQCA3NxcR6ctJyfHMc0vKyvLMS0wMzPTcb5WRkYG6urqAADp6emw2WxQVRXp6elQVRU2mw3p6ekQoqGjlJGRAQCora1FZmYmAKCqqgpZWVkAgIqKCuTk5ABo6Dzm5uYCAEpLS5GXlwcAKC4uRkFBgeZ1AoC6ujrWyYk6nThxAsePH4cQwmfq5Is5uatOVqsVFRUVPlUnX8zJnXU6duwYKioqUFJS4jN18sWc3FUnIYTP1ckXc3JnnSoqKnyuTnrNyVkuTb1zp/vuuw/r16/HDz/8gK5duwIAtm3bhiFDhuD48eNNVtSbMWMGcnJysGHDhmbPY7VaYbVaHT+Xl5cjLi4OJSUlCAsLc/RODQbDGbcVRYGiKK1u2+12GAyGVreBhp7uqdtGo9HRu25pWwgBg8HQ6razZWedWCfWiXVinVgn1ol1Yp1YJ9apeZ3Ky8sRFhbm1NQ7r+goPfDAA1i3bh2+//57JCQkOPYfOXIEiYmJ2L17N1JSUhz7x4wZg7CwMKxevfqsz+1N5yipasPUm7CwMEdoJA/mLzfmT2wDcmP+cmP+3sNt5yi1NyEE7r//fnzyySf49ttvm3SSACAhIQFRUVHYtGmTY19dXR22bNmCwYMHe7q47aKmpkbrIpCGmL/cmD+xDciN+cuN+euPScsXv++++/Dee+/hs88+Q0hIiOO8o9DQUFgsFiiKglmzZmHBggVISkpCUlISFixYgMDAQEyaNEnLoreJwWBAbGys1sUgjTB/uTF/YhuQG/OXG/PXJ01HlF5//XWUlZUhNTUV0dHRjtuaNWscj5k9ezZmzZqFmTNnYtCgQcjNzcVXX32FkJAQDUveNqqqoqioyDHvkuTC/OXG/IltQG7MX27MX580HVFy5vQoRVEwb948zJs3z/0F8oD6+nqti0AaYv5yY/7ENiA35i835q8/XrGYgzt502IORERERESkHd0s5iAbVVVRWFjIYVdJMX+5MX9iG5Ab85cb89cndpSIiIiIiIhOw6l3REREREQkBU6981KqqiIvL4/DrpJi/nJj/sQ2IDfmLzfmr0/sKHmYn5+f1kUgDTF/uTF/YhuQG/OXG/PXH069IyIiIiIiKXDqnZdSVRW5ubkcdpUU85cb8ye2Abkxf7kxf31iR8nDLBaL1kUgDTF/uTF/YhuQG/OXG/PXH069IyIiIiIiKXDqnZdSVRU5OTkcdpUU85cb8ye2Abkxf7kxf31iR8mDFEVBSEgIFEXRuiikAeYvN+ZPbANyY/5yY/76ZNK6ADJRFAVhYWFaF4M0wvzlxvyJbUBuzF9uzF+fOKLkQaqqIjs7m8OukmL+cmP+xDYgN+YvN+avT+woeZCiKOjYsSOHXSXF/OXG/IltQG7MX27MX5849c6DFEXhynsSY/5yY/7ENiA35i835q9PHFHyIFVVceTIEQ67Sor5y435E9uA3Ji/3Ji/PrGj5EGKoiAiIoLDrpJi/nJj/sQ2IDfmLzfmr0+ceudBiqIgODhY62KQRpi/3Jg/sQ3IjfnLjfnrE0eUPMhut+Pw4cOw2+1aF4U0wPzlxvyJbUBuzF9uzF+f2FHyIIPBgNjYWBgMfNtlxPzlxvyJbUBuzF9uzF+fOPXOgxRFgcVi0boYpBHmLzfmT2wDcmP+cmP++sRurQfZ7XYcOnSIw66SYv5yY/7ENiA35i835q9PihBCaF0IdyovL0doaCjKyso0X79eCIG6ujr4+/tz1RMJMX+5MX9iG5Ab85cb8/cervQNOPXOgxRFQUBAgNbFII0wf7kxf2IbkBvzlxvz1ydOvfMgu92OAwcOcNhVUsxfbsyf2Abkxvzlxvz1iVPvPEgIAZvNBpPJxGFXCTF/uTF/YhuQG/OXG/P3Hq70DTii5GFcFlJuzF9uzJ/YBuTG/OXG/PWHiXmQqqpIT0+HqqpaF4U0wPzlxvyJbUBuzF9uzF+fOPXOg4QQUFUVBoOBw64SYv5yY/7ENiA35i835u89OPXOi/GbBLkxf7kxf2IbkBvzlxvz1x92lDxIVVVkZGTwgyIp5i835k9sA3Jj/nJj/vrEqXdERERERCQFTr3zUkIIWK1W+HjflFrB/OXG/IltQG7MX27MX5/YUfIgVVWRnZ3NYVdJMX+5MX9iG5Ab85cb89cnTr0jIiIiIiIpcOqdlxJCoKamhsOukmL+cmP+xDYgN+YvN+avT+woeZCqqsjNzeWwq6SYv9yYP7ENyI35y4356xOn3hERERERkRQ49c5LCSFQWVnJYVdJMX+5MX9iG5Ab85cb89cndpQ8SAiBwsJCfkgkxfzlxvyJbUBuzF9uzF+fOPWOiIiIiIikwKl3XkoIgfLycn6bICnmLzfmT2wDcmP+cmP++sSOkgcJIVBSUsIPiaSYv9yYP7ENyI35y4356xOn3hERERERkRQ49c5LCSFQWlrKbxMkxfzlxvyJbUBuzF9uzF+f2FHyICEEKioq+CGRFPOXG/MntgG5MX+5MX994tQ7IiIiIiKSAqfeeSlVVVFcXAxVVbUuCmmA+cuN+RPbgNyYv9yYvz6xo+RhNTU1WheBNMT85cb8iW1Absxfbsxffzj1joiIiIiIpMCpd15KVVUUFRVx2FVSzF9uzJ/YBuTG/OXG/PWJHSUPq6+v17oIpCHmLzfmT2wDcmP+cmP++sOpd0REREREJAVOvfNSqqqisLCQw66SYv5yY/7ENiA35i835q9P7CgRERERERGdhlPviIiIiIhICpx656VUVUVeXh6HXSXF/OXG/IltQG7MX27MX5/YUfIwPz8/rYtAGmL+cmP+xDYgN+YvN+avP5x6R0REREREUtDN1Lvvv/8eN9xwA2JiYqAoCtatW9fkfiEE5s2bh5iYGFgsFqSmpmL//v3aFLYdqKqK3NxcDrtKivnLjfkT24DcmL/cmL8+adpRqqqqwgUXXIBXXnmlxfsXL16MJUuW4JVXXsHOnTsRFRWFYcOGoaKiwsMlbT8Wi0XrIpCGmL/cmD+xDciN+cuN+euP10y9UxQFn376KcaOHQugYTQpJiYGs2bNwpw5cwAAVqsVkZGRWLRoEe655x6nnpdT74iIiIiICNDR1LszyczMRH5+PoYPH+7YFxAQgKFDh2Lbtm2t/p7VakV5eXmTGwDHUKeqqk5tN/YfW9u22+1n3BZCNNtWVRVHjx5FfX09ADj2N243vn5r286W3ZN1aqkerFPL5bXZbDh69Khjny/UyRdzcled7HY7cnJyUF9f7zN18sWc3Fknm82GnJwc2Gw2n6mTL+bkrjqpqors7GxHGXyhTr6Yk7vqVF9fj5ycHMf/B32hTnrOyVle21HKz88HAERGRjbZHxkZ6bivJQsXLkRoaKjjFhcXBwAoLCwEABQVFaGoqAgAUFBQgOLiYgBAXl4eSktLAQC5ubkoKysDAOTk5Dim+mVlZaGqqgpAQ0eutrYWAJCRkYG6ujoAQHp6uuOPYHp6OlRVhc1mQ3p6OhRFgdlsRmZmJgCgtrbWsV1VVYWsrCwAQEVFBXJycgAAZWVlyM3NBQCUlpYiLy8PAFBcXIyCggLN6wQAdXV1yMjIYJ3OUqeTJ0/CZrNBURSfqZMv5uSuOtXX1yMkJASHDx/2mTr5Yk7urNPx48cREhKCsrIyn6mTL+bkrjopioLKykpUV1f7TJ18MSd31Sk7OxshISGorq72mTrpOSdnee3Uu23btmHIkCE4fvw4oqOjHY+bMWMGcnJysGHDhhafx2q1wmq1On4uLy9HXFwcSkpKEBYW5uidGgyGM24rigJFUVrdttvtMBgMrW4DDT3dU7eNRqOjd93SthACBoOh1W1ny846sU6sE+vEOrFOrBPrxDqxTqxT8zqVl5cjLCzMqal3XttROnLkCBITE7F7926kpKQ4HjdmzBiEhYVh9erVTj2vN52jpKoqcnJyEBcX5wiN5MH85cb8iW1Absxfbszfe/jEOUoJCQmIiorCpk2bHPvq6uqwZcsWDB48WMOStZ2iKOjYsSMURdG6KKQB5i835k9sA3Jj/nJj/vpk0vLFKysrcfjwYcfPmZmZ2LNnD8LDw9GtWzfMmjULCxYsQFJSEpKSkrBgwQIEBgZi0qRJGpa67RRF0XxUi7TD/OXG/IltQG7MX27MX580HVHatWsXUlJSHFPrHn74YaSkpOCpp54CAMyePRuzZs3CzJkzMWjQIOTm5uKrr75CSEiIlsVuM1VVceTIEce8S5IL85cb8ye2Abkxf7kxf33ymnOU3MWbzlESQqCqqgpBQUEcepUQ85cb8ye2Abkxf7kxf+/hSt9A06l3slEUBcHBwVoXgzTC/OXG/IltQG7MX27MX5+8djEHX2S323H48GHHBbpILsxfbsyf2Abkxvzlxvz1iR0lDzIYDIiNjeWykJJi/nJj/sQ2IDfmLzfmr0+ceudBiqLAYrFoXQzSCPOXG/MntgG5MX+5MX99YrfWg+x2Ow4dOsRhV0kxf7kxf2IbkBvzlxvz1yeueudBQgjU1dXB39+fK55IiPnLjfkT24DcmL/cmL/34Kp3XkpRFAQEBGhdDNII85cb8ye2Abkxf7kxf33i1DsPstvtOHDgAIddJcX85cb8iW1Absxfbsxfnzj1zoOEELDZbDCZTBx2lRDzlxvzJ7YBuTF/uTF/7+FK34AjSh7GZSHlxvzlxvyJbUBuzF9uzF9/mJgHqaqK9PR0qKqqdVFIA8xfbsyf2Abkxvzlxvz1iVPvPEgIAVVVYTAYOOwqIeYvN+ZPbANyY/5yY/7eg1PvvBi/SZAb85cb8ye2Abkxf7kxf/1hR8mDVFVFRkYGPyiSYv5yY/7ENiA35i835q9PnHpHRERERERS4NQ7LyWEgNVqhY/3TakVzF9uzJ/YBuTG/OXG/PWJHSUPUlUV2dnZHHaVFPOXG/MntgG5MX+5MX994tQ7IiIiIiKSAqfeeSkhBGpqajjsKinmLzfmT2wDcmP+cmP++sSOkgepqorc3FwOu0qK+cuN+RPbgNyYv9yYvz5x6h0REREREUmBU++8lBAClZWVHHaVFPOXG/MntgG5MX+5MX99YkfJg4QQKCws5IdEUsxfbsyf2Abkxvzlxvz1iVPviIiIiIhICpx656WEECgvL+e3CZJi/nJj/sQ2IDfmLzfmr0/sKHmQEAIlJSX8kEiK+cuN+RPbgNyYv9yYvz5x6h0REREREUmBU++8lBACpaWl/DZBUsxfbsyf2Abkxvzlxvz1iR0lDxJCoKKigh8SSTF/uTF/YhuQG/OXG/PXJ069IyIiIiIiKXDqnZdSVRXFxcVQVVXropAGmL/cmD+xDciN+cuN+esTO0oeVlNTo3URSEPMX27Mn9gG5Mb85cb89YdT74iIiIiISAqceuelVFVFUVERh10lxfzlxvyJbUBuzF9uzF+f2FHysPr6eq2LQBpi/nJj/sQ2IDfmLzfmrz+cekdERERERFLg1DsvpaoqCgsLOewqKeYvN+ZPbANyY/5yY/76xI4SERERERHRaTj1joiIiIiIpMCpd15KVVXk5eVx2FVSzF9uzJ/YBuTG/OXG/PWJHSUP8/Pz07oIpCHmLzfmT2wDcmP+cmP++sOpd0REREREJAVOvfNSqqoiNzeXw66SYv5yY/7ENiA35i835q9P7Ch5mMVi0boIpCHmLzfmT2wDcmP+cmP++sOpd0REREREJAVOvfNSqqoiJyeHw66SYv5yY/7ENiA35i83mfMvKCjAmjVr8NZbb2HNmjUoKCjQukhOM2ldAJkoioKQkBAoiqJ1UUgDzF9uzJ/YBuTG/OUmY/5paWlYsGAB1q5dC5vN5thvMpkwbtw4PP7440hOTtawhGfHqXdERERERNRuNm7ciLFjx8JmszXpJDUymUwwmUxYt24dRowY4dGyceqdl1JVFdnZ2VIOuxLzlx3zJ7YBuTF/ucmUf1paGsaOHQur1dpiJwkAbDYbrFYrxo4di7S0NA+X0HnsKHmQoijo2LGjVMOu9AfmLzfmT2wDcmP+cpMp/wULFsBms+Fsk9aEELDZbFi4cKGHSuY6Tr0jIiIiIqJzVlBQgK5du7Y6ktQSk8mE3NxcREREuLFkf+DUOy+lqiqOHDkixbArNcf85cb8iW1AbsxfbrLkv3nzZpc6SUDDNLzNmze7p0DniB0lD1IUBREREVIMu1JzzF9uzJ/YBuTG/OUmS/4VFRVt+r3y8vJ2Lkn74PLgHqQoCoKDg7UuBmmE+cuN+RPbgNyYv9xkyT8kJKRNv+etp8dwRMmD7HY7Dh8+DLvdrnVRSAPMX27Mn9gG5Mb85SZL/qmpqTCZmo7DhKIDLkcyhmEgLkcyQtG0U2QymZCamurBUjqPI0oeZDAYEBsbC4OB/VMZMX+5MX9iG5Ab85ebLPlHRkZi3LhxWLt2LWJtkZiARFyPMnRBJYywww4jTiAB/0EoPkQGck0FGD9+vMcWcnAVV70jjygtBdLSgJoawGIBkpOBsDCtS0VERERE7SktLQ3TB07GfJsJ3VCCMlhwEhbYYIAJKjqhBqGowVF0xJMmG97a/S6Sk5M9Vj5X+gYcUfIgu92OjIwMJCYmwmg0al0cj8jKAj78EFi/HigqAmw2wGQCOncGRo0CJkwA4uO1LqVnyJg//YH5E9uA3Ji/3GTKv2NZCF60dESHihwcRmeo+GMBCxuMKEAwTiAICTiJFy1x6FjWtvOaPIEjSh5SeOgQ9r33Huqqq+EfGIj+kyYholcvzcrjCbt3A3PmAEePAhH+J9Cz5gf426tQZwzCYcvlKKzrgm7dgEWLgIEDtS6te8mYPzUlhEBdXR38/f19ftWjlhzelY7vX12PurIq+IcG4cr7RqHnoCSti+VRMreBwkOHkPbee6grK4N/aCiSJTsGNrZ/m7UWpgCzdO1f9s+/bPlvvn4xun+9HHlB8aiqqUGttbbZY8wBZgRZLIiuykL2sOlIXf+Ix8rnSt+AHSU3O7BhA3Y/+ii67t2LTgCMAOwATgI4dsEFGPj88+h93XUeL5e7ZWUBM2YAtQcOYnjZa7iy4ht0RrljfmoROuD7kGvwVehMmHufh2XLfHNkqTH/2L370Rl+MEKBHQJFqEfuBf18Nn+iRltWf4Pdj76FlPzfT5ujHoxfo/pg4PPTMXTKNVoXk9xE1r+BjWRv/6y/fPUvyy5F5gVj4FdfjargKAAN15Cqq6uDEAKKosDf399xrlZwZT7q/IKQsO8zhHYL9UgZecFZL7HtlVeQe/31GLR3L8wAsk0mHLztNmSbTDADuHDvXuRefz22vfKK1kVtdx9+CBj3/oinjt2GCRVrYYEVR9EZ6YjBUXSGBVZMqPgYTx27Dca9P+Cjj7Qucfvb9sorOHb9Lbhw71FYkIhs00U4eNsDyDZdBAsSceHeozh2/S0+mT81Z7fbceDAAZ9f8ehUa+asRN3U2bg+fzsCUYejCEU6OuMoQhGIOlyfvx11U2djzZyVWhfVI2RrA6f/DTwKIP1///r630CgefvPNoXj4G03I9sULkX7l/3zL2v+Wf9OQ3BNEaotnRz7DAYDzGYzLBYLzGZzkwUtqiydEFxzAtlfpGlR3LPSRUfptddeQ0JCAsxmMy688EJs3bpV6yKd1YENG1Dz4IOIFgLpAArQcOVh5bPPYLPZUADgMIBoIVDz4IM4sGGDtgVuR6WlwNZ/7cdfTzyMWBThMKJRgI6wwQRAgQ0mFKAjDiMKsSjCX0/8Fd+/vR9lZVqXvP0c2LAB1Q8+gRjRCYeRjAJ0gc1WC+WzH2Cz1aIAXXAYyYgRnVD94BM+lT+1zGAwIDEx0edXPGq0ZfU3CF/8CmJRgsPohAIEwwYjGo4BDXPUD6MTYlGC8MWvYMvqb7QustvJ1Aaa/w00wYaOADrDho4ogMln/wYCrbR/G6B89jNsNvh8+5f98y9z/raKGhiEDcLg3DIIwmCCQdhRX1bt5pK1jdcfrdesWYNZs2bhiSeewK+//oorrrgCI0eOxNGjR7Uu2hntfvRRxAmBTADqqXfU1zs2VQCZALoKgV8ee8yzBXSjtDRgQPq76IZcZCISaivNTIUBmYhEN+TigvT3kOadXya0yU8PPYs4EYJMxENFDQArAAHU2xv+hRUqapCJeMSJDvjp4ee0LTB5hAz/QW60+9G30A3FyET4WY4B4eiGYux+7C0Pl1AbsrSBP/4GmqEiEcBFAPoDOP9//14EFYnIhNnn/gYCZ2j/9X+MJvpy+5f98y9z/qYQC1TFBEW1OfV4RbVBVYzwCw10c8naxuuP2EuWLMG0adMwffp09OnTB0uXLkVcXBxef/11rYvWqsJDh9B1716U4rROkskEMW5cw7Jv/6MCKAfQdc8enDh82KPldJcT6Vm4xvo1ymBp9QDZSIUB5bDgWusmnDic7aESulfhoUPoceAYytARKk45gdFkhBg3GDD9sdqNilqUIxQJv+f4TP7UMlVVkZ6eDlVVz/5gnTu8Kx0p+b+7dAxIyfsdR37N8FAJtSFLG/jjb6AFKi4AEIeGs5NqAFT9718jgDiouADlsPjU38BW23+LfwN8r/3L/vmXPf/4G5JRaemMwJqTTj0+qOYkKi1d0H2055YHd4VXd5Tq6urwyy+/YPjw4U32Dx8+HNu2bWvxd6xWK8rLy5vcADj+MKmq6tR24xoXrW3b7fZWt/e99x7CARSZTBBoGD8QJlPD2tgff4zG1TOEokCYTCgC0ElRsO/99xv2C+Eoy6nbzpbdHXUSQkAI0Wy7sYynbpf+uA5dUIwipQOEqaGJCYPS6vYJYxi6oBilP6zz2jqdut1SNqdu71jxIToZTTiJ4Ia6Gv9XV1UFPvkJsNkhjAYIQ8OqVyeMoehsMGHnig+9tk56aXveXCdFUZCUlOSogy/UqbWcvn91PTqjEkWm4IYyKoD4338Omm8bUAQLOhuqsOXVL7y2Tu2REwAkJSU5yuALdWqpjPveew/hRjNOoi+EwQJhqAZghTAq/zvuCQhj/f/2m3HCdD46GQKQ9t57XlsnV7a/f3U9OptqUGRo+IZcmP53vLfZgXXbIRrLYDJCKGho/6YabHnlC6+tkytt7/tX16PL/z7/4n+LOwqT8ZT/DzXdLoIFnZUqbHntj8+/t9XJlbb3/ev/q79iaXqsA6Cs3QZhtzf9P5ApCF1QiS2vfOG1dXKl7XXoFoqSy0chsL4cQqh//B/IoEAYlabbqh1meyVKLh+F0G6hHq2Ts7z6OkpFRUWw2+2IjIxssj8yMhL5+fkt/s7ChQvx/9u787io6vWB45+ZYR82EdkEBW/uK2Jxxf3mnqa3cktz1/ypJGKaSy6YilIp1xCuS2naVbstZi6VZKmZqVyVq1fJLRBccUE22Wbm/P4gSXCtdAY4z/v14uWZ7xyG5/h8z5nzzPc750RERNzV/sILL+Do6Eh0dDTz58/n6tWr+Pv7M3r0aKZPn05+fj5Dhw7Fzs6OuLg4rKysWLduHTNmzCAtLY169eoxaNAg3nzzTaysrOjUqRPVq1dnzZo15OTksHLlStasWcOPP/6Iw40bROh0zPHwgEuXaG1lhWurNmzd9QNobBn+t2fZ+80WTllb4eLlxbjUVMbb2JC3YQOjfXzw9fXlnXfewd7enunTp7NlyxaOHTuGra0t0dHRhIaGcuvWLdq2bUunTp148803S7Z9+/bt7N27FycnJ6Kjo5kwYQKFhYU0bdqUrl27smjRIvLy8pgyZQrJycn8+9//xsbGhu3bt9O3b1/y8/Np3bo1QUFBLF26FK1Wy0svvYTBYGDTpk3k5uayZcsWwsPDOX36NG3atGHQoEG8+uqrODo6MmLECE7kHeYjLpBjk83ENs+x4dvtXLZR8PP0otc5LTF2l8Dehk7XbVG83Pg27wqOGVcZrL9EaGgoaWlpeHl5MWXKFMLDw8nPz6dPnz74+fmxePFidDody5cvZ/HixZw8eRJ/f3/CwsKYOHEiVlZWtG3bloYNGxIXF0dubi7R0dFs3bqV+Ph4fH19Wbt2LZ06dcLBwYEePXrg4uLC+vXrURSF0aNHc/r0ab7//ns0Gg1btmzhxRdfJDMzk969exMYGMicOXNwcHBg4sSJ7N27l4SEBBRFITY2lmnTpnF87xFa+zpQ71wOH+qvg07LizedSKvuxMHsVCjMZVpgV2LP7icr4ya1q/vTV+fBrDVrWZLwPePHjyc7O5sPP/wQa2trPv30U1599VWuX79OUFAQnTt3ZuHCheh0Onr06IGjoyMbN24kJyeHDRs2EBkZydGjRwkMDOT1119n4MCBODo6MnDgQDIyMti+fTtGo5Hp06ezbds2EhMTcXd3Z9myZfTr14+ioiKGDx+Ora0tcXFx2NvbM3fuXD766CNOnTpFlSpVeOuttxg/fjwFBQV0796dwMBA5s2bh1arJTo6mrVr13LkyBE8PDyYN28eY8aMQafTERwcTMuWLYmOjubWrVvMnTuXgwcPsmXLFqpUqcKnn35Kz5490Wq1dOzYkZo1a/L++++j0WgYNGgQ169fZ/v27eTn5/Ptt9/yyiuvcPHiRbp06UKXLl2YNGkSer2esWPHcvToUfbu3YvBYGDp0qW8/fbbpKam0qBBA0JDQ/m///s/jEYjr7/+OikpKXz66afY2dk9sWPEihUr+OCDD9i/fz9169YlKiqKXr16odfrefHFF1EUhU2bNmEymQgLC2P//v389NNP6PV6NmzYQM+ePSksLGTAgAHl/hjxd5d2vKW7RrKHwl8u5tPR0Y/ljlfQXMqgm1NNcmt7sufQQRQbaya26cGGb7eRZVOE/tg2ul/pz9ChQzEYDIwbN47s7GzWrVuHnZ0dCxcuZNmyZRX2GLFs2TKmTp3KjRs3aNGiBS+//DLh4eGYTCZmzZpFQkICX3/9Nba2tsTGxjJ16lSysrKoV68e/fr1IyIigry8vHJ/jCi6eJHZzYKZc+h/KI62/NWlOt5p2WxyygajiWHZVdnva0PStYs45hYR1roPoQdPkvfhh/TU6Sr8MeLAiYOsdsmmIC+TCbe8+bS6iQtXLuOdr6FP+54s/X4zmoIi2jcLxvbUJb7JOod1NRd6Xr/MzJkz+fHHHyv0MeLMgSSe5Tp1/taKD77bhsZgpGebTlzbm8g+wzXwqMKUK+68TyrXPRz4y8U8hjs68d7uD9j0/E5GjBhBSkoKO3fuxGAwEBkZyerVqzl16hQ1atQgIiKiXB8jWtUIYi03+F6Xgb23O6+nubDIOoWiKnqa3dBQz82TjbpLaNIz6eten1NVrDh98iI3d66i44UXmDZtWsm+XFGPEbkZOXS2saGa7QU+LMpCd/MWoQ3a8lVOKmdSU3B1q0JYvTbM37WBgipV6NogD92ePWY7j9i8efMj1yLl+vLgFy9epHr16uzbt4+WLVuWtM+fP59169bx888/3/U7BQUFFBQUlDzOysrCz8+PjIwMXF1dS6pTrVb7wGWNRoNGo7nvstFoRKvV3nP5+7fewisiglQrKwwGa8AHrLwBG5SXnoEvDqHNz0HRXAXdFawN2fhpNFyJiODZmTNLql2tVltq+VFjfxLbdHtevclkKrWs0+lKPkm4vbxz7ly857zHOU19DLoiNAZT8adpWs09l600dtQ0HudyxGt0ePPNcrlNdy7fKzd3Lq8YvJSQ9TGkGR0p0lqDRoPGaMJkYwUvtkTzyT5QFFAUNCYFnc5EDSWLfQNDGblmXLncporS98rzNimKwpkzZ6hVqxbW1taVYpvul6cPR77HM2tiSLVyw2go/jQVnQ6NwXiPZS3WhiL8tFkkDBvPsFUTyuU2PY48GQyGUjecrAzbdK88bZ06D/93PiDN6EmRtrhdY1KKP1n+9bh357LOypoapoukvjGKrvOmlstt+j3La0Ys5ZmP4kg1OWE0aYtHD0wKaLUofVvBJz+iLTIWjzYYjVgrRvysckgYPJZh708ol9v0e/rehyPfI3hNDOes3DAYQaP8Oops+PX7OWWWrQ2F+GmySBgxnmErJ5TLbfo9fW/1yH8Q/EEM5zQuGHQ2Jcc6xdYGegfDZ/vQKMpv50BaEzUNGRwcHsqQlaHlcpv+SN87uSGR3NemUiXvPAUaPTl6dxSdFdqiIhwKbmBvyCbD3hf79xZSb0CgWbcpKysLV1fXin8fpcLCQhwcHPjkk0/4+9//XtI+YcIEEhMT2b1790NfwxL3UUo/dYqf69bFDnuu0ASwBwxAIb8OvgI2FA/o5eHFUW6RR4PTp6n21FNmifFJSj91iqS6TbHHjyv4AA+6kokDXlzkFmk0OH2sUmz/tn8lox/UFwcyuYLrQ9f34ia5uJC34VO69a/55AMU4gk785/TnH+6Dw4UcgXHh67vRQ652FDj8GfUCvyLGSIUT5Laj4Fq7/+y/ere/jud35vCmchPcP1xG455V9EqRkwaHTn21bjZ6jmemtYH39b+Zo+r0txHycbGhqCgIOLj40u1x8fHExISYqGoHs6jTh3O1muJC25osaH4y6sFKBoFxcUBRVN81TPIRYsNzlTll/ohlaJIgOLtv9C0Li6kouUWoAdsKS4Q+fVfW0CPlls4k8r5ZvUqzfa3ei6A3Q5tcaEA7R2X81A0/Jr/39bVYsKZAnY5tCOkW8U/QRD3pygKBQUFv2tudEX1VIvaHPGqjwt5pfaBeyneB/I44l2/0p0klKWWPmBdLYCdNl1wIfcR83+Lb226oHOrHMfA+/X/+78HVK7+r/b9X+35v5Nva3/ab5tMwNHN3FqynIyIpdxaspyAo5tpv22yRYqk36tcF0oA4eHhrFq1ig8++ICkpCQmTpxIamoqY8aMsXRoD5TWMppUahBA8m87ik4HHZsV/0vxDhJAMqn4cb7lEssF+wQ0X7iQ85pCAkhEyzmK78VuT3HRZA8Y0XKOABI5rykkKDLSovE+Tq6uoOnXj1SqE8CVh+T/CqlUR9e/Ly7muSG1sBCTycS5c+dKphNUds0XjiQVNwK4cd+TpeJ94AapuNE8cqSZIzQ/tfQBe3vY4/UKqXiWPgaW8dsx0IMfvF7BoXxeHfgPuWf/v+d7QOXs/2rf/9We/7JcarjQZGxrgqZ1psnY1rjUqDgnPOV66t1tsbGxREVFcenSJRo1asSSJUto27btI/2uJabe3bwJvXpBleTdjEsLowbpZOHANZwwoMMKI+5k40weqVRjmV80N2u1Y/NmKtXJ8r6YGPJeew0/RSETLddwvmP7s3DGxHmNBvulSwkZP97S4T5WKSkwvetOhp2cSQ1SycL+Pvmvweq6b7Hg62fx97d01EI8Xh+/sRq3qBhqcOPXfcAeA1qsMOFO3q/7gBs3poyn36Jhlg5XPCbyHlhM7f1ftl/d21+e/Z7aoEIUSn+GJQqlH36AMWPAzw9cMo7RKGkVbbK/w12ThVLVCc31bK4pzvzg9Df+V38kmVUak5YGy5dD69ZmCdFsfv76aw5Nm4ZvYiLuFN85wwhcA843a0ZQZCT1una1bJBPyOHD8NaI/9HwxEY6Fm7HXXMDpaojmus5XFPc+NamO8cb9Gfm+41o3tzS0YonTVEU8vPzsbOzQ6PRPPwXKondH+7k8LRVBF5Koho56DBiRMdVHDniXZ/mkSNpN+RZS4dpFmrqA1FR8P770MT5GE1OFr8HViMTHSaMaLmKCz84/Y2jdUdyNKsxI0fC5MmWjvrxu7P/u2ty73gP0Kui/6t9/1d7/ssrKZTuYIlCaccOCA2F2rXh9nuhLiedalf28ZfgAs4esOWqZwhGRw+g+AJoZ87A0qVQ5pZRlcbVM2c4tn49BRkZ2FapQuOXX64030l6kJQU+OQT2PHJVVzO/4eGbfM5vseOTN8WdO5TjT59kJEklTAajSQnJxMQEIBOp3v4L1Qyvxw5y+6YrRRk5GBbxZF243tUyjn5D6KmPpCSAqNGwfnzEBAANnnpVL24D+uiHIqsHbnuE0KhvQfJyeDrCytXVu5j4S9HzrI7ditWnrYYrhTQbqy6+r/a93+157+8kULpDpYeUfr1KsAPVFREpR1REsUyM+HYMbh1CxwcoHHjyjXFRAghyjp8GN54A1JTwdkZqlaF2/dev34dsrKgRg1YtAgZVRdCmM3vqQ3K9Q1nK6rGjcHdvfiNwMvrzmcUnJ1zycrS89sV4IrXq1at+PdE5eTiAq1aKeTm5qLX6yv9tBtxN0WR/Kud2vpA8+bFI0WffALbthV/IGg0Fn+XvVo16NsXVY2qqy3/ojTJf8UkhdIT4OoKzz1XPD+7WrWSC5yg1SpUr55OTo4/JlPxTmI0Fn+q1revjDBUdoqikJ6ejr+/vxwkVUjyL9TYB/z9i797NHq0jKqrMf/iN5L/ikmm3j0hZedn32s6utGIauZnCyGEEEIIYWmV5oazFZm/f/G8a19fOHsWLl8Gg0HB1TULg0Hh8uXidl/f4vWkSKr8FEUhKyur0t9sUtyb5F9IH1A3yb+6Sf4rJimUnqDb87NHjgS9Hi5cUIAMLlxQ0OuL21eulC+xqoWiKGRkZMhBUqUk/0L6gLpJ/tVN8l8xydQ7M5GrngkhhBBCCGFZMvWuHLp91bNnnrlJq1aKFEkqpCgKN2/elE+TVEryL6QPqJvkX90k/xWTFEpmpCgK2dnZspOolORf3ST/QvqAukn+1U3yXzHJ1DshhBBCCCGEKsjUu3LKZDJx48YNTCaTpUMRFiD5VzfJv5A+oG6Sf3WT/FdMUiiZWV5enqVDEBYk+Vc3yb+QPqBukn91k/xXPDL1TgghhBBCCKEKMvWunDKZTFy7dk2GXVVK8q9ukn8hfUDdJP/qJvmvmKRQMrOioiJLhyAsSPKvbpJ/IX1A3ST/6ib5r3hk6p0QQgghhBBCFWTqXTllMplIT0+XYVeVkvyrm+RfSB9QN8m/ukn+KyYplIQQQgghhBCiDJl6J4QQQgghhFCF31MbWJkpJou5XQdmZWVZOJLfhl09PDzQamUwT20k/+om+RfSB9RN8q9ukv/y43ZN8ChjRZW+UMrOzgbAz8/PwpEIIYQQQgghyoPs7GxcXFweuE6ln3pnMpm4ePEiTk5OaDQai8aSlZWFn58faWlpMg1QhST/6ib5F9IH1E3yr26S//JDURSys7Px8fF56OhepR9R0mq1+Pr6WjqMUpydnWUnUTHJv7pJ/oX0AXWT/Kub5L98eNhI0m0ySVIIIYQQQgghypBCSQghhBBCCCHKkELJjGxtbZk9eza2traWDkVYgORf3ST/QvqAukn+1U3yXzFV+os5CCGEEEIIIcTvJSNKQgghhBBCCFGGFEpCCCGEEEIIUYYUSkIIIYQQQghRhhRKQgghhBBCCFGGFEpmFBsbS0BAAHZ2dgQFBfHDDz9YOiRhBpGRkTz99NM4OTnh4eFB7969OXnypKXDEhYSGRmJRqMhLCzM0qEIM7lw4QKDBg2iatWqODg40KxZMw4dOmTpsIQZGAwG3nzzTQICArC3t6dWrVrMnTsXk8lk6dDEE7Bnzx569uyJj48PGo2GL774otTziqIwZ84cfHx8sLe3p3379hw/ftwywYpHIoWSmXz88ceEhYUxY8YMjhw5Qps2bejWrRupqamWDk08Ybt372bcuHHs37+f+Ph4DAYDnTt3Jjc319KhCTNLSEhgxYoVNGnSxNKhCDPJyMigVatWWFtb89VXX3HixAneffddXF1dLR2aMINFixbxz3/+k5iYGJKSkoiKiuLtt9/mvffes3Ro4gnIzc2ladOmxMTE3PP5qKgoFi9eTExMDAkJCXh5edGpUyeys7PNHKl4VHJ5cDMJDg6mefPmxMXFlbTVr1+f3r17ExkZacHIhLldvXoVDw8Pdu/eTdu2bS0djjCTnJwcmjdvTmxsLPPmzaNZs2ZER0dbOizxhE2dOpUff/xRZhCoVI8ePfD09OT9998vaXvxxRdxcHBg3bp1FoxMPGkajYZNmzbRu3dvoHg0ycfHh7CwMN544w0ACgoK8PT0ZNGiRbz66qsWjFbcj4womUFhYSGHDh2ic+fOpdo7d+7Mvn37LBSVsJTMzEwA3NzcLByJMKdx48bx3HPP0bFjR0uHIszoyy+/pEWLFvTp0wcPDw8CAwNZuXKlpcMSZtK6dWt27tzJqVOnAPjvf//L3r176d69u4UjE+aWnJzM5cuXS50L2tra0q5dOzkXLMesLB2AGly7dg2j0Yinp2epdk9PTy5fvmyhqIQlKIpCeHg4rVu3plGjRpYOR5jJxo0bOXz4MAkJCZYORZjZL7/8QlxcHOHh4UyfPp2DBw/y2muvYWtry+DBgy0dnnjC3njjDTIzM6lXrx46nQ6j0cj8+fMZMGCApUMTZnb7fO9e54Lnzp2zREjiEUihZEYajabUY0VR7moTldv48eM5evQoe/futXQowkzS0tKYMGECO3bswM7OztLhCDMzmUy0aNGCBQsWABAYGMjx48eJi4uTQkkFPv74Yz766CPWr19Pw4YNSUxMJCwsDB8fH4YMGWLp8IQFyLlgxSKFkhm4u7uj0+nuGj1KT0+/65MFUXmFhoby5ZdfsmfPHnx9fS0djjCTQ4cOkZ6eTlBQUEmb0Whkz549xMTEUFBQgE6ns2CE4kny9vamQYMGpdrq16/PZ599ZqGIhDlNnjyZqVOn0r9/fwAaN27MuXPniIyMlEJJZby8vIDikSVvb++SdjkXLN/kO0pmYGNjQ1BQEPHx8aXa4+PjCQkJsVBUwlwURWH8+PF8/vnnfPfddwQEBFg6JGFGzz77LMeOHSMxMbHkp0WLFgwcOJDExEQpkiq5Vq1a3XU7gFOnTlGzZk0LRSTM6datW2i1pU+1dDqdXB5chQICAvDy8ip1LlhYWMju3bvlXLAckxElMwkPD+eVV16hRYsWtGzZkhUrVpCamsqYMWMsHZp4wsaNG8f69evZvHkzTk5OJSOLLi4u2NvbWzg68aQ5OTnd9X00vV5P1apV5XtqKjBx4kRCQkJYsGABffv25eDBg6xYsYIVK1ZYOjRhBj179mT+/PnUqFGDhg0bcuTIERYvXszw4cMtHZp4AnJycjhz5kzJ4+TkZBITE3Fzc6NGjRqEhYWxYMECateuTe3atVmwYAEODg68/PLLFoxaPIhcHtyMYmNjiYqK4tKlSzRq1IglS5bI5aFV4H5zj1evXs3QoUPNG4woF9q3by+XB1eRrVu3Mm3aNE6fPk1AQADh4eGMGjXK0mEJM8jOzmbmzJls2rSJ9PR0fHx8GDBgALNmzcLGxsbS4YnHbNeuXXTo0OGu9iFDhrBmzRoURSEiIoLly5eTkZFBcHAwy5Ytkw/NyjEplIQQQgghhBCiDPmOkhBCCCGEEEKUIYWSEEIIIYQQQpQhhZIQQgghhBBClCGFkhBCCCGEEEKUIYWSEEIIIYQQQpQhhZIQQgghhBBClCGFkhBCCCGEEEKUIYWSEEIIIYQQQpQhhZIQQgizunz5Mp06dUKv1+Pq6nrfNo1GwxdffPFIrzlnzhyaNWv2ROI1l/bt2xMWFmbpMIQQQvxKCiUhhBBAcbESGhpKrVq1sLW1xc/Pj549e7Jz587H+neWLFnCpUuXSExM5NSpU/dtu3TpEt26dXuk13z99dcfe5xr1qwpKdoex3pCCCEqFitLByCEEMLyUlJSaNWqFa6urkRFRdGkSROKior45ptvGDduHD///PNj+1tnz54lKCiI2rVrP7DNy8vrkV/T0dERR0fHxxajEEIIISNKQgghGDt2LBqNhoMHD/LSSy9Rp04dGjZsSHh4OPv37y9ZLzU1lV69euHo6IizszN9+/blypUrpV5ry5YtBAUFYWdnR61atYiIiMBgMADg7+/PZ599xtq1a9FoNAwdOvSebXD31Lvz58/Tv39/3Nzc0Ov1tGjRggMHDgD3nnq3evVq6tevj52dHfXq1SM2NrbkuZSUFDQaDZ9//jkdOnTAwcGBpk2b8tNPPwGwa9cuhg0bRmZmJhqNBo1Gw5w5cx7p//J2LOvWrcPf3x8XFxf69+9PdnZ2yTq5ubkMHjwYR0dHvL29effdd+96ncLCQqZMmUL16tXR6/UEBweza9cuAPLz82nYsCGjR48uWT85ORkXFxdWrlz5SHEKIYR4MBlREkIIlbtx4wZff/018+fPR6/X3/X87WlliqLQu3dv9Ho9u3fvxmAwMHbsWPr161dyAv/NN98waNAgli5dSps2bTh79mzJyfzs2bNJSEhg8ODBODs7849//AN7e3sKCwvvaisrJyeHdu3aUb16db788ku8vLw4fPgwJpPpntu0cuVKZs+eTUxMDIGBgRw5coRRo0ah1+sZMmRIyXozZszgnXfeoXbt2syYMYMBAwZw5swZQkJCiI6OZtasWZw8eRLgd41YnT17li+++IKtW7eSkZFB3759WbhwIfPnzwdg8uTJfP/992zatAkvLy+mT5/OoUOHShV7w4YNIyUlhY0bN+Lj48OmTZvo2rUrx44do3bt2vzrX/8iODiY7t2707NnT1555RU6dOjAqFGjHjlOIYQQD6AIIYRQtQMHDiiA8vnnnz9wvR07dig6nU5JTU0taTt+/LgCKAcPHlQURVHatGmjLFiwoNTvrVu3TvH29i553KtXL2XIkCGl1rlXG6Bs2rRJURRFWb58ueLk5KRcv379nrHNnj1badq0acljPz8/Zf369aXWeeutt5SWLVsqiqIoycnJCqCsWrXqrm1JSkpSFEVRVq9erbi4uNz7P+MOZdebPXu24uDgoGRlZZW0TZ48WQkODlYURVGys7MVGxsbZePGjSXPX79+XbG3t1cmTJigKIqinDlzRtFoNMqFCxdK/a1nn31WmTZtWsnjqKgoxd3dXQkNDVW8vLyUq1evPjReIYQQj0ZGlIQQQuUURQGKp7o9SFJSEn5+fvj5+ZW0NWjQAFdXV5KSknj66ac5dOgQCQkJJSMnAEajkfz8fG7duoWDg8MfijExMZHAwEDc3Nweuu7Vq1dJS0tjxIgRpUZXDAYDLi4updZt0qRJybK3tzcA6enp1KtX7w/FeZu/vz9OTk6lXjs9PR0oHm0qLCykZcuWJc+7ublRt27dkseHDx9GURTq1KlT6nULCgqoWrVqyeNJkyaxefNm3nvvPb766ivc3d3/VNxCCCF+I4WSEEKoXO3atdFoNCQlJdG7d+/7rqcoyj2LqTvbTSYTERERvPDCC3etZ2dn94djvNd0vPu5PR1v5cqVBAcHl3pOp9OVemxtbV2yfOc2/Fl3vu7t1779urcL0wcxmUzodDoOHTp0V8x3TgFMT0/n5MmT6HQ6Tp8+TdeuXf907EIIIYrJxRyEEELl3Nzc6NKlC8uWLSM3N/eu52/evAkUjx6lpqaSlpZW8tyJEyfIzMykfv36ADRv3pyTJ0/y1FNP3fWj1f7xt5wmTZqQmJjIjRs3Hrqup6cn1atX55dffrkrhoCAgEf+mzY2NhiNxj8c8/089dRTWFtbl7pIRkZGRsll0QECAwMxGo2kp6fftQ13Xg1w+PDhNGrUiLVr1zJlyhROnDjx2OMVQgi1khElIYQQxMbGEhISwjPPPMPcuXNp0qQJBoOB+Ph44uLiSEpKomPHjjRp0oSBAwcSHR1dcjGHdu3a0aJFCwBmzZpFjx498PPzo0+fPmi1Wo4ePcqxY8eYN2/eH45vwIABLFiwgN69exMZGYm3tzdHjhzBx8en1BS22+bMmcNrr72Gs7Mz3bp1o6CggP/85z9kZGQQHh7+SH/T39+fnJwcdu7cSdOmTXFwcPjDUwfv5OjoyIgRI5g8eTJVq1bF09OTGTNmlCok69Spw8CBAxk8eDDvvvsugYGBXLt2je+++47GjRvTvXt3li1bxk8//cTRo0fx8/Pjq6++YuDAgRw4cAAbG5s/HacQQqidjCgJIYQgICCAw4cP06FDByZNmkSjRo3o1KkTO3fuJC4uDvjtct1VqlShbdu2dOzYkVq1avHxxx+XvE6XLl3YunUr8fHxPP300/z1r39l8eLF1KxZ80/FZ2Njw44dO/Dw8KB79+40btyYhQsX3jUt7baRI0eyatUq1qxZQ+PGjWnXrh1r1qz5XSNKISEhjBkzhn79+lGtWjWioqL+1Dbc6e2336Zt27Y8//zzdOzYkdatWxMUFFRqndWrVzN48GAmTZpE3bp1ef755zlw4AB+fn78/PPPTJ48mdjY2JLvjC1btoybN28yc+bMxxanEEKomUZ5lMnSQgghhBBCCKEiMqIkhBBCCCGEEGVIoSSEEEIIIYQQZUihJIQQQgghhBBlSKEkhBBCCCGEEGVIoSSEEEIIIYQQZUihJIQQQgghhBBlSKEkhBBCCCGEEGVIoSSEEEIIIYQQZUihJIQQQgghhBBlSKEkhBBCCCGEEGVIoSSEEEIIIYQQZfw/AduULCp0BasAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "idx = np.arange(len(beta_true))\n", + "plt.figure(figsize=(10, 6))\n", + "\n", + "# True Parameters\n", + "plt.scatter(idx, beta_true, \n", + " color='black', label='True', s=80, marker='o')\n", + "\n", + "# Elastic Net Estimates\n", + "plt.scatter(idx, beta_EN / scaler.scale_, \n", + " color='blue', label='Elastic Net', s=80, alpha=0.7)\n", + "\n", + "# Adaptive Elastic Net Estimates\n", + "plt.scatter(idx, beta_ENAL / scaler.scale_, \n", + " color='red', label='Adaptive Elastic Net', s=80, alpha=0.7)\n", + "\n", + "plt.axhline(0, color='k', linestyle='--', linewidth=0.5)\n", + "plt.xlabel('Coefficient Index')\n", + "plt.ylabel('Coefficient Value')\n", + "plt.legend()\n", + "plt.title('Adaptive Elastic Net vs Elastic Net: Coefficient Estimation')\n", + "plt.grid(True, linestyle=':', alpha=0.5)\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ReHLine", + "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.10.18" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/rehline/_base.py b/rehline/_base.py index fe6ba31..d9c7bd1 100644 --- a/rehline/_base.py +++ b/rehline/_base.py @@ -195,7 +195,7 @@ def call_ReLHLoss(self, score): if self.H > 0: rehu_input = (self._S.T * score[:, np.newaxis]).T + self._T return np.sum(_relu(relu_input), 0) + np.sum(_rehu(rehu_input, self._Tau), 0) - + @abstractmethod def fit(self, X, y, sample_weight): """Fit model.""" @@ -286,7 +286,7 @@ def ReHLine_solver( T=None, A=None, b=None, - rho=0.0, + rho=None, Lambda=None, Gamma=None, xi=None, @@ -307,6 +307,8 @@ def ReHLine_solver( A = np.empty(shape=(0, 0)) if b is None: b = np.empty(shape=(0)) + if rho is None: + rho = np.empty(shape=(0)) if Lambda is None: Lambda = np.empty(shape=(0, 0)) if Gamma is None: @@ -330,6 +332,7 @@ def ReHLine_solver( X, A, b, + rho, U, V, S, @@ -337,7 +340,6 @@ def ReHLine_solver( Tau, max_iter, tol, - rho, shrink, verbose, trace_freq, @@ -773,32 +775,32 @@ def _cast_sample_weight(U, V, Tau, S, T, C=1.0, sample_weight=None): # U_new = np.zeros((self.L+2, n+d)) # V_new = np.zeros((self.L+2, n+d)) # ## Block 1 -# if len(self.U): -# U_new[:self.L, :n] = self.U -# V_new[:self.L, :n] = self.V +# if len(self._U): +# U_new[:self.L, :n] = self._U +# V_new[:self.L, :n] = self._V # ## Block 2 # U_new[-2,n:] = l1_pen # U_new[-1,n:] = -l1_pen -# if len(self.S): +# if len(self._S): # S_new = np.zeros((self.H, n+d)) # T_new = np.zeros((self.H, n+d)) # Tau_new = np.zeros((self.H, n+d)) -# S_new[:,:n] = self.S -# T_new[:,:n] = self.T -# Tau_new[:,:n] = self.Tau +# S_new[:,:n] = self._S +# T_new[:,:n] = self._T +# Tau_new[:,:n] = self._Tau -# self.S = S_new -# self.T = T_new -# self.Tau = Tau_new +# self._S = S_new +# self._T = T_new +# self._Tau = Tau_new # ## fake X # X_fake = np.zeros((n+d, d)) # X_fake[:n,:] = X # X_fake[n:,:] = np.identity(d) -# self.U = U_new -# self.V = V_new +# self._U = U_new +# self._V = V_new # self.auto_shape() # return X_fake diff --git a/rehline/_class.py b/rehline/_class.py index 0da8ac0..d61cf1a 100644 --- a/rehline/_class.py +++ b/rehline/_class.py @@ -493,7 +493,7 @@ class plqERM_ElasticNet(_BaseReHLine, BaseEstimator): .. math:: - \min_{\mathbf{\beta} \in \mathbb{R}^d} C \sum_{i=1}^n \text{PLQ}(y_i, \mathbf{x}_i^T \mathbf{\beta}) + \text{l1_ratio} \| \mathbf{\beta} \|_1 + \frac{1}{2} (1 - \text{l1_ratio}) \| \mathbf{\beta} \|_2^2, \ \text{ s.t. } \ + \min_{\mathbf{\beta} \in \mathbb{R}^d} C \sum_{i=1}^n \text{PLQ}(y_i, \mathbf{x}_i^T \mathbf{\beta}) + \text{l1_ratio} \sum_{j=1}^d \omega_j | \beta_j | + \frac{1}{2} (1 - \text{l1_ratio}) \| \mathbf{\beta} \|_2^2, \ \text{ s.t. } \ \mathbf{A} \mathbf{\beta} + \mathbf{b} \geq \mathbf{0}, The function supports various loss functions, including: @@ -527,6 +527,9 @@ class plqERM_ElasticNet(_BaseReHLine, BaseEstimator): The ElasticNet mixing parameter, with 0 <= l1_ratio < 1. For l1_ratio = 0 the penalty is an L2 penalty. For 0 < l1_ratio < 1, the penalty is a combination of L1 and L2. + omega : array of shape (n_features, ), default=np.empty(shape=(0, 0)) + Weight coefficients for adaptive lasso. + verbose : int, default=0 Enable verbose output. Note that this setting takes advantage of a per-process runtime setting in liblinear that, if enabled, may not work @@ -584,6 +587,7 @@ def __init__( constraint=None, C=1.0, l1_ratio=0.5, + omega=None, U=None, V=None, Tau=None, @@ -602,8 +606,8 @@ def __init__( self.constraint = constraint if constraint is not None else [] self.C = C self.l1_ratio = l1_ratio - self.rho = l1_ratio / (1 - l1_ratio) self.C_eff = C / (1 - l1_ratio) + self.omega = omega if omega is not None else np.empty(shape=(0, 0)) self._U = U if U is not None else np.empty(shape=(0, 0)) self._V = V if V is not None else np.empty(shape=(0, 0)) self._S = S if S is not None else np.empty(shape=(0, 0)) @@ -678,6 +682,25 @@ def fit(self, X, y, sample_weight=None): self._xi = np.empty(shape=(0, 0)) self._mu = np.empty(shape=(0, 0)) + if self.l1_ratio == 0: + self.rho = None + if self.omega.size > 0: + warnings.warn( + f"Omega will be ignored since l1_ratio=0.", + UserWarning, + stacklevel=2 + ) + else: + if self.omega.size not in (0, d): + raise ValueError( + f"Omega length {self.omega.size} must be 0 or {d} (n_features)" + ) + if not np.all(self.omega > 0): + raise ValueError( + "All elements in omega must be strictly positive." + ) + self.rho = np.full(d, self.l1_ratio / (1 - self.l1_ratio)) * (self.omega if self.omega.size == d else 1.0) + result = ReHLine_solver( X=X, U=U_weight, diff --git a/rehline/_internal.pyi b/rehline/_internal.pyi index 5e03550..2486f8d 100644 --- a/rehline/_internal.pyi +++ b/rehline/_internal.pyi @@ -19,6 +19,7 @@ def rehline_internal( X: npt.NDArray[np.float64], A: npt.NDArray[np.float64], b: npt.NDArray[np.float64], + rho: npt.NDArray[np.float64], U: npt.NDArray[np.float64], V: npt.NDArray[np.float64], S: npt.NDArray[np.float64], @@ -26,7 +27,6 @@ def rehline_internal( Tau: npt.NDArray[np.float64], max_iter: int, tol: float, - rho: float = ..., shrink: int = ..., verbose: int = ..., trace_freq: int = ..., diff --git a/src/rehline.cpp b/src/rehline.cpp index 5754f7c..0fc5aba 100644 --- a/src/rehline.cpp +++ b/src/rehline.cpp @@ -19,15 +19,15 @@ using ReHLineResult = rehline::ReHLineResult; void rehline_internal( ReHLineResult& result, - const MapMat& X, const MapMat& A, const MapVec& b, + const MapMat& X, const MapMat& A, const MapVec& b, const MapVec& rho, const MapMat& U, const MapMat& V, const MapMat& S, const MapMat& T, const MapMat& Tau, - int max_iter, double tol, double rho = 0.0, int shrink = 1, + int max_iter, double tol, int shrink = 1, int verbose = 0, int trace_freq = 100 ) { - rehline::rehline_solver(result, X, A, b, U, V, S, T, Tau, - max_iter, tol, rho, shrink, verbose, trace_freq); + rehline::rehline_solver(result, X, A, b, rho, U, V, S, T, Tau, + max_iter, tol, shrink, verbose, trace_freq); } PYBIND11_MODULE(_internal, m) { diff --git a/src/rehline.h b/src/rehline.h index f388978..981ee84 100644 --- a/src/rehline.h +++ b/src/rehline.h @@ -80,6 +80,7 @@ void reset_fv_set(std::vector>& fvset, std::size_t n, st // * S, T, Tau: [H x n] // * A : [K x d] // * b : [K] +// * rho : [d] // - Pre-computed // * r: [n] // * p: [K] @@ -139,6 +140,7 @@ class ReHLineSolver const Index m_L; const Index m_H; const Index m_K; + const Index m_W; // Input matrices and vectors RMatrix m_X; @@ -149,7 +151,7 @@ class ReHLineSolver ConstRefMat m_Tau; RMatrix m_A; ConstRefVec m_b; - Scalar m_rho; // l1_ratio / (1 - l1_ratio) + ConstRefVec m_rho; // Pre-computed Vector m_gk_denom; // ||a[k]||^2 @@ -174,7 +176,7 @@ class ReHLineSolver // =================== Initialization functions =================== // // Compute the primal variable beta from dual variables - // beta = A'xi - U3 * vec(Lambda) - S3 * vec(Gamma) + // beta = A'xi - U3 * vec(Lambda) - S3 * vec(Gamma) + 2 * Mu - rho // A can be empty, one of U and V may be empty inline void set_primal() { @@ -193,9 +195,9 @@ class ReHLineSolver if (m_H > 0) LHterm.noalias() += m_S.cwiseProduct(m_Gamma).colwise().sum().transpose(); m_beta.noalias() -= m_X.transpose() * LHterm; - - if (m_rho > 0) - m_beta.noalias() += Scalar(2.0) * m_mu - m_rho * Vector::Ones(m_d); + // L1 related + if (m_W > 0) + m_beta.noalias() += Scalar(2.0) * m_mu - m_rho; } // =================== Evaluating objection function =================== // @@ -224,7 +226,7 @@ class ReHLineSolver // Quadratic term result += Scalar(0.5) * m_beta.squaredNorm(); // L1 penalty term - result += m_rho * m_beta.template lpNorm<1>(); + result += m_beta.cwiseAbs().cwiseProduct(m_rho).sum(); return result; } @@ -251,8 +253,8 @@ class ReHLineSolver } // 2 * Mu - rho, [d x 1] Vector MuR = Vector::Zero(m_d); - if (m_rho > 0) - MuR = Scalar(2.0) * m_mu - m_rho * Vector::Ones(m_d); + if (m_W > 0) + MuR = Scalar(2.0) * m_mu - m_rho; // Compute dual objective function value Scalar obj = Scalar(0); // If K = 0, all terms that depend on A, xi, or b will be zero @@ -261,7 +263,7 @@ class ReHLineSolver // 0.5 * ||Atxi||^2 - Atxi' * U3L - Atxi' * S3G + Atxi' MuR + xi' * b const Scalar Atxi_U3L = (m_L > 0) ? (Atxi.dot(U3L)) : Scalar(0); const Scalar Atxi_S3G = (m_H > 0) ? (Atxi.dot(S3G)) : Scalar(0); - const Scalar Atxi_MuR = (m_rho > 0) ? (Atxi.dot(MuR)) : Scalar(0); + const Scalar Atxi_MuR = (m_W > 0) ? (Atxi.dot(MuR)) : Scalar(0); obj += Scalar(0.5) * Atxi.squaredNorm() - Atxi_U3L - Atxi_S3G + Atxi_MuR + m_xi.dot(m_b); } // If L = 0, all terms that depend on U, V, or Lambda will be zero @@ -269,7 +271,7 @@ class ReHLineSolver { // 0.5 * ||U3L||^2 + U3L' * S3G - U3L' * MuR - tr(Lambda * V') const Scalar U3L_S3G = (m_H > 0) ? (U3L.dot(S3G)) : Scalar(0); - const Scalar U3L_MuR = (m_rho > 0) ? (U3L.dot(MuR)) : Scalar(0); + const Scalar U3L_MuR = (m_W > 0) ? (U3L.dot(MuR)) : Scalar(0); obj += Scalar(0.5) * U3L.squaredNorm() + U3L_S3G - U3L_MuR - m_Lambda.cwiseProduct(m_V).sum(); } @@ -277,14 +279,15 @@ class ReHLineSolver if (m_H > 0) { // 0.5 * ||S3G||^2 - S3G' * MuR + 0.5 * ||Gamma||^2 - tr(Gamma * T') - const Scalar S3G_MuR = (m_rho > 0) ? (S3G.dot(MuR)) : Scalar(0); + const Scalar S3G_MuR = (m_W > 0) ? (S3G.dot(MuR)) : Scalar(0); obj += Scalar(0.5) * S3G.squaredNorm() - S3G_MuR + Scalar(0.5) * m_Gamma.squaredNorm() - m_Gamma.cwiseProduct(m_T).sum(); } - // If rho = 0, all terms that depend on rho, or Mu will be zero - if (m_rho > 0) - obj += Scalar(2.0) * m_mu.squaredNorm() - Scalar(2.0) * m_rho * m_mu.sum() + - Scalar(0.5) * m_d * m_rho * m_rho; + // If W = 0, all terms that depend on rho or Mu will be zero + if (m_W > 0) + // 2.0 * ||Mu||^2 - 2.0 * rho' * Mu + 0.5 * ||rho||^2 + obj += Scalar(2.0) * m_mu.squaredNorm() - Scalar(2.0) * m_rho.dot(m_mu) + + Scalar(0.5) * m_rho.squaredNorm(); return obj; } @@ -368,7 +371,7 @@ class ReHLineSolver // Update mu and beta inline void update_mu_beta() { - if (m_rho <= 0) + if (m_W <= 0) return; // Save original Mu @@ -598,7 +601,7 @@ class ReHLineSolver // Overloaded version based on free variable set inline void update_mu_beta(std::vector& fv_set, Scalar& min_pg, Scalar& max_pg) { - if (m_rho <= 0) + if (m_W <= 0) return; // Permutation @@ -618,12 +621,13 @@ class ReHLineSolver for (auto j: fv_set) { const Scalar mu_j = m_mu[j]; + const Scalar rho_j = m_rho[j]; // Compute g_j const Scalar g_j = m_beta[j]; // PG and shrink Scalar pg; - const bool shrink = pg_mu(mu_j, g_j, m_rho, lb, ub, pg); + const bool shrink = pg_mu(mu_j, g_j, rho_j, lb, ub, pg); if (shrink) continue; @@ -632,7 +636,7 @@ class ReHLineSolver min_pg = std::min(min_pg, pg); // Compute new mu_j const Scalar candid = mu_j - g_j * Scalar(0.5); - const Scalar newmu = std::max(Scalar(0), std::min(m_rho, candid)); + const Scalar newmu = std::max(Scalar(0), std::min(rho_j, candid)); // Update mu and beta m_mu[j] = newmu; m_beta[j] += Scalar(2.0) * (newmu - mu_j); @@ -647,8 +651,9 @@ class ReHLineSolver ReHLineSolver(ConstRefMat X, ConstRefMat U, ConstRefMat V, ConstRefMat S, ConstRefMat T, ConstRefMat Tau, ConstRefMat A, ConstRefVec b, - Scalar rho) : - m_n(X.rows()), m_d(X.cols()), m_L(U.rows()), m_H(S.rows()), m_K(A.rows()), + ConstRefVec rho) : + m_n(X.rows()), m_d(X.cols()), m_L(U.rows()), m_H(S.rows()), m_K(A.rows()), + m_W(rho.rows()), // check if l1 penalty is implemented m_X(X), m_U(U), m_V(V), m_S(S), m_T(T), m_Tau(Tau), m_A(A), m_b(b), m_rho(rho), m_gk_denom(m_K), m_gli_denom(m_L, m_n), m_ghi_denom(m_H, m_n), @@ -691,10 +696,10 @@ class ReHLineSolver // Gamma.fill(std::min(1.0, 0.5 * tau)); } - // Each element of Mu satisfies 0 <= mu_j <= rho, + // Each element of Mu satisfies 0 <= mu_j <= rho_j, // and we use min(0.5 * rho, 1) to initialize (rho can be Inf) - if (m_rho > 0) - m_mu.setConstant( std::min(Scalar(0.5) * m_rho, Scalar(1)) ); + if (m_W > 0) + m_mu.noalias() = (m_rho * Scalar(0.5)).cwiseMin(Scalar(1)); // Set primal variable based on duals set_primal(); @@ -746,15 +751,15 @@ class ReHLineSolver } - if (m_rho > 0) + if (m_W > 0) { // Check shape of warmstart parameters if (mu_ws.size() != m_d){ throw std::invalid_argument("mu_ws must have size d"); } // Check values of warmstart parameters - if ((mu_ws.array() < Scalar(0)).any() || (mu_ws.array() > m_rho).any()){ - throw std::invalid_argument("mu_ws must be in [0, rho]"); + if ((mu_ws.array() < Scalar(0)).any() || (mu_ws.array() > m_rho.array()).any()){ + throw std::invalid_argument("mu_ws_j must be in [0, rho_j]"); } m_mu = mu_ws; } @@ -928,10 +933,10 @@ template void rehline_solver( ReHLineResult& result, const Eigen::MatrixBase& X, const Eigen::MatrixBase& A, - const Eigen::MatrixBase& b, + const Eigen::MatrixBase& b, const Eigen::MatrixBase& rho, const Eigen::MatrixBase& U, const Eigen::MatrixBase& V, const Eigen::MatrixBase& S, const Eigen::MatrixBase& T, const Eigen::MatrixBase& Tau, - Index max_iter, typename DerivedMat::Scalar tol, typename DerivedMat::Scalar rho = 0, + Index max_iter, typename DerivedMat::Scalar tol, Index shrink = 1, Index verbose = 0, Index trace_freq = 100, std::ostream& cout = std::cout ) diff --git a/tests/test_elastic_net.py b/tests/test_elastic_net.py index eace518..bd25377 100644 --- a/tests/test_elastic_net.py +++ b/tests/test_elastic_net.py @@ -31,7 +31,7 @@ # --------------------------------------------------------------------------- # Helper # --------------------------------------------------------------------------- - +import pytest def _regression_dataset(n, n_features, n_informative, seed=42): X, y = make_regression( @@ -223,3 +223,121 @@ def test_dual_variable_mu(): rho = clf.rho assert np.all(clf._mu >= 0), "mu should be non-negative" assert np.all(clf._mu <= rho + 1e-10), f"mu should be <= rho={rho}" + + +def test_different_omegas(): + """ElasticNet should fit successfully for all tested omega values.""" + n, n_features, C, l1_ratio = 2000, 10, 0.01, 0.5 + + X, y = make_regression( + n_samples=n, + n_features=n_features, + noise=0.1, + random_state=42, + n_informative=6, + ) + scaler = StandardScaler() + X_scaled = scaler.fit_transform(X) + + for i in range(5): # conduct 5 tests with different omega + rng = np.random.default_rng(seed=42+i) + omega = rng.uniform(low=0.1, high=0.2+i, size=n_features) + clf = plqERM_ElasticNet( + loss={"name": "mse"}, + C=C, + l1_ratio=l1_ratio, + omega=omega, + max_iter=5000, + tol=1e-4, + ) + clf.fit(X_scaled, y) + assert clf.coef_ is not None, f"Fit failed for omega={omega}" + assert clf.coef_.shape == (n_features,), f"Wrong coef_ shape for omega={omega}: {clf.coef_.shape}" + + +def test_with_omega_vs_without_omega(): + """ElasticNet with omega=(1, 1, ..., 1) should exactly match that without omega.""" + n, n_features, C, l1_ratio = 2000, 10, 0.01, 0.5 + + X, y = make_regression( + n_samples=n, + n_features=n_features, + noise=0.1, + random_state=42, + n_informative=6, + ) + scaler = StandardScaler() + X_scaled = scaler.fit_transform(X) + + clf_with_omg = plqERM_ElasticNet( + loss={"name": "mse"}, + C=C, + l1_ratio=l1_ratio, + omega=np.ones(n_features), + max_iter=5000, + tol=1e-4, + ) + clf_with_omg.fit(X_scaled, y) + + clf_without_omg = plqERM_ElasticNet( + loss={"name": "mse"}, + C=C, + l1_ratio=l1_ratio, + max_iter=5000, + tol=1e-4, + ) + clf_without_omg.fit(X_scaled, y) + + assert np.array_equal(clf_with_omg.coef_.flatten(), clf_without_omg.coef_.flatten()), \ + "ElasticNet with omega=(1, 1, ..., 1) should exactly match that without omega." + + +def test_omega_validation(): + """Test omega related validations raise appropriate warnings or errors""" + n, n_features, C, l1_ratio = 2000, 10, 0.01, 0.5 + + X, y = make_regression( + n_samples=n, + n_features=n_features, + noise=0.1, + random_state=42, + n_informative=6, + ) + scaler = StandardScaler() + X_scaled = scaler.fit_transform(X) + + # Test invalid omega shape (must align with n_features) + with pytest.raises(ValueError, match="Omega length"): + clf = plqERM_ElasticNet( + loss={"name": "mse"}, + C=C, + l1_ratio=l1_ratio, + omega=np.ones(n_features + 1), + max_iter=5000, + tol=1e-4, + ) + clf.fit(X_scaled, y) + # Test invalid omega value (all elements must be strictly positive) + with pytest.raises(ValueError, match="All elements in omega must be strictly positive"): + omega = np.ones(n_features) + omega[0] = -1 + clf = plqERM_ElasticNet( + loss={"name": "mse"}, + C=C, + l1_ratio=l1_ratio, + omega=omega, + max_iter=5000, + tol=1e-4, + ) + clf.fit(X_scaled, y) + # Test ineffective omega (when omega provided but l1_ratio==0) + with pytest.warns(UserWarning, match="Omega will be ignored since l1_ratio=0"): + clf = plqERM_ElasticNet( + loss={"name": "mse"}, + C=C, + l1_ratio=0.0, + omega=np.ones(n_features), + max_iter=5000, + tol=1e-4, + ) + clf.fit(X_scaled, y)