notebooks/OPT_robust_learning_uncertainty_set_deprecated.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ロバスト最適化における不確実集合の推定\n",
    "\n",
    "参考\n",
    "* [Learning for Robust Optimization](https://arxiv.org/abs/2305.19225)：これが一番いいかも？\n",
    "* [Data-Driven Robust Optimization](https://arxiv.org/abs/1401.0212)：仮説検定を使ってChance-constraint Robust Optimizationを実現する手法？Uncertainty set自体は推定しない．\n",
    "\n",
    "ロバスト最適化では不確実集合全てにおいて制約を満たすような最適化を目指しますが，不確実集合の選択の仕方次第で性能がだいぶ変わります．\n",
    "不確実集合は小さければ小さいほど嬉しいですが，どれくらい小さくできるでしょうか？\n",
    "\n",
    "今回は，集めたデータから不確実集合を推定する手法について見てみます．\n",
    "また，通常のロバスト最適化では単一の問題に対して不確実集合をチューニングしますが，今回は似た問題（パラメータが異なるfamily）に対して不確実集合を考えます．\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "次の問題を考えましょう．\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& x(\\theta, y) \\in \\operatorname{argmin} \\quad f(x, y) \\\\\n",
    "& \\text { subject to } g(x, u, y) \\leq 0 \\quad \\forall u \\in \\mathcal{U}(\\theta) \\text {, } \\\\\n",
    "&\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "ここで，\n",
    "\n",
    "* $y \\in \\mathbf{R}^p$はfamily parameter\n",
    "    * 問題固有のパラメータです．\n",
    "* $x \\in \\mathbf{R}^n$は最適化変数\n",
    "* $u \\in \\mathbf{R}^m$は不確実パラメータ\n",
    "    * これは$\\theta$によってパラメータ化された凸集合$\\mathcal{U}(\\theta) \\subseteq \\mathbf{R}^m$から値を取ります．\n",
    "* $f: \\mathbf{R}^n \\times \\mathbf{R}^p \\to \\mathbf{R}$は目的関数\n",
    "* $g: \\mathbf{R}^n \\times \\mathbf{R}^m \\times \\mathbf{R}^p \\to \\mathbf{R}^d$は不確実制約\n",
    "    * これは$g(x, u, y)=\\max _{l=1, \\ldots, L} g_l(x, u, y)$の形式を取るとして，$g_l$は凹関数とします．  \n",
    "\n",
    "このときに，解$x(\\theta, y)$が様々な$y$のもとで，いい感じになるような$\\theta$を見つけることがゴールです．\n",
    "より正確には，次を解くことを目指します：\n",
    "\n",
    "$$\n",
    "\\mathbf{P}_{(u, y)}(g(x(\\theta, y), u, y) \\leq 0) \\geq 1-\\eta\n",
    "$$\n",
    "\n",
    "\n",
    "ここで，$\\mathbf{P}_{(u, y)}$は$u$と$y$についての同時分布です．\n",
    "つまり，解$x(\\theta, y)$が制約を満たす確率が$1-\\eta$以上です．\n",
    "\n",
    "この制約は次のValue at riskが非正になるのと等価です．\n",
    "\n",
    "$$\n",
    "\\operatorname{VaR}(g(x(\\theta, y), u, y), \\eta)=\\inf \\left\\{\\gamma \\mid \\mathbf{P}_{(u, y)}(g(x(\\theta, y), u, y) \\leq \\gamma) \\geq 1-\\eta\\right\\} \\leq 0\n",
    "$$\n",
    "\n",
    "つまり，解$x(\\theta, y)$が$g(x(\\theta, y), u, y) \\leq \\gamma$を満たす確率が$1-\\eta$以上の制約のうち，最小の$\\gamma$の値です．\n",
    "これが非正なら，結局もとの確率の制約を満たすのと同じです．\n",
    "\n",
    "実はこれは計算することがほぼできません（Chance constraintの評価はNP困難です）．\n",
    "その代わりに，もっと扱いやすいConditional value at riskを考えましょう．\n",
    "\n",
    "$$\n",
    "\\mathbf{C V a R}(g(x(\\theta, y), u, y), \\eta)=\\inf _\\alpha\\left\\{\\mathbf{E}_{(u, y)}\\left((1 / \\eta)(g(x(\\theta, y), u, y)-\\alpha)_{+}\\right)+\\alpha\\right\\}\n",
    "$$\n",
    "\n",
    "ここで，$(a)_{+}=\\max \\{a, 0\\}$です．\n",
    "上で定義したものには次の関係があります：\n",
    "\n",
    "$$\n",
    "\\operatorname{CVaR}(g(x, u, y), \\eta) \\leq 0 \\Longrightarrow \\operatorname{VaR}(g(x, u, y), \\eta) \\leq 0 \\Longleftrightarrow \\mathbf{P}(g(x, u, y) \\leq 0) \\geq 1-\\eta\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 簡単な例（新聞販売の例）\n",
    "\n",
    "最初に簡単な問題設定を通して，全体像を把握しましょう．\n",
    "新聞販売を例に出して，どのようなロバスト最適化を考えるべきか見てみます．\n",
    "\n",
    "* 毎朝，新聞販売員は$n=2$個の商品$x \\in \\mathbf{R}^n$を販売することを考えます．$x$のそれぞれの要素は在庫を表します．\n",
    "* 商品のコストは$k \\in \\mathbf{R}^n$です．\n",
    "* この商品は価格$p\\in \\mathbf{R}^n$で販売され，ここで$p > k$です．\n",
    "* 商品の販売は需要$u$か在庫$x$がなくなるまで行われます．\n",
    "\n",
    "このとき，目的関数はコストから収益を引いた値になります：\n",
    "\n",
    "$$\n",
    "k^T x-p^T \\min \\{x, u\\}\n",
    "$$\n",
    "\n",
    "さて，需要が特定の値ではなく，不確実集合のどれかの値を取る場合を考えましょう．\n",
    "つまり，次の問題を考えます：\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\operatorname{minimize} & \\tau \\\\\n",
    "\\text { subject to } & k^T x+\\max \\left\\{-p^T x,-p^T u\\right\\} \\leq \\tau \\quad \\forall u \\in \\mathcal{U}(\\theta) \\\\\n",
    "& x \\geq 0\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "さらに，需要$u$は次の平均と分散の対数正規分布に従うとします．\n",
    "\n",
    "$$\n",
    "\\mu=\\left[\\begin{array}{l}\n",
    "0.9 \\\\\n",
    "0.7\n",
    "\\end{array}\\right], \\quad \\Sigma=\\left[\\begin{array}{rr}\n",
    "0.6 & -0.4 \\\\\n",
    "-0.3 & 0.1\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "\n",
    "この値は未知ですが，$N=50$個の$u$のサンプルがある場合を考えます．\n",
    "このときに，$\\mathcal{U}$を推定し，適切なロバスト最適化がしたいわけです．\n",
    "\n",
    "---\n",
    "\n",
    "**family parameterについて**\n",
    "\n",
    "$y = (k, p)$としましょう．上の新聞販売の問題設定は，コスト$k$と価格$p$によってパラメータ化されています．\n",
    "\n",
    "このパラメータは有限のサポートを持つとします．特に，$y_j=\\left(k_j, p_j\\right)$ for $j=1, \\ldots, 8$の８パターンがあり得るとしましょう．この$y$は次の手順で生成されるとします：\n",
    "\n",
    "* $k_j$は$[2, 6]$から一様分布に従って生成されます．\n",
    "* $p_j=k_j + r$とし，$r$は$[2, 4]$の一様分布に従って生成されます．\n",
    "\n",
    "---\n",
    "\n",
    "**普通のUncertainty setの構築方法**\n",
    "\n",
    "最初に，通常のUncertainty setの構築方法を考えましょう．\n",
    "\n",
    "$$\n",
    "\\mathcal{U}(\\theta)=\\left\\{\\hat{\\mu}+\\hat{\\Sigma}^{1 / 2} u \\mid\\|u\\|_2 \\leq \\rho\\right\\}=\\left\\{u \\mid\\left\\|A^{\\mathrm{st}} u+b^{\\mathrm{st}}\\right\\|_2 \\leq \\rho\\right\\}\n",
    "$$\n",
    "\n",
    "ここで，$\\rho$はハイパーパラメータです．\n",
    "$\\rho$をいじることで，Uncertainty setのサイズを調節できます．\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.19725539 0.95420074]\n",
      "[[ 1.76857532 -0.25195458]\n",
      " [-0.25195458  1.28694971]]\n",
      "[[0.75953395 0.06897572]\n",
      " [0.06897572 0.891385  ]]\n",
      "[-0.97517281 -0.93314178]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_79745/1429582079.py:31: RuntimeWarning: covariance is not symmetric positive-semidefinite.\n",
      "  u = np.random.multivariate_normal(norm_mu, norm_Sigma, size=N)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff6a621fcd0>]"
      ]
     },
     "execution_count": 268,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGwCAYAAACOzu5xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABF1ElEQVR4nO3deXhU5f338c/MJJmsM0lIQgIkEDYhsssqIC6oUItLW7dqxa2tFvf2V7XPryJtH9H2abVWi2vRSnGplVJthbpCUZAlokSQNRCWLCQhM1nIJJk5zx8hKYHs25mTvF/XNVeYycnMNwPDfObc3/u+bYZhGAIAALAAu9kFAAAAtBbBBQAAWAbBBQAAWAbBBQAAWAbBBQAAWAbBBQAAWAbBBQAAWEaI2QV0RCAQ0JEjRxQTEyObzWZ2OQAAoBUMw1Bpaan69esnu71t51AsHVyOHDmi1NRUs8sAAADtcPDgQQ0YMKBNP2Pp4BITEyOp9hd3uVwmVwMAAFrD6/UqNTW1/n28LSwdXOqGh1wuF8EFAACLaU+bB825AADAMgguAADAMgguAADAMgguAADAMgguAADAMgguAADAMgguAADAMgguAADAMgguAADAMiy9ci7M5Q8Y2phdrILSSiXFhGtyerwcdja7BAB0HYIL2mVVVq4Wvb1duZ7K+ttS3OFaOC9Dc0almFgZAKAnY6gIbbYqK1e3L8tsEFokKc9TqduXZWpVVq5JlQEAejqCC9rEHzC06O3tMhr5Xt1ti97eLn+gsSMAAOgYggvaZGN28WlnWk5mSMr1VGpjdnH3FQUA6DUILmiTgtKmQ0t7jgMAoC0ILmiTpJjwTj0OAIC2ILigTSanxyvFHa6mJj3bVDu7aHJ6fHeWBQDoJQguaBOH3aaF8zIk6bTwUnd94bwM1nMBAHQJggvabM6oFC25foKS3Q2Hg5Ld4Vpy/QTWcQEAdBkWoEO7zBmVogszklk5FwDQrQguaDeH3aZpQ/qYXQYAoBdhqAgAAFgGwQUAAFgGwQUAAFgGwQUAAFgGwQUAAFiGqcFl0KBBstlsp10WLFhgZlkAACBImTodetOmTfL7/fXXs7KydOGFF+rKK680sSoAABCsTA0uiYmJDa4/+uijGjJkiGbNmtXo8T6fTz6fr/661+vt0voAAEBwCZoel6qqKi1btkw333yzbLbGV19dvHix3G53/SU1NbWbqwQAAGayGYZhmF2EJL3xxhv67ne/q5ycHPXr16/RYxo745KamiqPxyOXy9VdpQIAgA7wer1yu93tev8OmiX/X3zxRc2dO7fJ0CJJTqdTTqezG6sCAADBJCiCy4EDB/T+++/rrbfeMrsUAAAQxIKix2Xp0qVKSkrSJZdcYnYpAAAgiJkeXAKBgJYuXar58+crJCQoTgABAIAgZXpwef/995WTk6Obb77Z7FIAAECQM/0Ux0UXXaQgmdgEAACCnOlnXAAAAFqL4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACzD9OBy+PBhXX/99erTp48iIiI0evRobd682eyyAABAEAox88GPHTum6dOn67zzztO7776rxMRE7d69W3FxcWaWBQAAgpSpweWxxx5Tamqqli5dWn9benq6iRUBAIBgZupQ0T/+8Q9NnDhRV155pZKSkjR+/Hg9//zzTR7v8/nk9XobXAAAQO9hanDZt2+flixZomHDhmn16tW6/fbbddddd+nll19u9PjFixfL7XbXX1JTU7u5YgAAYCabYRiGWQ8eFhamiRMn6tNPP62/7a677tKmTZu0fv360473+Xzy+Xz1171er1JTU+XxeORyubqlZgAA0DFer1dut7td79+mnnFJSUlRRkZGg9tGjhypnJycRo93Op1yuVwNLgAAoPcwNbhMnz5dO3fubHDbrl27NHDgQJMqAgAAwczU4HLvvfdqw4YNeuSRR7Rnzx4tX75czz33nBYsWGBmWQAAIEiZGlwmTZqkFStW6NVXX9WoUaP0y1/+Uk888YSuu+46M8sCAABBytTm3I7qSHMPAAAwh2WbcwEAANqC4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACzD1ODy8MMPy2azNbiMGDHCzJIAAEAQCzG7gDPPPFPvv/9+/fWQENNLAgAAQcr0lBASEqLk5ORWHevz+eTz+eqve73erioLAAAEIdN7XHbv3q1+/fpp8ODBuu6665STk9PksYsXL5bb7a6/pKamdmOlAADAbDbDMAyzHvzdd99VWVmZzjjjDOXm5mrRokU6fPiwsrKyFBMTc9rxjZ1xSU1Nlcfjkcvl6s7SAQBAO3m9Xrnd7na9f5saXE5VUlKigQMH6ne/+51uueWWFo/vyC8OAADM0ZH3b9OHik4WGxur4cOHa8+ePWaXAgAAglBQBZeysjLt3btXKSkpZpcCAACCkKnB5Sc/+YnWrFmj/fv369NPP9UVV1whh8Oha6+91syyAABAkDJ1OvShQ4d07bXXqqioSImJiZoxY4Y2bNigxMREM8sCAABBytTg8tprr5n58AAAwGKCqscFAACgOQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGQQXAABgGUETXB599FHZbDbdc889ZpfSa/kDhtbvLdLKrYe1fm+R/AHD7JIAAGggxOwCJGnTpk169tlnNWbMGLNL6bVWZeVq0dvbleuprL8txR2uhfMyNGdUiomVAQDwX6afcSkrK9N1112n559/XnFxcWaX0yutysrV7csyG4QWScrzVOr2ZZlalZVrUmUAADRkenBZsGCBLrnkEs2ePbvFY30+n7xeb4MLOsYfMLTo7e1qbFCo7rZFb29n2AgAEBRMDS6vvfaaMjMztXjx4lYdv3jxYrnd7vpLampqF1fY823MLj7tTMvJDEm5nkptzC7uvqIAAGiCacHl4MGDuvvuu/WXv/xF4eHhrfqZBx98UB6Pp/5y8ODBLq6y5ysobTq0tOc4AAC6kmnNuVu2bFFBQYEmTJhQf5vf79fatWv11FNPyefzyeFwNPgZp9Mpp9PZ3aX2aEkxrQuNrT0OAICuZFpwueCCC7Rt27YGt910000aMWKE7r///tNCC7rG5PR4pbjDleepbLTPxSYp2R2uyenx3V0aAACnMS24xMTEaNSoUQ1ui4qKUp8+fU67HV3HYbdp4bwM3b4sUzapQXixnfi6cF6GHHZbIz8NAED3Mn1WEcw3Z1SKllw/QcnuhsNBye5wLbl+Auu4AACChs0wDMvOc/V6vXK73fJ4PHK5XGaXY3n+gKGN2cUqKK1UUkzt8BBnWgAAna0j799BsXIugoPDbtO0IX3MLgMAgCYxVAQAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyD4AIAACyjXcElEAg0eXtOTk6HCgIAAGhKm4KL1+vVVVddpaioKPXt21cPPfSQ/H5//fePHj2q9PT0Ti8SAABAauMmiz//+c/1xRdf6JVXXlFJSYl+9atfKTMzU2+99ZbCwsIkSRbebBoAAAS5Np1x+fvf/65nn31W3/nOd3Trrbdq8+bNOnr0qObNmyefzydJstlsXVIoAABAm4LL0aNHNXDgwPrrCQkJev/991VaWqpvfOMbqqio6PQCAQAA6rQpuKSlpWnHjh0NbouJidG///1vHT9+XFdccUWnFgcAAHCyNgWXiy66SEuXLj3t9ujoaK1evVrh4eGdVhgAAMCp2tScu2jRIh05cqTR78XExOi9995TZmZmpxQGAABwqjYFl7i4OMXFxekXv/hFs8fNmjWrQ0UBAAA0pk3Bpc6KFSsaXK+urlZ2drZCQkI0ZMgQPfTQQ51SHAAAwMnaFVw+//zz027zer268cYbadAFAABdptP2KnK5XFq0aJF+/vOfd9ZdAgAANNCpmyx6PB55PJ7OvEsAAIB67RoqevLJJxtcNwxDubm5euWVVzR37txOKQwAAOBU7Qoujz/+eIPrdrtdiYmJmj9/vh588MFOKQwAAOBU7Qou2dnZnV0HAABAizq1xwUAAKArEVwAAIBlEFwAAIBlEFwAAIBlEFwAAIBlEFwAAIBlmBpclixZojFjxsjlcsnlcmnatGl69913zSwJAAAEMVODy4ABA/Too49qy5Yt2rx5s84//3xddtll+uqrr8wsCwAABCmbYRiG2UWcLD4+Xr/5zW90yy23tHis1+uV2+2Wx+ORy+XqhuoAAEBHdeT9u10r53YFv9+vv/71ryovL9e0adMaPcbn88nn89Vf93q93VUeAAAIAqY3527btk3R0dFyOp267bbbtGLFCmVkZDR67OLFi+V2u+svqamp3VwtAAAwk+lDRVVVVcrJyZHH49Gbb76pF154QWvWrGk0vDR2xiU1NZWhIgAALKQjQ0WmB5dTzZ49W0OGDNGzzz7b4rH0uKAn8AcMbcwuVkFppZJiwjU5PV4Ou83ssgCgy/SIHpc6gUCgwVkVoCdblZWrRW9vV66nsv62FHe4Fs7L0JxRKSZWBgDBydTg8uCDD2ru3LlKS0tTaWmpli9fro8//lirV682syygW6zKytXtyzJ16inPPE+lbl+WqSXXTyC8AMApTA0uBQUFuuGGG5Sbmyu3260xY8Zo9erVuvDCC80sC+hy/oChRW9vPy20SJIhySZp0dvbdWFGMsNGAHASU4PLiy++aObDA6bZmF3cYHjoVIakXE+lNmYXa9qQPt1XGAAEOdOnQwO9UUFp06GlPccBQG9BcAFMkBQT3qnHAUBvQXABTDA5PV4p7nA11b1iU+3sosnp8d1ZFgAEPYILYAKH3aaF82oXWTw1vNRdXzgvg8ZcADgFwQUwyZxRKVpy/QQluxsOByW7w5kKDQBNCLoF6IDeZM6oFF2YkczKuQDQSgQXwGQOu40pzwDQSgwVAQAAyyC4AAAAyyC4AAAAyyC4AAAAyyC4AAAAyyC4AAAAyyC4AAAAyyC4AAAAyyC4AAAAyyC4AAAAy2DJf1iOP2Cwtw8A9FIEF1jKqqxcLXp7u3I9lfW3pbjDtXBeBrspA0AvwFARLGNVVq5uX5bZILRIUp6nUrcvy9SqrFyTKgMAdBeCCyzBHzC06O3tMhr5Xt1ti97eLn+gsSMAAD0FwQWWsDG7+LQzLSczJOV6KrUxu7j7igIAdDuCCyyhoLTp0NKe4wAA1kRzLtrMMAz5agIq89WorLKm9quvRuW+GlXVBFTlD6jGb6jaHzhxqf2z3zBkt9nksNlks6n2z3ab7DbJZrMpzGFXeJhDEaEOhYfaT3ytvdhtrZs1lBQT3sW/PQDATAQXSKoNI0fLfDp87LjyPJUqLPOpsKxKReU+FZVVqaisSoVlPhVXVKmsskY1QdhLEuV0aEeuV0XlPiVGO5UYU3uJdobI1srgAwAIbgSXXqTaH9CBonLtzi/TvsJyHTp2XIeOVejwseM6XHJcvppAm+8zKsyh6PAQRTtDFOUMkTPErlCHXSEOu8IcNoXY7QoNsSvUbpPdbpNhSAHDOHGp/bNhGPIHDFX7DVVW+3W82q/jVX5VVvtVWR2ovV7tV1UL9ZX7/PrFO9tPuz081K7EGKeSXeHqFxuh/rER6h9X+3VAXIT6xUYoMoyXAgBYAf9b90CGYSinuEJfHvJoV36p9hSUaXdBmfYXljd7psRmk5Jd4Upxhysxxqk+0U4lRIXVfo12qk90mOKjwhRTF1TCQmTvxoXffDV+rfz8sH6zeqeOllXV3+4KD9Gs4YmKjwpTYVmVjpb6dLTMp6OlPpX5alRZHdDB4uM6WHxc0rFG7zsuMlSp8ZFKT4jSoD5RGpxY+zU9MUqu8NBu+g0BAC2xGYYRfOf8W8nr9crtdsvj8cjlcpldjmnyPJXaerBE2w6X6MtDHn15yCPP8epGj40Kc2hoUrSGJEYrNT5S/eNqzzoMiI1UsjtcYSHB36/dlpVzK6pqVFhapaNllcr1VNafXTpSclyHTvy5tLKm2cdLiA6rDzPD+8ZoRLJLw5OjlRjtZAgKANqhI+/fBBcLyvdWav3eotrLviLlFFecdkyYw66R/VzKSHFpaFK0hiVFa2hStFLc4bzZnsJbWa3Dx44rp7hC2YXl2l9Yrn2F5couLNfRUl+TPxcfFabhfaNrg0zfGI1IidHIZJciwhzdWD0AWA/BpYcHl8pqvz7dW6gPdhRo/d4i7Sssb/B9h92m4X1jNKa/W2NS3Ro7IFbD+8ZY4uxJsCutrNaBogrtKyzXnoIy7cor1c78Uu0vKldjrxy7TRqWFKPRA9wa3d+tUf3dykghzADAyQguPTC4FJdX6cOvC/T+9nyt3X1UFVX++u/ZbdKo/m5NHdxH0wb30aT0eEU7aVfqTser/NpTUKad+aXalV+qr/NKtf2IV4Vlp5+hcdhtGpYUrdH93ZowME5nDYzT0MTobu0PAoBgQnDpIcGl3FejVVl5euvzQ1q/t0gn99GmuMM1e2RfzRqeqEnp8XJH0DAajPK9lfrykEfbDnuUdbi236ixMOMKD6kNMWlxOmtQnMalxjKzCUCvQXCxcHAJBAxtyC7S37Yc1rtZuQ3OrIxMcenCjL66KKOvzuznojfFggzDUL7Xp22HPfriYIm2HDimrQdLdLza3+A4h92mjBSXpg6O19lDEnrFWbS2NFkD6FksG1wWL16st956S19//bUiIiJ09tln67HHHtMZZ5zRqp+3cnAp89XojU0HtfTT7BPTdGsN6hOpb00YoCvG91dqfKSJFaKr1PgD2pFbqi0HirX5wDFtOXDstH2YHHabxgxw6+whfTRtcILOGhjXo/pkVmXlatHb2xv83inucC2cl6E5o1JMrAxAd7BscJkzZ46uueYaTZo0STU1NfrZz36mrKwsbd++XVFRUS3+vBWDy6FjFXr50/16beNBlfpqp+HGhIfom2P66Ttn9deEtDjOrPRCR0qOa9P+Yq3fW6RP954+UyzMYdf4tFjNOiNR5w5P0siUGMv+O1mVlavbl2WettN33W+z5PoJhBegh7NscDnV0aNHlZSUpDVr1uicc8457fs+n08+33/7Bbxer1JTUy0RXHKKKvTE+7u08osj8p9oXhmSGKWbZ6TrW+MH9KhP0+i4Q8cq6qe8f7q3SHnehmdk+rqcmjU8UeeekaTpQxMs0/PkDxia8diHTe70bZOU7A7XuvvPZ9gI6ME6ElyCahDd4/FIkuLj4xv9/uLFi7Vo0aLuLKnDjpb69NSHu7V8Y46q/bWBZcbQBN0yI12zhicyswSNGhAXqSsnRurKiakyDEP7iyq0bvdRfbzzqD7dW6R8r09vbD6kNzYfksNu01lpcTpvRJIuOrOvhiRGm11+kzZmFzcZWiTJkJTrqdTG7GJNG9Kn+woDYBlBc8YlEAjo0ksvVUlJidatW9foMVY641JZ7dcfP96r59fuq2/EnDksQf9z8RkaMyDW3OJgaZXVfm3aX6yPdx7VxzsLtPdow3V9hiZF6+Iz++riM5M1ur87qIaUVm49rLtf29ricb+/ZpwuG9e/6wsCYIoeccZlwYIFysrKajK0SJLT6ZTT6ezGqtrn072F+tlb27S/qLZPYWxqrO6/+AydPTTB5MrQE4SHOjRzWKJmDkvUz7+ZoYPFFfp4Z4He21Gg9XsLtaegTHsKyvT0R3uV4g7XRRm1IWZyerxCHOYuSpgUE96pxwHofYLijMsdd9yhlStXau3atUpPT2/1zwVbc25JRZUe+dcOvbH5kKTaPoSF887U3FHJQfWptyswtTU4eI5X6+OdBfr3V/n6aGdBg+n1CdFhumR0ii4d18+0JvC6Hpc8T+VpzbkSPS5Ab2HZ5lzDMHTnnXdqxYoV+vjjjzVs2LA2/XwwBZfN+4u1YHmm8r21Q1nXT03TT+eM6BU7CzO1NThVVvv1yZ5Crf4qT+9tz9exiv9uvNk/NkLzxvbTvLEpykjp3jWC6mYVSWoQXphVBPQelg0uP/rRj7R8+XKtXLmywdotbrdbERERLf58MAQXwzC09JP9euRfO1QTMDQ4MUq//vYYTRzUeINxT8PUVmuo9ge0bk+h3t56RKu/ylP5SWdihiRG6Yrx/fWtCQPUL7bl111nIOwCvZtlg0tTn/KWLl2qG2+8scWfNzu4HK/y63/e/ELvfJkrSfrmmBQ99u0xiurhK57WYWqrNVVW+/Xh1wX6x9Yj+nBngapqApIkm02aOSxRV541QBdm9FV4aNdO0Wd4Eei9LBtcOsrM4FJaWa1bXtqsjfuLFWK36X8vGan5Zw/q8b0sJ1u/t0jXPr+hxeNe/f5UprYGKW9ltVZn5elvmYe0YV9x/e3uiFBdPq6frpyYynYTADpdj5hVZCXHyqt049KN+uKQRzHhIXpx/iRNTu8dQ0MnKyhtej2O9hyH7ucKD9WVE1N15cRUHSgq15tbDunNLYeU66nUy+sP6OX1B5SR4tL3pg3UZeP6sREkANOZOzfSgkoqqnTt8xv0xSGP4iJD9er3p/bK0CIxtbWnGdgnSj++6Aytu/98vXzzZH1zTIrCHHZtz/Xqwbe2acojH+gXb29XdmF5y3cGAF2Ej09tUO0P6Ed/ydTXeaVKinHqL7dO0bC+MWaXZZrJ6fFKcYe3OLW1twY7q3LYbZo1PFGzhifqWHmV/rrloJZtyFFOcYX+9Em2/vRJtmYOS9AN0wbp/BFJ9KUA6Fb0uLSSYRj6379n6S+f5SgqzKG3fjRdZyT33tBSh6mtvUMgYGjNrqN6ZcMBfbSzQHX/a6TGR+jWGYN15cQBDCMBaDWac7shuPzlswP6PyuyZLNJL9wwUReM7Nulj2clTG3tXXKKKvSXzw7o9c0HVXJibZjYyFDdMHWgbjh7kBKig391awDmIrh0cXA5dKxCFz2+VhVVfj0wd4RumzWkyx7Lqpja2vscr/Lrr1sO6oX/ZCunuHZ7C2eIXd85a4BunTlY6QlRJlcIIFgRXLowuBiGoVte3qwPvy7Q5EHxeu0HU9nRGTiJP2BoVVaenlu7V18cqt3h3WaT5o3pp7suGKqhSQypAmiI4NKFweWfX+ZqwfJMhTns+tfdM/hPGGiCYRj6LLtYz63dpw+/LpBUG2AuHdtPd54/TEOTok2uEECwYB2XLmIYhn7/wS5J0m3nDiG0AM2w2WyaOriPpg7uo6zDHj35wW79e3u+Vm49on98cYQAA6BTsI5LM9bsOqpd+WWKCnPolhmt37Ua6O1G9XfruRsm6p07Z+iijL4yDGnl1iO68PE1uu+NrTpSctzsEgFYFMGlGS/8J1uSdPWkNLkjev4uz0BnayzAvJV5WOf9v4/16Ltfy1tZ3fKdAMBJCC5NyPNUat2eQtls0k3TB5ldDmBpdQFm5YLpmpIeL19NQM+s2atZv/5If1qXXb/RIwC0hODShLW7jkqSxgyIVWp8pMnVAD3D2NRYvfaDqXrhhokamhStYxXV+sU72zX7d2u0+qs8WXiuAIBuQnBpwprdtcFl1vBEkysBehabzabZGX216u6ZeuSK0UqIdiqnuEI/fGWLbn5pkw4UsRcSgKYRXJqwfm+RJGnW8ASTKwF6phCHXd+dkqY1/3OufnTuEIU6bPpo51Fd+PhaPf7eLlVW+80uEUAQIrg0orSyWsXlVZKkEcldu5UA0NtFOUP00zkjtOqeczRjaIKqagL6/Qe7ddHja/XRifVgAKAOwaURR0pq99yJjQxVlJOlboDuMCQxWq/cMllPf3eCkl3hyimu0E0vbdK9r29VSUWV2eUBCBIEl0bUrTHRzx1hciVA72Kz2XTJmBR98ONZ+v7MdNlt0orPD+vCx9fq/e35ZpcHIAgQXBpRUVU7th7N2RbAFFHOEP2fSzL05u1na3BilI6W+nTrnzdz9gUAwaUx4aG1T4uvxrrNgf6AofV7i7Ry62Gt31skf4BpprCeCWlx+tddM/XDcwY3OPvyyZ5Cs0sDYBJOKTQiPNQhSaqstuaiWKuycrXo7e3K9VTW35biDtfCeRmaMyrFxMoa8gcMbcwuVkFppZJiwjU5PV4Odt7GKcJDHXrwGyN18ahk/c9fv9Deo+W6/sXPtODcobpn9jCFOPj8BfQmBJdGRIbVBhfPcestR74qK1e3L8vUqedX8jyVun1ZppZcPyEowotVwhWCx4S0OL1z50z94p2v9OrGg3rqoz36LLtIT147Xin0owG9Bh9VGjE4sXb32jxvpaXCiz9gaNHb208LLZLqb1v09nbTh43qwtXJoUX6b7halZVrUmUIdhFhDi3+1hg9ee14RTtDtGn/Mc39/X/0wQ4ad4HeguDSCHdEqPrH1n6C25lXanI1rbcxu/i0MHAyQ1Kup1Ibs4u7r6hTtCdc0a+DU106tp/euXOGRvd3q6SiWrf+ebOe/mgPWwYAvQBDRU0YkRyjwyXHtf2IR5PT480up1UKSpsOLe05riu0JVxNG9KHISU0aVBClN68fZp++c52LduQo9+s3qld+aV67Ntj6vvUAPQ8nHFpwoSBcZKkNSc2W7SCpJjwTj2uK7QlXDGkhJY4Qxz61eWj9cvLR8lht2nl1iO6+rkNKvCaF84BdC2CSxMuPjNZkrRuT6G8ldboc5mcHq8Ud7iampdjU+3ZCjPPILU2NCVEOzutX4ehpp7ve1MH6pVbJis2MlRfHCzRpU99YqlhXgCtR3BpwtCkaA1Nila137DMfikOu00L52VI0mnhpe76wnkZnTLluLEw0JqA0NpwJUOd0q+zKitXMx77UNc+v0F3v7ZV1z6/QTMe+5CzNT3Q2UMStHLBdA1Nilaet1JXP7deWw+WmF0WgE5mMyzczeb1euV2u+XxeORydf5miL/990794cM9mjksQa/cMqXT77+rdHVfSGP3HxsZKkkqqfjv2ammHrNuCEhSgzMqdWFmyfUT5KsJ6O7XtrZYy++vGafLxvVvss7Gpoaf/Dj0yfQ8JRVVuumlTfo8p0RRYQ49f8NEnT2UXd6BYNKR92+CSzMOFldo1m8+UsCQ/nXXTGX0s85O0S0t7tbexd+aCgONaS4gtBSu1u8t0rXPb2jxMV79/lRNG9LntNv9AUMzHvuwybM2NknJ7nCtu/98Fr3rgcp9NfrBK5v1yZ4ihYXY9dS143XRieFfAOYjuHRRcJGkO1/9XG9/cUSXj+unJ64Z3yWP0d3ae0ampTDQmOYCQnPhqe6x8jyVjYakloJHR4NPR7AicHDw1fh116ufa/VX+Qqx2/TcDWfp/BF9zS4LgDr2/k2PSwt+eM5gSdLbX+bqYHGFydV0XEdm6rQ0lbkxzfWiOOw2TRvSR5eN669pQ/o0eHPvaL+OWVPD6akJHs4Qh57+7gRdOrafagKGbl+Wqc/2FZldFoAOMjW4rF27VvPmzVO/fv1ks9n097//3cxyGjWqv1szhyXIHzD0f/+5w+xyOqSjK+u+tz2v3Y/dnoAwZ1SKllw/QcnuhjORkt3hLfanmDE1nOnbwSfEYddvrxqrC0YkyVcT0K0vb1bWYY/ZZQHoAFODS3l5ucaOHaunn37azDJa9L+X1H6yX/VVntZaaF2XU3VkZd1VWbn60yf72/3Y7Q0Ic0alaN395+vV70/V768Zp1e/P1Xr7j+/xaba7p4abpXtFnqjUIddT183QVPS41Xqq9ENf9qo7MJys8sC0E6mBpe5c+fqV7/6la644gozy2jRGckxuvHsQZKkh//xlXw1fnMLaqe2Dp/UTW9ekXlIP1uR1a7H7IyA0NyQUnM/011TwyVrbLfQm4WHOvTC/Ika3d+t4vIq/eDPm1XmqzG7LADtYKkeF5/PJ6/X2+DSXe6ePUwJ0U7tKyzXUx/u6bbH7UxtGT45uVfj3je+UHF5VZsfrysCQlt0ZKiprayw3UJXC/aF/mLCQ/XijRPV1+XU7oIy/fiNrQoEWY0AWmapvYoWL16sRYsWmfLYrvBQPXxphu5Y/rme+miPpg3po7OHWGttiLrhk5Zm6hwr92nB8s9bNeX5ZJFhDlVU/fdsVHIQ7Ck0Z1SKLsxI7vJZPlbYbqErWWVPqaSYcC25/ixd8+wGrf4qX0vW7NWC84aaXRaANgia6dA2m00rVqzQ5Zdf3uQxPp9PPp+v/rrX61VqamqXToc+1U/f/EJvbD6kvi6n3r37HMVHhXXL43aWlhZ/e/q7E/TLf25v8+whSfrLrVNkt9naFBB6ytThjk7ftjIrLvT32sYcPfDWNtlstf9urfYhBLC6XjMd2ul0yuVyNbh0t4cvPVNDEqOU7/XpJ3/9wnKnmlsaPomLCmtzaKnrY5k6uE+belF60tTh7u6pCRZWbUq+ZnKarpmUKsOQfvrmlyqn3wWwDEsFl2AQGRaip747QWEhdn34dYEeW/212SW1WXMzddo65bk1b8qN9T70xKnD3dlTEyys3JT8v9/MUP/YCB06dlyPrbLe6xjorUztcSkrK9OePf9tdM3OztbWrVsVHx+vtLQ0Eytr3sgUlx779mjd+/oXenbNPqXFR+q6KQPNLqtN6mbqnKw9U56T3eH6+SUZckeEaeXWw6cN9zTW+5DscqqyJtDkp3Sbaj+lX5iRbLkzFB3tqbHa0JmVm5KjnSH69XfG6LoXPtOf1x/Q3FEpnb6KMoDOZ2pw2bx5s84777z66/fdd58kaf78+XrppZdMqqp1rhg/QDlFx/X4+7v00Mqv1C82QuedkWR2We1Wd8q/NeKjQvXzb56pZFdtI++pPTF1TZmSGu19yPP61JyTP6Vb8Y2ksVDYGlZpcD2Z1ZuSpw9N0LWT0/Tqxhw9tDJLq+45J6iDIgCTh4rOPfdcGYZx2iXYQ0uduy4Yqm9PGCB/wNAdf8lUZs4xs0tqt7Ys5//IFaN1xfj+8hyv0oLlnzc53PPAW9vaPDPpZMH4Kb2rWHXorLsX+usKD8wdIXdEqHYXlGnF54fNLgdAC+hx6QCbzabF3xqt6UP7qLzKrxte3GjZ8NLakHDL9EGaMyqlxaZMQ1JJRXWHagrWT+mn6uj6JV3R4Npda6r0hKZkd0Sobj93iCTp8fd2WXaBSaC3sNQ6LsEoLMSu52+YqJuWbtJn2cW64cWN+vMtkzUhLc7s0tqktSFhdkay/AFDL32S3a4p061RN3U4mD+l1+mM4Z22NLi2Zgiqu4ec6pqST+tlCvJhrpPNnzZISz/J1uGS43pt40HNP7FSNoDgwxmXThAZFqKlN03SlPR4lflqdMOLG7XlgLXOvLT2lP+xcp9mPPahftlFG05a5VO61HnDO53Z4GrWkFN795QKFhFhDt1xYiG6lz/dryBZ3gpAIwguneTU8HL9C5/pgx35ZpfVag67TZeOTWm2J+XSsSmN9rS0h01SbGSokl3WnDrcmcM7ndXg2pGaOmNoqT17SgWTb00YoKgwh/YVlmvTfmt98ID5gn3Li56EoaJOVBdebluWqbW7jur7f96sX1w2StdPDf6p0quycvXc2uwmv3/rzHT944vcVjfb1gWTYxXVsqnxVXof/dboblmOvyt05vBOa7diaGnorL01WXE2U1eIcoZo7ugUvbnlkK5+br2yF19idkmwCF5D3YszLp0sMixEL86fqKsmDlDAkP7371l69N2vg2KF3aY+ETT3Sb3O3zIPtfpMS13sWPyt0XqmhQXZrPopvTOHdzqrwbU9NVl1NlNXmTmsdul/w5Dyvb1nVhvaj9dQ9+OMSxcIddj12LfHqH9spB5/f5eeWbNXB49V6DffGaPIMHOe8uY+EbgjWl7mv7i89TOETm3KtOpZleZ09volndHg2taaWhpasvJCgO116dh+uvu1rZKkNTuP6qpJqeYWhKDGa8gcBJcuYrPZdPfsYeofF6EH/val/vllrvbkl+mZ752l9ISobq2lqU3w6j4R3DR9UKc91s8vGakbp6c3eJG2d0G2YNZZwzsn6+iqu22tqbNnM/UENptNl43rp5Vbj+hAcbnZ5SDI8RoyB0NFXew7Zw3Q8u9PVWKMUzvzS3XpH9bp31+1bT+gjmhNw+bKrUdadV/xUWFNzjqSantaTg0t7RXsjW5dtX6Jw27T5PR4JcWEq6C09j+81v7uba3Jysv1d6XpJ3aK3rAv+PZXQnDhNWQOgks3mJwer3/eOUOTBsWp1FejH7yyRb9e9bVq/IEuf+zWfCIoKq9SfFRos6Ek2eXUry4b1WwfTElFdZs3aWyMVXaN7opNFTv6u7elJqsv199VzhpUuwbTV0c8TItGs3gNmYOhom6S5ArX8u9P1eJ/fa0/fZKtP368V59lF+vxq8YprU9klz1ua5P+FeP660+f7D9tBlCdypqAJEOxkaFNrojbGeO5LQ1rBdtU6c7cVHF/Ybkef3/3ace09XdvbU1dMdzVE/SPjZAkVVYH5DlerdjIMJMrQrDiNWQOzrh0o1CHXQ/Ny9CT145XjDNEWw4c09zfr9Ubmw522Se71ib9frERunn6IEU6HY1+31NRrR8t/7zZZfxPHs9tj65Y+r47tHdm1KlnVxoLLVL7fvfW1NQTluvvCuGhDsVH1YaVrlodGj0DryFzEFxMcOnYfnr3npmanB6v8iq/fvq3L3Xbsi0qLq/q9MdqaUVcSbLbpF/+c4de/GS/yn2N79PSlqjQ3vHctjS6WV1TUyib0lW/e1cMd/UEfU8sjJjHlGi0gNdQ92OoyCQD4iL16ven6rm1+/S793Zq9Vf5ysxZq19edman/kOv+0Rw+7LMJoeBOvsERnvHc3tLo1tr1s1pSlf87h0d7uqJ6n5zu633PgdoPV5D3YvgYiKH3abbzx2imcMSdM/rW7WnoEy3LcvURRl9teiyM5XijuiUx2lqjRC7rXNDS0fHc3tLo1tLZ5aa09jvfnKfTHv/w+yJU9Y74nh17ZnHiNDGh06BU/Ea6j4ElyAwqr9b79w5Q09/tEdLPt6rf2/P16d7i/TTOWfouikDOyW1n/qJoLDU16GNEptaxr8j47m9pdGtPWdNmvrdO2Op8c4IPj3N8ara4BIZRnABgg09LkEiPNShH190hv5510xNSItVma9GD638St955lNtP+LtlMc4uWEzIcbZ5p+v2yH6j9/tmvHc3tLo1t4zRqf+7p2x1HhHpl8H+1o77eWr8auwzCdJiotiRhEQbGyGhRcq8Hq9crvd8ng8crlcZpfTaQIBQ3/ZmKPH3v1aZb4a2W3S1ZPS9JOLhqtPdNsDR2PW7y3Stc9vaPXxdW+XdeGkKz+l9/QNy/wBQzMe+7DVw0XxUaF65IrRDX73lu6j7gzNuvvPb/Lvpamp56f+XTf1sz317+jLQyW69KlPFBsZqs9/fqFs9LkAna4j798ElyCW56nUr/65Xe98WfvpN8YZorsuGKb5Zw9SWEjHTpbVvfE1NSxzqu5+U+rpwxersnJ127LMVh37+NXjdMX4/g1ua23wfPX7Uxsdd+9I8OlI4LGC5Z/l6GcrtmnmsAS9cssUs8sBeqSOvH8zVBTEkt3heuq7E/TGD6dpVH+XSn01+r//2qGLn1ir97fnd2jtl+aGZRrz80tGduubkVV3jW6tOaNSdO/sYa06Ntl1+tBSR2dgtXfquVXX2mmL9fuKJNX2ngEIPgQXC5icHq+VC2bo198eo4ToMGUXluvWP2/W1c9u0Kb97V/Xo262Ud9G3hhPZlPtOi9WfjMKRnecP6zRUFKnrqeosWbkjs7Aam/w6elr7ZT5auq3rbj4zGSTqwHQGIKLRTjsNl01KVUf/eRc/XDWYIWF2LVxf7GufGa9bly6UVmHPe263zmjUvTbK8c2e8ypb0Z1TZkrPj+sF/+zTysyD/Wo5szu4rDb9PClGbKp7c3ILS0s2FzokdoffKy+1k5LDcWrsvJUWR3Q4IQojR3AGRcgGDEd2mJiwkP14NyRuvHsQXrygz16Y/NBfbzzqD7eeVSXjE7RvRcO19Ck6DbdZ2G5r1XHFZRWNtqUWSfYmjODuU+mrjZfTUD3zB6uVzfmNFilNbmF57K5hQVbMwOrvVPPrbzWTksNxYZhaNmGA5KkK8b3pykXCFI051pc7cZ8u/SPL47IMGoXlZs3tp9+dO5QnZEc06r7aG2j572zh+mJ93c328xrU3A0ZwbzrJfGakt2OXXt5DQNSohqU8jqyO9Z12QrNR58Gvt7bKmpuzWzmczQmoZiyabblm1RRKhDa/7nXCW1MIQKoP2YVdSLg0udHble/fbfu/T+jvz62y7M6Ks7zhuqsamxzf5sa9+MDMNQnrflszMpJr9xBfOsl66orSNnltoTfNoTeMzUmhlUfV1ORYaFaF9hue48f6h+fNEZ3Vsk0MsQXAgu9bIOe/THj/fo3aw81f3NzhiaoAXnDdXUwfFNnv5u6c3ontnD9fj7u1pdR3PTcLty+KYz1jfpKsFaW3v+ToL5jNap2rJmUVxkqNb+9DzFhId2cVVA79aR9296XHqYUf3d+uN1Z2lPQZmWfLxXf996WOv2FGrdnkKNTY3VzdMH6RujUxTqaNiX3dR+RnW9Fr6aQJvqaKw5szve7Noy66W79xUJ1tras8eKlTaVa0uj8P9cPILQAgQ5gksPNTQpWr+9aqzumT1Mz/9nn17bdFBfHCzR3a9t1SP/2qEbpg3SdyenNVjSvLk3o/V7i9r0+Kc2ZzY1RFK3PH1nDS8E86yXYK6tPayyqVxrG4UnDYrTtZNTu7gawDqCdYIDwaWHS42P1C8uG6W7Lhim5Z/l6M/rDyjf69NvVu/Ukx/s1rcmDNBN0wdpeN/aRt6m3oxamoVyslOn4ba0aJlNtYuWXZiR3O4XRd0LrDV77EjmzHqx8owcK2vNv127TXrm+rOYSQScEMzDwazj0kskRDt11wXD9MkD5+l3V43Vmf1c8tUE9OrGHF30+Fpd/ex6rdx6WJXV/kZ//uSVdptj0+nTcLt60bKTNwp8Nyu/2WNbWt+kK3V07RUrM3NDxtasEv3jTtwHDLC6ztjAtStxxqWXcYY49K0JA3TF+P7atP+Yln6SrdVf5emz7GJ9ll2suMhQfXvCAF0zOe209WCa6oOp01Qa78ohkqaGoJpiyLwdpju69opVBcMnt+b+7S44b4gWnNe67ReAnq47zpB3FLOKoFzPcb2+6aBe33SwwX/qU9Lj9d0pabr4zGSFhzrqb68blsnzVqq4zKf4qDAluyOaHP/s6IaATWnrLstS7Vo0d88e3urju0IwvJF3l2Cbml7hq9Gtf96sT/cWySbpiWvG6bJx/Vv8OaC36Kr/r0/FrCJ0SIo7QvfMHq47zx+mNbsKtPyzg/rw6/z6szAxzhDNHZ2sK8YP0JQT4aQt/2Dbu0prS1oagmrMoISoNh3fFaw0I6cjgu2TW1GZTz98ZYs2HzimELtNj317DKEFOIUVJhEERY/L008/rUGDBik8PFxTpkzRxo0bzS6pV3LYbTp/RF+9MH+iPnngfN07e7j6x0ao1FejNzYf0rXPb9DMX3+k36z+WnsKytp0v031GHRkiKQ9L5xgaXzt6btfS8G1IeP2I15d8cdPtfnAMcWEh+jlmyfr22cN6PLHBazGCpMITA8ur7/+uu677z4tXLhQmZmZGjt2rC6++GIVFBSYXVqvluKO0N2zh+k/Pz1Pr/9gqq6ZlKqY8BAdLjmupz/aq9m/W6NLn1qnF/6zT4dLjrd4f3U9Bsnuhv/Yk93h7R4uaMsLpyc3vgarYPjk5g8YembNXl329DrlFFcoNT5CK350tqYPTeiyxwSszAqTCEzvcZkyZYomTZqkp556SpIUCASUmpqqO++8Uw888ECDY30+n3y+/y457/V6lZqaSo9LN6ms9uuDHQV6K/OQ1uw6qpqTZoaMS43VJaNTNGdUslLjI5u8j85cF6ClrQrqBOtS9D1dd42VN+XQsQr9+I0v9NmJMzoXZfTVo98eo/iT1i4CcLru2NbDskv+V1VVKTIyUm+++aYuv/zy+tvnz5+vkpISrVy5ssHxDz/8sBYtWnTa/RBcul9RmU/vfJmrf27L1ab9xTr5X9HYAW59Y3SKvjE6pdkQ0xmaeoGdrKc2vgY7szZkrKz260+fZOupD/eoosqvqDCHFs47U1dOHMA6LUArdfUkAssGlyNHjqh///769NNPNW3atPrbf/rTn2rNmjX67LPPGhzPGZfgVFBaqdVZefrntlxtzC7WyUt0DO8brQtG9tXskUkalxrXJb0cjb3A4qNCdcW4/pqdkdwjG1+tojs3ZDQMQ+9tz9ev/rlDOcUVkqSJA+P026vGamAf85uyAavpypVze82sIqfTKaeTRaKCTVJMuL43bZC+N22Qjpb6tPqrPP1rW64+yy7Wrvwy7cqv3TcpPipM556RqNkj+2rmsIRO2xOmt8zSsaKW9sDqzD2qnnh/lz49sTVFX5dTD8wdocvG9pedfwdAuwTrth6mBpeEhAQ5HA7l5zdc7TQ/P1/JyckmVYWOSIxx6vqpA3X91IEqqajSml1H9cGOAn28s0DF5VV6K/Ow3so8rFCHTZMGxWvGsASdMyxRGSmuDr3BBOsLDF0XLA3D0Lo9hfrDh3vqZyaFhdj1/Znp+tG5QxXltNTnMgCtFBTNuZMnT9Yf/vAHSbXNuWlpabrjjjtOa849FQvQWUe1P6DN+4/pgx35+uDrAmUXljf4fnxUmM4e0kfnDEvUjGEJ6hcbYVKlCHZVNQGt/ipPL6zL1hcHSyRJYQ67rpw4QLfNGtLlfVUAOs6yPS5S7XTo+fPn69lnn9XkyZP1xBNP6I033tDXX3+tvn37NvuzBBfr2ne0TP/ZXaj/7C7U+r2FKq9quEfS4MQoTRvcR5PT4zUlvc9p06jR++QUVWj5xhy9ueWgCsuqJEnhoXZdOzlNPzxnCP9GAAuxdHCRpKeeekq/+c1vlJeXp3HjxunJJ5/UlClTWvw5gkvPUO0PaOvBEv1nd6HW7T6qrQdLdOoefGnxkZqcHn8iyMQrLT6SGSK9QJmvRh/syNffMg9r7a6j9bcnxTh1zaRUfW/aICXG0PcGWI3lg0t7EVx6Js/xam3YV6SN2cXamF2sr454TgsyfV1OjU+N07i0WI1LjdXo/m56GnqIiqoafbCjQP/8Mlcf7SyQryZQ/72ZwxJ03ZSBumBkkkIdpq+fCaCdCC4Elx6ttLJaWw4cqw8yXxwqUbW/4T9bu00a3jdG408EmXGpcRqSGKUQ3tws4dCxCq3dVag1uwq0dlehjlf/d+gwPSFKl4xO0ZUTBzCtGeghCC4El16lstqvrQdLai85tV/zvKcvG+8MsWtEcowy+rmUkeJSRj+XRiS7ODMTBMp9Ndq0v7g+rOw92rBZOy0+Ut8ck6JLxqQoI8XFsCDQwxBcCC69Xp6n8r9h5uAxbTvkOa3hV5JsNmlQnyiNTInR8L4xGpoUraFJ0UpPiJIzxGFC5b3DkZLj2nzgmLbsL9bmA8e0I9fbYPjPYbdpfGqsZg1P1LlnJGlUf8IK0JMRXAguOEUgYOhAcYW2H/FqR65X23O92n7E2+iZGal2qCktPlJDk6I1JClaQxOjNTgxSqnxkUqMdvIm2kqBgKFDx45re27t874j16tthz2N7hLdPzZCM4claNbwRJ09NEHuiM5ZkBBA8CO4EFzQSkVlPu3ILdWOXK92F5RqT0GZ9hSUyVtZ0+TPhIfalRYfqbT4SKXGR2pgfKTS+kSqf2ykkl3hckWE9LpgU1nt18HiCmUXlutAUYWyi8q1K69UX+eVqsx3+nPpsNt0Zj+XJqTFaeKgOJ01ME4pbtbqAXorggvBBR1gGIaOlvm0p6BMe08EmT1Hy3SgqEJHSo6fNqPpVOGhdqW4I9TX5TzxNVzJLqeSXOGKjwpTfFSY4iLDFBcZaolm4UDAUGG5T/ken/K8lcrzVirfU6l8b6UOlxzX/sJy5Xor1dT/HGEOu4b1jdbIFJdGptT2F41NdSsyjN4iALV6zV5FQFew2WxKiglXUky4zh6S0OB7VTUBHSk5rpziCh0ortDB4grlFNX+OddzXCUV1aqsDii7sPy01YAb444IrQ8zsRGhinKGKMoZopjwEEWFhSg6PETRToeinCGKDHMo1GGvv4Q57Apx2Or/bLNJhiEFDEOGTnw1DBmG5DcM+aoDOl7tV2X9JaDKar+OV/tVVlmjkuPVKqmolud41Ymv1So5Xq1j5VWqaSmtSYpxhmhQQpQG9olUekKUhiTWhpXBiVFMVQbQZQguQDPCQuwalBClQQmNT8OtrPYrz3PirIS3UrmeytrrnkoVlvlUXFGl4vLaYCDVrlHjOV7dqpBjJrtNSoh2KtldG+iS3U4lu8KV7I5QekKkBvaJUp+osF43RAbAfAQXoAPCQx3NBps6Nf5A/dmM4hMXz/FqlflqVOarUbmvRmU+/0l/rlFltV/VfkPV/kDtpSagqhPXq2oCMmTIbrPJJtV+tdWePbLbaq+HhzrkDLUrItSh8FBH/dfwULuinSGKjQxTbGSo3BGhio0IlTsyVLERYYqLClVitNMSw1oAeh+CC9ANQhx2JUQ7lRDN8vQA0BF8pAIAAJZBcAEAAJZBcAEAAJZBcAEAAJZBcAEAAJZBcAEAAJZBcAEAAJZBcAEAAJZBcAEAAJZBcAEAAJZBcAEAAJZBcAEAAJZBcAEAAJZBcAEAAJYRYnYBHWEYhiTJ6/WaXAkAAGituvftuvfxtrB0cCktLZUkpaammlwJAABoq6KiIrnd7jb9jM1oT9wJEoFAQEeOHFFMTIxsNpvZ5ZjC6/UqNTVVBw8elMvlMrucoMPz0zyen6bx3DSP56d5PD/N83g8SktL07FjxxQbG9umn7X0GRe73a4BAwaYXUZQcLlcvDiawfPTPJ6fpvHcNI/np3k8P82z29veaktzLgAAsAyCCwAAsAyCi8U5nU4tXLhQTqfT7FKCEs9P83h+msZz0zyen+bx/DSvI8+PpZtzAQBA78IZFwAAYBkEFwAAYBkEFwAAYBkEFwAAYBkEFwt7+umnNWjQIIWHh2vKlCnauHGj2SUFjbVr12revHnq16+fbDab/v73v5tdUtBYvHixJk2apJiYGCUlJenyyy/Xzp07zS4raCxZskRjxoypXzhs2rRpevfdd80uK2g9+uijstlsuueee8wuJSg8/PDDstlsDS4jRowwu6ygcfjwYV1//fXq06ePIiIiNHr0aG3evLlN90FwsajXX39d9913nxYuXKjMzEyNHTtWF198sQoKCswuLSiUl5dr7Nixevrpp80uJeisWbNGCxYs0IYNG/Tee++purpaF110kcrLy80uLSgMGDBAjz76qLZs2aLNmzfr/PPP12WXXaavvvrK7NKCzqZNm/Tss89qzJgxZpcSVM4880zl5ubWX9atW2d2SUHh2LFjmj59ukJDQ/Xuu+9q+/bt+u1vf6u4uLi23ZEBS5o8ebKxYMGC+ut+v9/o16+fsXjxYhOrCk6SjBUrVphdRtAqKCgwJBlr1qwxu5SgFRcXZ7zwwgtmlxFUSktLjWHDhhnvvfeeMWvWLOPuu+82u6SgsHDhQmPs2LFmlxGU7r//fmPGjBkdvh/OuFhQVVWVtmzZotmzZ9ffZrfbNXv2bK1fv97EymBFHo9HkhQfH29yJcHH7/frtddeU3l5uaZNm2Z2OUFlwYIFuuSSSxr8P4Rau3fvVr9+/TR48GBdd911ysnJMbukoPCPf/xDEydO1JVXXqmkpCSNHz9ezz//fJvvh+BiQYWFhfL7/erbt2+D2/v27au8vDyTqoIVBQIB3XPPPZo+fbpGjRpldjlBY9u2bYqOjpbT6dRtt92mFStWKCMjw+yygsZrr72mzMxMLV682OxSgs6UKVP00ksvadWqVVqyZImys7M1c+ZMlZaWml2a6fbt26clS5Zo2LBhWr16tW6//Xbdddddevnll9t0P5beHRpAxyxYsEBZWVmMwZ/ijDPO0NatW+XxePTmm29q/vz5WrNmDeFF0sGDB3X33XfrvffeU3h4uNnlBJ25c+fW/3nMmDGaMmWKBg4cqDfeeEO33HKLiZWZLxAIaOLEiXrkkUckSePHj1dWVpaeeeYZzZ8/v9X3wxkXC0pISJDD4VB+fn6D2/Pz85WcnGxSVbCaO+64Q++8844++ugjDRgwwOxygkpYWJiGDh2qs846S4sXL9bYsWP1+9//3uyygsKWLVtUUFCgCRMmKCQkRCEhIVqzZo2efPJJhYSEyO/3m11iUImNjdXw4cO1Z88es0sxXUpKymnhf+TIkW0eSiO4WFBYWJjOOussffDBB/W3BQIBffDBB4zDo0WGYeiOO+7QihUr9OGHHyo9Pd3skoJeIBCQz+czu4ygcMEFF2jbtm3aunVr/WXixIm67rrrtHXrVjkcDrNLDCplZWXau3evUlJSzC7FdNOnTz9t6YVdu3Zp4MCBbbofhoos6r777tP8+fM1ceJETZ48WU888YTKy8t10003mV1aUCgrK2vwCSc7O1tbt25VfHy80tLSTKzMfAsWLNDy5cu1cuVKxcTE1PdFud1uRUREmFyd+R588EHNnTtXaWlpKi0t1fLly/Xxxx9r9erVZpcWFGJiYk7rh4qKilKfPn3ok5L0k5/8RPPmzdPAgQN15MgRLVy4UA6HQ9dee63ZpZnu3nvv1dlnn61HHnlEV111lTZu3KjnnntOzz33XNvuqOMTnGCWP/zhD0ZaWpoRFhZmTJ482diwYYPZJQWNjz76yJB02mX+/Plml2a6xp4XScbSpUvNLi0o3HzzzcbAgQONsLAwIzEx0bjggguMf//732aXFdSYDv1fV199tZGSkmKEhYUZ/fv3N66++mpjz549ZpcVNN5++21j1KhRhtPpNEaMGGE899xzbb4Pm2EYRicGKgAAgC5DjwsAALAMggsAALAMggsAALAMggsAALAMggsAALAMggsAALAMggsAALAMggsAALAMggsAALAMgguAoPXcc8/p3HPPlcvlks1mU0lJidklATAZwQVA0KqoqNCcOXP0s5/9zOxSAAQJggsA0wwaNEhPPPFEg9vGjRunhx9+WJJ0zz336IEHHtDUqVO7vzgAQYngAgAALIPgAgAALIPgAgAALIPgAsA0drtdhmE0uK26utqkagBYAcEFgGkSExOVm5tbf93r9So7O9vEigAEuxCzCwDQe51//vl66aWXNG/ePMXGxuqhhx6Sw+Go/35eXp7y8vK0Z88eSdK2bdsUExOjtLQ0xcfHm1U2ABMRXACY5sEHH1R2dra++c1vyu1265e//GWDMy7PPPOMFi1aVH/9nHPOkSQtXbpUN954Y3eXCyAI2IxTB5gBAACCFD0uAADAMgguAADAMgguAADAMgguAADAMgguAADAMgguAADAMgguAADAMgguAADAMgguAADAMgguAADAMgguAADAMv4/SbOrc7H+WUgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 結構めんどくさい．ギブアップ\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def _convert_lognormal_params(mu, cov):\n",
    "    \"\"\"Helper function to convert lognormal means and covariances to normal\n",
    "    means and covariances.\n",
    "    \"\"\"\n",
    "    sigma = np.diag(cov)\n",
    "    norm_sigma = np.log(1 + (sigma / (mu ** 2)))\n",
    "    norm_mu = np.log(mu) - 0.5 * norm_sigma\n",
    "    i, j = np.indices(cov.shape)\n",
    "    norm_cov = np.log(1 + (cov[i, j] / np.exp(norm_mu[i] + norm_mu[j] + 0.5 * (norm_sigma[i] + norm_sigma[j]))))\n",
    "    return norm_mu, norm_cov\n",
    "\n",
    "\n",
    "def psd_sqrt(P):\n",
    "    evalues, evectors = np.linalg.eig(P)\n",
    "    assert (evalues >= 0).all()\n",
    "    sqrt_matrix = evectors * np.sqrt(evalues) @ np.linalg.inv(evectors)\n",
    "    return sqrt_matrix\n",
    "\n",
    "\n",
    "mu = np.array([0.9, 0.7])\n",
    "Sigma = np.array([\n",
    "    [0.6, -0.4],\n",
    "    [-0.3, 0.1]\n",
    "])\n",
    "norm_mu, norm_Sigma = _convert_lognormal_params(mu, Sigma)\n",
    "\n",
    "N = 50\n",
    "u = np.random.multivariate_normal(norm_mu, norm_Sigma, size=N)\n",
    "u = np.exp(u)\n",
    "\n",
    "# 本当にこの式で合ってるかは微妙...\n",
    "mu_est = np.log(u).mean(axis=0)\n",
    "Sigma_est = np.zeros_like(Sigma)\n",
    "for u_ in u:\n",
    "    Sigma_est += (np.log(u_) - mu_est).reshape(2, 1) @ (np.log(u_) - mu_est).reshape(1, 2)\n",
    "Sigma_est = Sigma_est / N\n",
    "\n",
    "mean_est = np.exp(mu_est + np.diag(Sigma_est) / 2)\n",
    "Sigma_est = (mean_est.reshape(2, 1) @ mean_est.reshape(1, 2)) * (np.exp(Sigma_est) - 1)\n",
    "\n",
    "print(mean_est)\n",
    "print(Sigma_est)\n",
    "\n",
    "A_st = psd_sqrt(np.linalg.inv(Sigma_est))\n",
    "b_st = - A_st @ mean_est\n",
    "print(A_st)\n",
    "print(b_st)\n",
    "\n",
    "plt.scatter(u[:, 0], u[:, 1])\n",
    "plt.xlabel(\"u1\")\n",
    "plt.ylabel(\"u2\")\n",
    "\n",
    "\n",
    "k = np.random.choice(4, size=2) + 2\n",
    "p = k + np.random.choice(2, size=2) + 2\n",
    "\n",
    "def empirical_constraint_vio_prob(rho):\n",
    "\n",
    "\n",
    "\n",
    "theta = np.linspace(0, 2 * np.pi, 100)\n",
    "x = rho * np.cos(theta)\n",
    "y = rho * np.sin(theta)\n",
    "U = np.vstack((x.reshape(1, -1), y.reshape(1, -1)))\n",
    "U = mean_est.reshape(2, 1) + Sigma_est @ U\n",
    "\n",
    "plt.plot(U[0], U[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}