notebooks/CVX_minimum_volume_ellipsoids.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimum Volume Ellipsoids\n",
    "\n",
    "参考\n",
    "* [Minimum Volume Ellipsoids](https://www.amazon.co.jp/-/en/Michael-J-Todd/dp/1611974372)\n",
    "* [Determinants and Volumes](https://textbooks.math.gatech.edu/ila/determinants-volumes.html#:~:text=The%20theorem%20on%20determinants%20and,a%20number%20A%20a%20B%20.)\n",
    "\n",
    "## モチベーション\n",
    "\n",
    "まず適当な$m$点の凸包をイメージしてみましょう。\n",
    "例えば$X\\in \\mathbb{R}^n$を\n",
    "\n",
    "$$X := \\operatorname{conv} \\{x_1, x_2, \\dots, x_m\\}$$\n",
    "\n",
    "とします。\n",
    "これは$m$個の線形不等式についての多面体として表記することも可能です：\n",
    "$$X = \\{x\\in \\mathbb{R}^n: Ax \\leq b\\}$$\n",
    "\n",
    "このような多面体$X$上で、適当な線形関数$c^Tx$を最適化することを考えてみましょう。\n",
    "\n",
    "1. 前者の場合：（確か）線形関数は多面体の頂点でその最大（最小）値を取るため、シンプルに$x_1, \\dots, x_m$について$c^Tx_i$を調べれば、最適化できます。\n",
    "2. 後者の場合：$m$個の不等式についての線型計画法を解けばわかります。\n",
    "\n",
    "しかし、どちらも$m$について依存した計算量がかかってしまいます。なんとか$n$にだけ依存した計算量に変えられないでしょうか？\n",
    "\n",
    "一つの候補として、軸にそった直方体を使う方法があります。つまり、$X$を囲むような箱を考えます。\n",
    "これは良さそうに見えますが、この方法は軸の作り方によってでてくる直方体の形が全然変わってしまい、特定のデータセットをよく表すことができません。\n",
    "例えば次のように、「データを45度回転させたデータ上の直方体」$\\neq$「元のデータの直方体を45度回転させたもの」になってしまいます（アフィン不変性が保証されていません）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9gAAAH/CAYAAACsHHXWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAnEElEQVR4nO3dbZCV9Xn48eucBRYRdleSZdVxFxds4wB/bceHCDEVEmy0M0l5gUmnbYLWMOpgZhw6bSUzyjiZDG102k6No05owUni5EGLtmmMz2JS0YgJVVDMCOgiypN0d5HGXdhz/18gi4TdZYHfOffuns9n5n5xlgP35b3Ctd89D1vIsiwLAAAA4KQU8x4AAAAARgKBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAmUNbDvvvvuOO+886Kuri7q6upi5syZ8cgjj5TzlABABdn1AHBYIcuyrFx/+H/+539GTU1N/N7v/V5kWRb33Xdf3H777fHrX/86pk+fXq7TAgAVYtcDwGFlDey+TJw4MW6//fa49tprK3laAKBC7HoAqtWoSp2op6cnfvzjH8e+ffti5syZfd6nq6srurq6em+XSqXYs2dPfOxjH4tCoVCpUQGgX1mWxd69e+PMM8+MYtFbmXyUXQ/ASHBSuz4rs5dffjk79dRTs5qamqy+vj77r//6r37vu3Tp0iwiHA6Hw+EY8sfWrVvLvUKHDbve4XA4HCPxOJFdX/aniHd3d0dbW1t0dHTEAw88EMuXL4/Vq1fHtGnTjrrv735Xu6OjI1paWmLr1q1RV1dXzjEBYFA6Ozujubk52tvbo76+Pu9xhgS7HoCR5GR2fcVfgz137tyYOnVq3Hvvvce8b2dnZ9TX10dHR4elC8CQYDcdm10PwHB2Mrup4i8eK5VKR3znGgAYWex6AKpVWd/kbMmSJXHllVdGS0tL7N27N+6///545pln4tFHHy3naQGACrHrAeCwsgb2zp074ytf+Uq8++67UV9fH+edd148+uijcfnll5fztABAhdj1AHBYWQP7X//1X8v5xwMAObPrAeAwP8ATAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAFDVPtjfE7v2dsUH+3vyHoXw+WB4G5X3AAAAkIcX39wTy3++OR5/dUeUsohiIeLyaU2x8NNT4sKzJ+Y9XtXx+WAk8Ag2AABV57vPvxVfvGdNPPHazihlBz9WyiKeeG1nXHXPmvje82/lO2CV8flgpBDYAABUlRff3BO3PrQ+sojoOVRzH+opZZFFxC0PrY+1b+7JZb5q4/PBSCKwAQCoKst/vjmKxcKA9ykWC7H8F1sqNFF18/lgJPEabAAAqsYH+3t6X+M7kJ5SFo+8/E6cNXl+RM/+ygxXjWpGR82X/iUKxYEf9+spZfHYhu3xwf6eGDu6pkLDwfET2AAAVI29Hxw4ZlwfUigW493d7VH6v/ayzlTNiuMaovkYcX1IKTv4+RPYDGUCGwCAqjFh7KgoFmJQkZ2VStHQUBcx4ZTyD1atakZHViod8xHsiIPvKj5hrHxhaPN/KAAAVWPs6Jq4fFpTPPHazqPeUOujioWIGa0T4+qnnq/gdNVpxVObY0Nb+4Df9KgpFuLyaU0evWbI8yZnQFX7YH9P7NrbFR/s78l7FAAq5KufnhKlYzyEXcoiLps+qUITVbfLpk865jMKSqUsvnppa2UGgpPgEWygKr345p5Y/vPNvW90UyxEXD6tKRZ+ekpcePbEvMcDoIwuOntifGPejLjlofVRLBaOeCT70NPH589sjtam8TlOWT2mNI2P+TOb44E1W496+n5NsRClUhbfmDfDfmZYKOsj2MuWLYuLLrooJkyYEJMmTYp58+bF66+/Xs5TAhzTd59/K754z5p44rWdvUu8lEU88drOuOqeNfG959/Kd0AYRux6hqu/vGRy/Pj6mXH5tKbISqWIOPia6xktDfG1P/n9mHVuY84TVpdZ5zbG1/7k92NGS8MRn4/LpzXFj6+fGX95yeScJ4TBKWRZNsj3UTx+V1xxRfzZn/1ZXHTRRXHgwIH4+te/HuvXr49XX301Tj311GP+/s7Ozqivr4+Ojo6oq6sr15hAFXnxzT3xxXvWxED/8BUi4sfXz/SdcvpkNx3JrmckOGtya7y7uz0aGuriO15znbuFn7kk2ts744yPN8Tbb/nZ11Teyeymsj5F/Gc/+9kRt1euXBmTJk2Kl156Kf7oj/6onKcG6NPyn28+6umAv6tYLMTyX2wR2DAIdj0jQs/+gz+Ky7uFDw2HPh89x/4mHQw1FX0NdkdHR0RETJzY9xetXV1d0dXV1Xu7s7OzInNBNbjwwgtj+/bteY+Rr5rRUfOlfznmjwLpKWXxyMvvxFmT50f07K/QcEPT6aefHmvXrs17DIYRux6AalaxwC6VSnHTTTfFpz71qZgxY0af91m2bFncdtttlRoJqsr27dtj27ZteY+Rq+K4hmgexM/ZjIgoFIvx7u72g99BBwbFrgeg2lUssBctWhTr16+PX/ziF/3eZ8mSJbF48eLe252dndHc3FyJ8aBqFIvFaGhsynuMfNSMjqxUOuYj2BEH31iloaGuap8u2L5rR5Q+fJMZGCy7HoBqV5HAvvHGG+MnP/lJPPvss3HWWWf1e7/a2tqora2txEhQtRoam+I7q1/Ke4zcrHhqc2xoax/w520WCxEzWifG1VX8RjcLL7sg9ux4N+8xGEbsegAo84/pyrIsbrzxxli1alU89dRT0drqh8MD+bps+qQB4zri4I/sumz6pMoMBMOcXQ8Ah5U1sBctWhTf+9734v77748JEybE9u3bY/v27fHb3/62nKcF6NeUpvExf+bBp6MWC0f+2qHb82c2R2vT+ApPBsOTXQ8Ah5X1KeJ33313RETMnj37iI+vWLEirr766nKeGqBfs85tjDNOOyVWb9gZr7S1R5ZFFAoRM1oa4rLpk8Q1HAe7HgAOK2tgZ9kxnocJkJPWpvHR2jQ+ug+Uomt/T9SOrokxo8r6pB4Ykex6ADisoj8HG2CoGTOqKKwBAEjCV5UAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAEaA7gOl2Pvb/dF9oJT3KFC1RuU9AAAAcOI273g/Vm/YGevb2iPLIgqFiBktDTF7+qRobRqf93hQVQQ2AAAMU/+9cVc8uGZrFAsRWXbwY1kWsaGtPV55qz3mz2yOWec25jskVBFPEQcAgGFo847348E1WyMiopQd+WuHbj+wZmts2fF+hSeD6iWwAQBgGFq9YWcUCwPfp1g4eD+gMgQ2AAAMM90HSrG+rf2oR65/VymLeKWt3RufQYUIbAAAGGa69vf0vub6WLLs4P2B8hPYAAAwzNSOronCMZ4efkihcPD+QPkJbAAAGGbGjCrGjJaGQb0G+/+1NMSYUb7sh0rwNw0AAIahy6ZPGtRrsC+bPqkyAwECGwAAhqMpTeNj/szmiIijHsk+dHv+zOZobRpf4cmgeo3KewAAAODEzDq3Mc447ZRYvWFnvNLWHll28DXXM1oa4rLpk8Q1VJjABgCAYay1aXy0No2P7gOl6NrfE7Wja7zmGnIisAEAYAQYM6oorCFn/gYCAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYQFXrPlCKvb/dH90HSnmPMrTUjI7iuIaImtF5TwIAMGyMynsAgDxs3vF+rN6wM9a3tUeWRRQKETNaGmL29EnR2jQ+7/Fyc+i6jP+Lb8eEYjGyUimu++7aWPjpKXHh2RPzHg8AYEgr6yPYzz77bHz+85+PM888MwqFQjz00EPlPB3AoPz3xl3x7Z/+JjZ8GNcREVkWsaGtPe786W/iuY278h0wJx+9LoXiwfVQKBbjidd2xlX3rInvPf9WzhMyFNn1AHBYWQN73759cf7558ddd91VztMADNrmHe/Hg2u2RkREKTvy1w7dfmDN1tiy4/0KT5avga5LTymLLCJueWh9rH1zT+WHY0iz6wHgsLI+RfzKK6+MK6+8spynADguqzfsjGLh6Ij8qGLh4P2q6anig7ouxUIs/8UWTxXnCHY9ABw2pF6D3dXVFV1dXb23Ozs7c5wGGGm6D5R6X3M9kFIW8Upbe3QfKMWYUSP/vSAHe116Slk8tmF7fLC/J8aOrqnMcIw4dj0AI9mQ+spx2bJlUV9f33s0NzfnPRIwgnTt7zlmRB6SZQfvXw2O57qUsoi9Hxwo70CMaHY9ACPZkArsJUuWREdHR++xdevWvEcCRpDa0TVRKAzuvoXCwftXg+O5LsVCxISxQ+rJTwwzdj0AI9mQ+iqptrY2amtr8x4DGKHGjCrGjJaG2NDWfszXYM9oaaiKp4dHDP661BQLcfm0Jk8P56TY9QCMZNXx1SPAhy6bPmnAiIw4+DToy6ZPqsxAQ8Sgrkspi69e2lqZgQAAhqGyBvb7778f69ati3Xr1kVExJYtW2LdunXR1tZWztMC9GtK0/iYP/Pgaz6Lv/O06EO3589srqp3EI8Y+LrUFAtRiIhvzJvhHcQ5il0PAIeV9Snia9eujTlz5vTeXrx4cURELFiwIFauXFnOUwP0a9a5jXHGaafE6g0745UP3z278OHTwi+bPqnq4vqQj16X/9myJwrFYmSlUlw+48z46qWt4po+2fUAcFhZA3v27NmRDfataQEqqLVpfLQ2jY/uA6Xo2t8TtaNrquY11wM5dF0WfubPo729M874eEPc860teY/FEGbXA8BhQ+pNzgAqbcyoorDuS8/+KP1fe0TPqXlPAgAwbPiqEgAAABIQ2AAAAJCAwAYAAIAEBDYAAAAkILABAAAgAYENAAAACQhsAAAASEBgAwAAQAICGwAAABIQ2AAAAJCAwAYAAIAEBDYAAAAkILABAAAgAYENAAAACQhsAAAASEBgAwAAQAICGwAAABIQ2AAAAJCAwAYAAIAEBDYAAAAkILABAAAgAYENAAAACQhsAAAASEBgAwAAQAICGwAAABIQ2AAAAJCAwAYAAIAEBDYAAAAkILABAAAgAYENAAAACQhsAAAASEBgAwAAQAICGwAAABIQ2AAAAJCAwAYAAIAEBDYAAAAkILABAAAgAYENAAAACQhsAAAASEBgAwAAQAICGwAAABIQ2AAAAJCAwAYAAIAEBDYAAAAkILABAAAgAYENAAAACQhsAAAASEBgAwAAQAICGwAAABIQ2AAAAJCAwAYAAIAEBDYAAAAkILABAAAgAYENAAAACQhsAAAASEBgAwAAQAICGwAAABIQ2AAAAJCAwAYAAIAEBDYAAAAkILABAAAgAYENAAAACQhsAAAASEBgQ7WoGR3FcQ0RNaPzngQAKIPuA6XY+9v90X2glPcoULVG5T0AUF4vvrknlv98c9R86V+iuViMrFSKFU9tjtnTJ0Vr0/i8xwMATtLmHe/H6g07Y31be2RZRKEQMaOlwa6HHFTkEey77rorzj777Bg7dmx88pOfjF/+8peVOC1Uve8+/1Z88Z418cRrO6NQPPjXvVAsxoa29rjzp7+J5zbuynlCYKSw6yEf/71xV3z7p7+JDR/GdUREloVdDzkpe2D/8Ic/jMWLF8fSpUvjV7/6VZx//vnxuc99Lnbu3FnuU0NVe/HNPXHrQ+sji4ieUnbErx26+cCarbFlx/uVHw4YUex6yMfmHe/Hg2u2RsTh3X6IXQ/5KHtg/+M//mMsXLgwrrnmmpg2bVrcc889MW7cuPi3f/u3cp8aqtryn2+OYrEw4H2KhYjVG3wBDJwcux7ysXrDzjjGqrfrocLK+hrs7u7ueOmll2LJkiW9HysWizF37txYs2bNUffv6uqKrq6u3tudnZ3lHA9GrA/298Tjr+446rvZv6uURfzPlj2x8DN/HtGzvzLDMSy079qR9wgME3Y9I0n7rh2x8LIL8h5jcGpGx/i/+HbvS8D6Mxx3vR3EcFbWwN69e3f09PREU1PTER9vamqKjRs3HnX/ZcuWxW233VbOkaAq7P3gwDHj+pBCsRjt7Z1R+r/2ss4EjEx2PSNJqVSKPTvezXuMQSmOa4gJx4jrQ+x6qJwh9S7iS5YsicWLF/fe7uzsjObm5hwnguFpwthRUSwc/XqsvmSlUpzx8YaInlPLPhfDz+mnn573CIwwdj1D0bD8t65mdGSl0jEfwY4Yvrt+WH5eqHplDeyPf/zjUVNTEzt2HPk0jx07dvT5F6a2tjZqa2vLORJUhbGja+LyaU3xxGs7j3qDs4+qKRbi8hlnxj3f2lLB6YCRxK5nJFi7dm3eI5yQ67671q6HIaasb3I2ZsyYuOCCC+LJJ5/s/VipVIonn3wyZs6cWc5TQ9X76qenROkYD2GXSll89dLWCk0EjER2PeTHroehp+zvIr548eL4zne+E/fdd1+89tprccMNN8S+ffvimmuuKfepoapddPbE+Ma8GVGIg9+9/qiaYiEKEfGNeTPiwrMn5jIfMHLY9ZAPux6GnrK/BvtLX/pS7Nq1K2699dbYvn17/MEf/EH87Gc/O+rNUID0/vKSyXHu6RNi+S+2xGMbtkcpO/jjOi6f1hRfvbTVwgWSsOshP3Y9DC2FLMsG+V7DldfZ2Rn19fXR0dERdXV1eY8Dw9oH+3ti7wcHYsLYUTF2dE3e48CwZTel5XpCOnY9pHEyu2lIvYs4UD5jR9dYtgAwgtn1kL+yvwYbAAAAqoHABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkEDZAvub3/xmzJo1K8aNGxcNDQ3lOg0AkCP7HgAOK1tgd3d3x1VXXRU33HBDuU4BAOTMvgeAw0aV6w++7bbbIiJi5cqV5ToFAJAz+x4ADitbYJ+Irq6u6Orq6r3d2dmZ4zQAQGp2PQAj2ZB6k7Nly5ZFfX1979Hc3Jz3SABAQnY9ACPZcQX2zTffHIVCYcBj48aNJzzMkiVLoqOjo/fYunXrCf9ZAMCJKee+t+sBGMmO6ynif/3Xfx1XX331gPeZMmXKCQ9TW1sbtbW1J/z7AYCTV859b9cDMJIdV2A3NjZGY2NjuWYBAIYA+x4ATkzZ3uSsra0t9uzZE21tbdHT0xPr1q2LiIhzzjknxo8fX67TAgAVZN8DwGFlC+xbb7017rvvvt7bf/iHfxgREU8//XTMnj27XKcFACrIvgeAwwpZlmV5D9Gfzs7OqK+vj46Ojqirq8t7HACwmxJzPQEYak5mNw2pH9MFAAAAw5XABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgATKFthvvvlmXHvttdHa2hqnnHJKTJ06NZYuXRrd3d3lOiUAUEF2PQAcaVS5/uCNGzdGqVSKe++9N84555xYv359LFy4MPbt2xd33HFHuU4LAFSIXQ8ARypkWZZV6mS333573H333bF58+ZB3b+zszPq6+ujo6Mj6urqyjwdAByb3TQwux6A4e5kdlPZHsHuS0dHR0ycOLHfX+/q6oqurq7e252dnZUYCwBIxK4HoJpV7E3O3njjjbjzzjvjuuuu6/c+y5Yti/r6+t6jubm5UuMBACfJrgeg2h13YN98881RKBQGPDZu3HjE79m2bVtcccUVcdVVV8XChQv7/bOXLFkSHR0dvcfWrVuP/78IADgpdj0AnJjjfg32rl274r333hvwPlOmTIkxY8ZERMQ777wTs2fPjksuuSRWrlwZxeLgm97rsgAYaqphN9n1AFSzir4Gu7GxMRobGwd1323btsWcOXPiggsuiBUrVhzXwgUA8mHXA8CJKdubnG3bti1mz54dkydPjjvuuCN27drV+2unn356uU4LAFSIXQ8ARypbYD/++OPxxhtvxBtvvBFnnXXWEb9WwZ8MBgCUiV0PAEcq2/O4rr766siyrM8DABj+7HoAOJIXSgEAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASENgAAACQgMAGAACABAQ2AAAAJCCwAQAAIAGBDQAAAAkIbAAAAEhAYAMAAEACAhsAAAASKGtgf+ELX4iWlpYYO3ZsnHHGGfHlL3853nnnnXKeEgCoILseAA4ra2DPmTMnfvSjH8Xrr78eDz74YGzatCnmz59fzlMCABVk1wPAYYUsy7JKnew//uM/Yt68edHV1RWjR48+5v07Ozujvr4+Ojo6oq6urgITAsDA7KaB2fUADHcns5tGlWmmo+zZsye+//3vx6xZs/pduF1dXdHV1dV7u6OjIyIO/gcCwFBwaCdV8PvTw4ZdD8BIcFK7Piuzv/3bv83GjRuXRUR2ySWXZLt37+73vkuXLs0iwuFwOByOIX9s2rSp3Ct02LDrHQ6HwzESjxPZ9cf9FPGbb745/uEf/mHA+7z22mtx7rnnRkTE7t27Y8+ePfHWW2/FbbfdFvX19fGTn/wkCoXCUb/vd7+r3d7eHpMnT462traor68/njFHvM7Ozmhubo6tW7d6St1HuC79c2365rr0zXXpX0dHR7S0tMT//u//RkNDQ97jlIVdPzT4e9g316V/rk3fXJf+uTZ9O5ldf9yBvWvXrnjvvfcGvM+UKVNizJgxR3387bffjubm5njuuedi5syZxzyX12X1z7Xpm+vSP9emb65L31yX/lXDtbHrhwbXpm+uS/9cm765Lv1zbfpW0ddgNzY2RmNj4/H+toiIKJVKERFHfOcaABha7HoAODFle5OzF154IV588cW49NJL47TTTotNmzbFLbfcElOnTh3Ud7QBgKHNrgeAI5Xt52CPGzcu/v3f/z0++9nPxic+8Ym49tpr47zzzovVq1dHbW3toP6M2traWLp06aDvX01cm765Lv1zbfrmuvTNdemfa3OYXV9erk3fXJf+uTZ9c13659r07WSuS0V/DjYAAACMVGV7BBsAAACqicAGAACABAQ2AAAAJCCwAQAAIIFhFdhf+MIXoqWlJcaOHRtnnHFGfPnLX4533nkn77Fy9eabb8a1114bra2tccopp8TUqVNj6dKl0d3dnfdoufvmN78Zs2bNinHjxkVDQ0Pe4+TqrrvuirPPPjvGjh0bn/zkJ+OXv/xl3iPl7tlnn43Pf/7zceaZZ0ahUIiHHnoo75GGhGXLlsVFF10UEyZMiEmTJsW8efPi9ddfz3usIeHuu++O8847L+rq6qKuri5mzpwZjzzySN5jjTh2/dHs+oHZ9wfZ9Uez6/tm1/cvxa4fVoE9Z86c+NGPfhSvv/56PPjgg7Fp06aYP39+3mPlauPGjVEqleLee++NDRs2xD/90z/FPffcE1//+tfzHi133d3dcdVVV8UNN9yQ9yi5+uEPfxiLFy+OpUuXxq9+9as4//zz43Of+1zs3Lkz79FytW/fvjj//PPjrrvuynuUIWX16tWxaNGieP755+Pxxx+P/fv3xx//8R/Hvn378h4td2eddVb8/d//fbz00kuxdu3a+MxnPhN/+qd/Ghs2bMh7tBHFrj+aXT8w+96u749d3ze7vn9Jdn02jD388MNZoVDIuru78x5lSPnWt76Vtba25j3GkLFixYqsvr4+7zFyc/HFF2eLFi3qvd3T05OdeeaZ2bJly3KcamiJiGzVqlV5jzEk7dy5M4uIbPXq1XmPMiSddtpp2fLly/MeY0Sz6/tm1x+tmve9XX9sdn3/7PqBHe+uH1aPYH/Unj174vvf/37MmjUrRo8enfc4Q0pHR0dMnDgx7zEYArq7u+Oll16KuXPn9n6sWCzG3LlzY82aNTlOxnDR0dEREeHflN/R09MTP/jBD2Lfvn0xc+bMvMcZsez6/tn1HGLXc7Ls+r6d6K4fdoH9d3/3d3HqqafGxz72sWhra4uHH34475GGlDfeeCPuvPPOuO666/IehSFg9+7d0dPTE01NTUd8vKmpKbZv357TVAwXpVIpbrrppvjUpz4VM2bMyHucIeGVV16J8ePHR21tbVx//fWxatWqmDZtWt5jjTh2/cDsej7Krudk2PVHO9ldn3tg33zzzVEoFAY8Nm7c2Hv/v/mbv4lf//rX8dhjj0VNTU185StfiSzLcvwvKI/jvS4REdu2bYsrrrgirrrqqli4cGFOk5fXiVwX4MQsWrQo1q9fHz/4wQ/yHmXI+MQnPhHr1q2LF154IW644YZYsGBBvPrqq3mPNeTZ9X2z6/tn30Nl2PVHO9ldX8hy3li7du2K9957b8D7TJkyJcaMGXPUx99+++1obm6O5557bsQ9Re94r8s777wTs2fPjksuuSRWrlwZxWLu3zspixP5/2XlypVx0003RXt7e5mnG3q6u7tj3Lhx8cADD8S8efN6P75gwYJob2/3qNCHCoVCrFq16ohrVO1uvPHGePjhh+PZZ5+N1tbWvMcZsubOnRtTp06Ne++9N+9RhjS7vm92ff/s+8Gz6wfHrj+aXT84x7vrR5V5nmNqbGyMxsbGE/q9pVIpIiK6urpSjjQkHM912bZtW8yZMycuuOCCWLFixYheuCfz/0s1GjNmTFxwwQXx5JNP9i6UUqkUTz75ZNx44435DseQlGVZfO1rX4tVq1bFM888Y+EeQ6lUGpE7KDW7vm92ff/s+8Gz6zledv3xOd5dn3tgD9YLL7wQL774Ylx66aVx2mmnxaZNm+KWW26JqVOnjrjvaB+Pbdu2xezZs2Py5Mlxxx13xK5du3p/7fTTT89xsvy1tbXFnj17oq2tLXp6emLdunUREXHOOefE+PHj8x2ughYvXhwLFiyICy+8MC6++OL453/+59i3b19cc801eY+Wq/fffz/eeOON3ttbtmyJdevWxcSJE6OlpSXHyfK1aNGiuP/+++Phhx+OCRMm9L5+r76+Pk455ZScp8vXkiVL4sorr4yWlpbYu3dv3H///fHMM8/Eo48+mvdoI4Zd3ze7fmD2vV3fH7u+b3Z9/5Ls+rK8l3kZvPzyy9mcOXOyiRMnZrW1tdnZZ5+dXX/99dnbb7+d92i5WrFiRRYRfR7VbsGCBX1el6effjrv0SruzjvvzFpaWrIxY8ZkF198cfb888/nPVLunn766T7//1iwYEHeo+Wqv39PVqxYkfdoufurv/qrbPLkydmYMWOyxsbG7LOf/Wz22GOP5T3WiGLX982uH5h9f5BdfzS7vm92ff9S7PrcX4MNAAAAI8HIfgEPAAAAVIjABgAAgAQENgAAACQgsAEAACABgQ0AAAAJCGwAAABIQGADAABAAgIbAAAAEhDYAAAAkIDABgAAgAQENgAAACQgsAEAACCB/w8in6+udVbgjQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1200x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.path import Path\n",
    "import matplotlib.patches as patches\n",
    "\n",
    "fig = plt.figure(figsize=(12, 6))\n",
    "ax = fig.add_subplot(121)\n",
    "\n",
    "# 集合の描画\n",
    "x = [0, 1, 2, 1]\n",
    "y = [0, 2, 1, 1]\n",
    "\n",
    "# minimum volume enclosing box の描画\n",
    "\n",
    "x_min, x_max = min(x), max(x)\n",
    "y_min, y_max = min(y), max(y)\n",
    "\n",
    "verts = [\n",
    "    (x_min, y_min), # left, bottom\n",
    "    (x_min, y_max), # left, top\n",
    "    (x_max, y_max), # right, top\n",
    "    (x_max, y_min), # right, bottom\n",
    "    (0, 0), # ignored\n",
    "    ]\n",
    "\n",
    "codes = [Path.MOVETO,\n",
    "         Path.LINETO,\n",
    "         Path.LINETO,\n",
    "         Path.LINETO,\n",
    "         Path.CLOSEPOLY,\n",
    "         ]\n",
    "\n",
    "path = Path(verts, codes)\n",
    "patch = patches.PathPatch(path, facecolor='lightblue', lw=2)\n",
    "ax.add_patch(patch)\n",
    "plt.scatter(x, y, s=50)  # 集合\n",
    "plt.ylim(-3, 3)\n",
    "plt.xlim(-3, 3)\n",
    "\n",
    "ax = fig.add_subplot(122)\n",
    "\n",
    "# 45度回転させます\n",
    "\n",
    "R = np.array([[np.cos(np.pi / 4), -np.sin(np.pi / 4)],\n",
    "              [np.sin(np.pi / 4),  np.cos(np.pi / 4)]])\n",
    "x2, y2 = [], []\n",
    "for _x, _y in zip(x, y):\n",
    "    v = np.dot(R, np.array([_x, _y]))\n",
    "    x2.append(v[0])\n",
    "    y2.append(v[1])\n",
    "x, y = x2, y2\n",
    "\n",
    "# minimum volume enclosing box の描画\n",
    "\n",
    "x_min, x_max = min(x), max(x)\n",
    "y_min, y_max = min(y), max(y)\n",
    "\n",
    "verts = [\n",
    "    (x_min, y_min), # left, bottom\n",
    "    (x_min, y_max), # left, top\n",
    "    (x_max, y_max), # right, top\n",
    "    (x_max, y_min), # right, bottom\n",
    "    (0, 0), # ignored\n",
    "    ]\n",
    "\n",
    "codes = [Path.MOVETO,\n",
    "         Path.LINETO,\n",
    "         Path.LINETO,\n",
    "         Path.LINETO,\n",
    "         Path.CLOSEPOLY,\n",
    "         ]\n",
    "\n",
    "path = Path(verts, codes)\n",
    "patch = patches.PathPatch(path, facecolor='lightblue', lw=2)\n",
    "ax.add_patch(patch)\n",
    "plt.scatter(x, y, s=50)  # 集合\n",
    "plt.ylim(-3, 3)\n",
    "plt.xlim(-3, 3)\n",
    "plt.show()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 楕円\n",
    "\n",
    "直方体の代替案として楕円を考えてみましょう。\n",
    "\n",
    "$$\n",
    "\\mathcal{E}(H, \\bar{x}):=\\{x\\in \\mathbb{R}^n : (x-\\bar{x})^TH(x-\\bar{x}) \\leq n\\}\n",
    "$$\n",
    "\n",
    "ここで、\n",
    "* $\\bar{x}$は楕円の中心\n",
    "* $H$は$n$次の正則な対称行列です。つまり、$H^T=H$かつ$v^THv > 0$です。\n",
    "\n",
    "###  楕円の性質\n",
    "\n",
    "楕円は円をアフィン変換したものと同じです。\n",
    "楕円について別の解釈をするため、$H$をコレスキー分解してみます。つまり、$H=L^TL$となる下三角行列$L$で分解します（これが正則行列に対して存在することは帰納的に示せます）。\n",
    "よって、$\\|L^T(x-\\bar{x})\\|\\leq \\sqrt{n}$として楕円を表すこともできます。つまり、\n",
    "\n",
    "$$\n",
    "\\mathcal{E}(H, \\bar{x})=\\{x=\\bar{x} + (\\sqrt{n}L^{-T})z : z\\in \\mathbb{R}^{n}, \\|z\\|\\leq 1\\}\n",
    "$$\n",
    "\n",
    "これは$\\bar{x}$を中心とした半径$1$の円が$\\sqrt{n}L^{-T}$によってアフィン変換されていることを表しています。\n",
    "また、楕円は下三角行列の$n(n+1)/2$個の要素で表せることもわかります。\n",
    "\n",
    "さらに、このような楕円ではアフィン不変性が保証されることに注意しましょう。（教科書参照）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqEAAAH5CAYAAACxl+dZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABl30lEQVR4nO3de3wU5b0/8M9uIBtuu0kkZBMa5WaBCBIMJkapUggkwqnSY08BoVyKUKlgNZwK6Sk3UUOE46EChUpFQEDQHlEQmxYDyFEj0CBVIPIrNIJANggxuxAkhOz8/kh3ySZ7mb3M7Fw+79drX5rd2WUmIZ95mHm+38cgCIIAIiIiIiIZGaO9A0RERESkPxyEEhEREZHsOAglIiIiItlxEEpEREREsuMglIiIiIhkx0EoEREREcmOg1AiIiIikl2baO9ANDidTpw/fx6dOnWCwWCI9u4QkQoJgoDLly8jNTUVRqN2/j3PfCSicInNR10OQs+fP4+0tLRo7wYRacDXX3+N733ve9HejYhhPhJRpATKR10OQjt16gSg6ZtjNpujvDdEpEYOhwNpaWnuPNEK5iMRhUtsPupyEOq6xWQ2mxmyRBQWrd2yZj4SUaQEykftTGQiIiIiItXgIJSIiIiIZMdBKBERERHJjoNQIiIiIpIdB6FEREREJDsOQomIiIhIdhyEEhEREZHsOAglIiIiItlxEEpEREREsuMglIiIiIhkx0EoEREREcmOg1AiIiIikh0HoUREREQkuzbR3gEipWh0CjhYWYMLl6+hS6c4ZHVPRIzRIPp1ItIGf7/rzAGiyJF0ELp//34sXboU5eXlqKqqwvbt2zF69Gi/79m3bx8KCgpw7NgxpKWl4be//S0mT57ssc2qVauwdOlS2Gw2DBgwACtWrEBWVpZ0B0Ka4O/kUXK0Cot2HkeV/Zp7+xRLHBb8KB35/VICvk5E6uMtE3Yft/n8XQfAHCCKIEkHoXV1dRgwYAB+/vOf49///d8Dbl9ZWYlRo0bh8ccfx+bNm1FaWorHHnsMKSkpyMvLAwBs27YNBQUFWLNmDbKzs7F8+XLk5eXhxIkT6NKli5SHQwoX6iATAGZsOgyhxefZ7NcwY9NhTL+/O17ZX+nz9dUT7kJ+vxReISFSKLGDzfj2bVF7taHV+232a3h802Gvn908B4anW5kBREEwCILQ8twqzR9kMAS8Ejpnzhzs2rULR48edT83duxY1NbWoqSkBACQnZ2Nu+++GytXrgQAOJ1OpKWlYdasWZg7d66ofXE4HLBYLLDb7TCbzaEfFClGKINM16nB4uPE49rGYACcPn5LDACsljjMG9UXi3dV+LxCwgGq9mg1R7R2XN6ywddgM1SGf32mqY0RNke9+3lmAOmV2BxR1JzQsrIy5ObmejyXl5eHp556CgBw/fp1lJeXo7Cw0P260WhEbm4uysrKfH5ufX096utvBoPD4YjsjlNUlRyt8nsl09K+bavXALif83cyEgD4+2eaAKDKfg2/3PJZq9eaX0nd8fcq3sIjkpmvbIjkABRoyoFvfVxBZQYQ+aao6nibzYbk5GSP55KTk+FwOPDdd9/h4sWLaGxs9LqNzWbz+blFRUWwWCzuR1pamiT7T9JqdAooO3UJ7x45h7JTl9DoFNDoFLBo53Gfg0wBkT/hiOX68/+wv9Lj5APcPDmVHK2Kyr4RaVHzjPj45EUs3HHMazbIhRlA5J+iroRKpbCwEAUFBe6vHQ4HB6Iq4+t2+9i7b20V7mogoOkW3qKdxzE83QoAvFVHFAZvGaFkzTNgaJ9klJ/+lr//pDuKGoRarVZUV1d7PFddXQ2z2Yx27dohJiYGMTExXrexWq0+P9dkMsFkMkmyzyQ9f7fb/+eD/yfpnx1oTmg4XLfyV+75B7Ye+pq36ohC5CsjlM6VAfcUfYCaupt3bPj7T3qhqNvxOTk5KC0t9Xhu9+7dyMnJAQDExsYiMzPTYxun04nS0lL3NqRuLW+5X7/h9Hu7PVwGAAnt27r/v+VrADDtB92bBqM+Xg/X/3zwD96qIwqSKyu2f3YOv9n+haQDUIOP/4+U5gNQgL//pB+SXgm9cuUKTp486f66srISR44cQWJiIm699VYUFhbi3Llz2LhxIwDg8ccfx8qVK/HMM8/g5z//Ofbs2YM333wTu3btcn9GQUEBJk2ahEGDBiErKwvLly9HXV0dpkyZIuWhkAy83U5L7BCLmrrrIX+mq2r126sNMMBz4Oo6mRT9e38Arfv/WZtdjRh4a4LX1+eNSsfiXcdhs1+L6Emw5e163pojukmqW++ujGhZPW8N0Cf0u4ZG2K82RCwD+PtPeiHpIPRvf/sbfvjDH7q/ds3LnDRpEtavX4+qqiqcOXPG/Xr37t2xa9cuPP300/jd736H733ve/jjH//o7hEKAGPGjME333yD+fPnw2azISMjAyUlJa2KlUhdfN1OC2YAGs4gE4DfHn/5/VJ8vm40NrWA8vbnh3NSct2qO1hZg6zuiZwzSoTI3HoPNNj0lwXeXtt93BbxDODvP+mBbH1ClURrffDUrtEpYHDxnrCuajyde3vAeZVS9unzVTj10IAUvLK/EkDoJ6ep93XD+0dtnDOqMFrNESUfVySyArj5+xPJ5vLBZoBY/P0nNRKbIxyEKixktc7bQPBgZQ3Grf00pM9zNYv/aM5QANGtMPc1yPVX2R9qYZXrqFyrNZH8tJojSjwu1+/WxycvYuXek4Hf0IwrI5b9ZAAu1tVLmg3BZMAtHWJxKcSpRvz9J6VTZbN60jZfg7EH+/nubOCPK4gX/CjdfULJ6XlLuLsZshijweuf7+tWPgBsPXTG53xSf5X5nDNGehHO/M/mGXHf7Z0ju2NeBJMBmbcl4IGle/n7T7rGQSjJwl+bpXUffyXqMxI7tPWoIm05p1PJfJ2cFvwo3e9cMjGrNX36z0swGgycL0aaE+78TyVlhLcMiMTvP+eMkppxEEqSC7Sqkdj12T/89Q8119A5v18KVk+4y2vR1IP9rKIG6E9sPoza79hjkLTFX274YkBTR43fjuoLq6Wd4jMiEr//Hxy3oeDNI5wzSqrEOaEKmfOkZWWnLome8+mrwl3rc58iOVdWL9+zaNNqjijluILJDUDdf++lmCsPqPN7QdrAOaEUNS0D1Wb/TtT7vFWBKul2mpS83arL6p6IFEtc0D1IOV+M1K7RKeDjk98E9R41Z0Uov/+cM0pawEEoRZSvhvNi5KZb8ZtR6Zzb9C8xRoPPOWOBNJ8vFs1iLaJgBVuINPOHPXFfryTNZYW/3/9g54wyA0ipFLVsJ6mbq4ig5cnj2wBtSAxomsPkOonk9LwFD2d0RU7PWzR1UgmFa86Y1RLn8Xz8v5YaDcTmuOaxDGqjr4m3RArgK0O8ceXG08N7azYrfP3+Wy1x+Pl93UR9BjOAlIxXQikiAhUf+eKtzRJ58tbexekUMP7VAwHfu/i9Yx4dBViwQEoVTCGSnnLDV4u3g5U1ogqXmAGkZByEUkQcrKwRdfVCzW2WoqnlnLFGpyBqvmjz7zXQ1BJrxqbDLFggxRGbIYD+ciOcOePMAFIy3o6nkDQ6BY9bPGKLj+b92x14Y9o9+N3YDLwx7R58NGcogzAErvliwM2rQi7+rgu5TlaLdh7nbTkVWbVqFbp164a4uDhkZ2fj4MGDPrd9++23MWjQIMTHx6NDhw7IyMjA66+/LuPeBseVJX8+WiVq+5k/7MXcADOAtIFXQilo4RQfWc1xnCQfIb56DCYGWA6QBQvqsm3bNhQUFGDNmjXIzs7G8uXLkZeXhxMnTqBLly6ttk9MTMR//dd/oU+fPoiNjcV7772HKVOmoEuXLsjLy4vCEfgWympI9/XqrPlb8GIxA0jt2CdUQ/395OBrBZNA1dvN13jnCSSyvLXEevrNvwd83+/GZuDhjK4y7KE2yZUj2dnZuPvuu7Fy5UoAgNPpRFpaGmbNmoW5c+eK+oy77roLo0aNwuLFiwNuK9dxBbsaEjPEN2YAKQ37hFLEsfhImVrOFys7dUnU+7p0ivPaJJs/I+W4fv06ysvLUVhY6H7OaDQiNzcXZWVlAd8vCAL27NmDEydOoLi42Os29fX1qK+vd3/tcDjC3/EAgl0NiRniHzOA1IqDUBKNxUfqIKbJtdUSh2/r6jG4eA+X+1OwixcvorGxEcnJyR7PJycn48svv/T5Prvdjq5du6K+vh4xMTH4/e9/j+HDh3vdtqioCIsWLYrofgcSTBESwAwJFjOA1IKDUBLtwmVxJ415/3YHrOY4/ss6SgI1uQaAhwak4Iktn7U6QbFyVhs6deqEI0eO4MqVKygtLUVBQQF69OiBIUOGtNq2sLAQBQUF7q8dDgfS0tIk3T+xWTIx5zY82C+FGRIkZgCpBavjyaeWFfCdO5hEvc9VfMSG89Hjr8n1qkfvwo6/V/mdVsHKWWXo3LkzYmJiUF1d7fF8dXU1rFarz/cZjUb06tULGRkZmD17Nn7yk5+gqKjI67Ymkwlms9njITWxWfJgvxRmSIiYAaQGvBJKXnmrWrWa4xDfvi3sVxv83uLJ6p4o236Sb/6aXPu7FcrKWeWIjY1FZmYmSktLMXr0aABNhUmlpaWYOXOm6M9xOp0e8z6jqeRoFRbuOO53G2ZJZDADSOk4CKVWfFWtVjtuzi/ydYuHhQPK4q3JtdhboWK3I2kVFBRg0qRJGDRoELKysrB8+XLU1dVhypQpAICJEyeia9eu7iudRUVFGDRoEHr27In6+nq8//77eP3117F69epoHgYAcRXxzJLIYgaQknEQSh4CVcAb0LRuuamNETbHzSsrLBxQjy6d4gJvFMR2JK0xY8bgm2++wfz582Gz2ZCRkYGSkhJ3sdKZM2dgNN6cWVVXV4df/vKXOHv2LNq1a4c+ffpg06ZNGDNmTLQOAYD4ivhkswkLH7qDWSIhZgApBfuEsk+oh7JTlzBu7acBt9v8WDaMBgOLj1So0SlgcPGegJWzy34yABfr6vnz9UGrOSLVcQWTLff16hyxP5daE5sBH80ZCgBs4URBY59QConY2y8Xr9SzybFKBaqcFQB819CI8a8ecD/Pti0UrmCyhaQlpnp+wY/Ssfu4rVVtALOAIonV8eSBt2n0wVflbHz7tgCA2qsNHs+72raUiFzfm6glZouy+KueXz3hLgDAjE2HWxUwMQsoknglVOdarpaReVuCqCbHrFpVv5aVs507mjD7zSNet3XNB1608ziGp1t5O46C0ugU4BQExLdri9rvGrxuw2yRn6/qeQAYXLzHb20As4AigYNQHfPWhinFEoeHBqTglf2VrIDXgeaVs2WnLnkUm7XEti0UCm850xKzJXq8Vc+XnbrEFk4kC96O1ylXqxRvt1pe2V+J6fd393mbhnOBtIltWyjSfOVMS8wWZWEWkFx4JVSHxLRh2vH3Knz46x+i/PS3rIrUCc7Zo0gS05Ipvn1brBp3F+7hqkiKwiwguXAQqkNiV8soP/0tb7XoSFb3RFHzgTNvS0DZqUv8xwn5FShngKYCOKPRwL8/CiM2C7K6J7aqK2AeUDA4CNUh3mohb8S0bXloQAoeWLqXLVsoIOaMerGFE8mFc0J1iLdayBd/bVum398dr+yvZMsWEoU5o25s4URy4JVQHQrmVgvpj7e2LZm3JeCBpXvZsoVEY86oH1s4kdR4JVQnGp0Cyk5dwrtHzuFgZQ3mjeoL4OatFRe2SiHgZtuWhzO6IqfnLSg//a3oli1EQNPfoXmj+vocgALMGTVomQUxRoPougLmAQXCK6E64Ksf6PT7u2PH36s8nrdyPg95wfl9FKySo1VYvKvC62vMGXVjHlCkcBCqca4+fS2vRrj6ga569C4kdIhlZSP5xfl9FAxfueMybxQHoGrGPKBI4e14DQvUDxQAFu86jqzuiR63Wohacs3v8/W3w4Cmq+uc30eB+oMa0JQ7jU5/HURJyZgHFCkchGoY5+1QpLhatgC+5xHPG5WOg5U1ePfIOZSdusRBhk4xd7RPTB645vs2r0dgLlBLvB2vYZy3Q5HkatnScn6x1RKHhwakYPEu9gsk5o5e+MsD1++9r3oE5gK5cBCqYZy3Q5HmrWXLt3X1eGLLZ17nHc/YdJhrgusMc0c/fLVwijEa/NYjMBfIhYNQDWOfPpKCq2UL0DT/j/0CqTnmjr40zwOXQPUIzAVy4ZxQDQtm3g5RKDj/j1pi7hBzgcTiIFRDvE0AD7T0Gm+HUDg4/49aanQKsLSLxc/v64aEDm09XmPu6ANzgcTi7XiNCDQB3Ne8HaJwcP4fNecthxI7xGJ0RiqGp1uZOzrBXCCxZLkSumrVKnTr1g1xcXHIzs7GwYMHfW47ZMgQGAyGVo9Ro0a5t5k8eXKr1/Pz8+U4FEVyTQBvefvDNQG85GiV16XXiMLFfoHk4iuHvq27jtc+/gr2764zd3SCuUBiST4I3bZtGwoKCrBgwQIcPnwYAwYMQF5eHi5cuOB1+7fffhtVVVXux9GjRxETE4P/+I//8NguPz/fY7s33nhD6kNRJDEN6RftZGNokgbn/xHAHCJPzAUSS/JB6EsvvYRp06ZhypQpSE9Px5o1a9C+fXusW7fO6/aJiYmwWq3ux+7du9G+fftWg1CTyeSxXUJCgtSHokicAE7RFmje8fB0K5tVaxxziFoSW4/AZvb6Jumc0OvXr6O8vByFhYXu54xGI3Jzc1FWVibqM1599VWMHTsWHTp08Hh+37596NKlCxISEjB06FA899xzuOWWW7x+Rn19Perr691fOxyOEI5GmTgBnJTA17zj3cdtGFy8h82qNY45RN4EqkdgM3uS9EroxYsX0djYiOTkZI/nk5OTYbPZAr7/4MGDOHr0KB577DGP5/Pz87Fx40aUlpaiuLgYH374IR588EE0NjZ6/ZyioiJYLBb3Iy0tLfSDUhhOACelaDnvePdxW8C5yqQNzCHyxVc9gphaBtI+RbdoevXVV9G/f39kZWV5PD927Fg89NBD6N+/P0aPHo333nsPhw4dwr59+7x+TmFhIex2u/vx9ddfy7D38uAEcFIizhHUF+YQBYP5QC6SDkI7d+6MmJgYVFdXezxfXV0Nq9Xq9711dXXYunUrpk6dGvDP6dGjBzp37oyTJ096fd1kMsFsNns8tIITwEmJOEdQX5hDFAzmA7lIOgiNjY1FZmYmSktL3c85nU6UlpYiJyfH73vfeust1NfXY8KECQH/nLNnz+LSpUtISdHPHJLmk7kt7WKx6tGBbEhPisE5gvozPN2Kp3Jvh6UdG9STf8wHcpG8WX1BQQEmTZqEQYMGISsrC8uXL0ddXR2mTJkCAJg4cSK6du2KoqIij/e9+uqrGD16dKtioytXrmDRokV45JFHYLVacerUKTzzzDPo1asX8vLypD4cRfA1mXveqHQkdIhlQ3qKOs4R1BdvmRTfri2m3NcdM4f2Yg6RB+YDuUg+CB0zZgy++eYbzJ8/HzabDRkZGSgpKXEXK505cwZGo+cF2RMnTuCjjz7CX//611afFxMTg88//xwbNmxAbW0tUlNTMWLECCxevBgmk0nqw4k612TuljNlbPZreGLLYayecBcezugalX0jcnHNEbTZr3md92VA0xUyzhFUP1+ZZP+uAcs/+H/obe3Iq6DkgflALgZBEHQ389fhcMBiscBut6tqfmijU2jV7qY51y/uR3OG8soDRZ1rcALA40Tj+pup9lu0as2RQII5LmYShUrr+aB3YnNE0dXx5ImTuUlN2MRe+5hJFCoxzezZyF77JL8dT5HDydykNmxir23MJAqHv2b2bGSvDxyEqggnc5MauZpVu/ib1zxj02HehlMRZhKFq2U+AMwIPeHteBVhQ2hSOzap1hZmEkUaM0JfOAhVETaEJrXjHEJtYSZRpDEj9IWDUIVrOTF7eLo14GRuIqXiHELtye+XglWPDkRCBzapp/AxI/SFc0IVzN/E7I/mDPU6mZtIyTiHUHtKjlZh8a4K1NQ1uJ9L7BCLeaNYQELBY0boC6+EKpRrYnbL2xKuidm7j9uQ0/MWPJzRFTk9b+EAlFSBcwi1xVdOfVt3HU9sOYySo1VR2jNSK2aEvnAQqkCcmE1axTmE2sGcIikwI/SFg1AF4sRs0jIxTapJ+ZhTJBVmhH5wTqgCcWI2aZ2vJtUAUHbqEuc6qwBziqTkr5G9S6NTYG2EynEQqkCcmE164K2JPVdIUQ/mFEnNWyN7F+aFNvB2vAJxYjbpTaBCPBa4KA9ziqKFeaEdHIQqECdmk56wwEWdmFMUDcwLbeEgVKE4MZv0ggUu6sWcIrkxL7SFc0IVpOUk6+Hp1oATs4nUjgUu6jY83YpOcW1RduoSAAE5PTrjHvYuJokwL7SFg1CF4CRr0isWuKiXt9z638PnmFskGeaFtvB2vAJwkjXpGQtc1Im5RdHAvNAWDkKjjJOsSe9Y4KI+zC2KFuaFtnAQGmWcZE3EAhe1YW5RNDEvtINzQqOMk6yJmohZIUWvVq1ahaVLl8Jms2HAgAFYsWIFsrKyvG67du1abNy4EUePHgUAZGZm4oUXXvC5fSiYWxRtzAtt4CA0yjjJmugmbyuk6H1pvm3btqGgoABr1qxBdnY2li9fjry8PJw4cQJdunRptf2+ffswbtw43HvvvYiLi0NxcTFGjBiBY8eOoWvXrhHZJ+YWKYG/FZUAZocaGARB0N2kHYfDAYvFArvdDrPZHNV9aXQKGFy8Bzb7Na/zqwxousXw0Zyh/OUh3VFy1wi5ciQ7Oxt33303Vq5cCQBwOp1IS0vDrFmzMHfu3IDvb2xsREJCAlauXImJEycG3F7McTG3SOmUnB16IDYfOSc0yjjJmsg7Vl8D169fR3l5OXJzc93PGY1G5ObmoqysTNRnXL16FQ0NDUhM9F4tXF9fD4fD4fEIhLlFSsbsUA8OQhWAk6yJPLH6usnFixfR2NiI5ORkj+eTk5Nhs9lEfcacOXOQmprqMZBtrqioCBaLxf1IS0sT9bnMLVIiZoe6cE6oQnCSNdFNwVRf+5sTpndLlizB1q1bsW/fPsTFeZ+fWVhYiIKCAvfXDocjqIEoc4uUhNmhLhyERomvCdP8pSBi9bVL586dERMTg+rqao/nq6urYbVa/b532bJlWLJkCT744APceeedPrczmUwwmUwR2V+iaGN2qAsHoVHACdNE/rH6uklsbCwyMzNRWlqK0aNHA2gqTCotLcXMmTN9vu/FF1/E888/j7/85S8YNGiQZPvHLCOlYXaoC+eEyowTpokC49J8NxUUFGDt2rXYsGEDKioqMGPGDNTV1WHKlCkAgIkTJ6KwsNC9fXFxMebNm4d169ahW7dusNlssNlsuHLlSkT3i1lGSsTsUBcOQmXECdNE4rD6+qYxY8Zg2bJlmD9/PjIyMnDkyBGUlJS4i5XOnDmDqqqbA77Vq1fj+vXr+MlPfoKUlBT3Y9myZRHbJ2YZKRWzQ13YJ1TGPqFlpy5h3NpPA273xrR7ODeUCMq+3aukfsORJOa4mGWkdErODj0Qm4+cEyojTpgmCg6rr5WJWUZKx+xQBw5CZcQJ00TBa9k1otEpoOzUJZ5YoohZRmrgr+MMl/RUBg5CZeSaMB1oqTtOmCbyjrfYlIFZRmrGHFEOFibJiBOmiULHamzlYJaRWjFHlIWDUJlxqTui4LEaW3mYZaQ2zBHl4e34KOCEaaLgcCk+ZWKWkZowR5SHg1CZeJsEzb/kROKwGpuIwsUcUR4OQmXASdBE4WE1tjIx20hNmCPKwzmhEuMkaKLwcSk+5WG2kdowR5SHg1AJcRI0UWSwGltZmG2kRswR5ZFlELpq1Sp069YNcXFxyM7OxsGDB31uu379ehgMBo9HXJznpXFBEDB//nykpKSgXbt2yM3NxT/+8Q+pDyNowUyCJiL/WI2tHMw2UivmiLJIPid027ZtKCgowJo1a5CdnY3ly5cjLy8PJ06cQJcuXby+x2w248SJE+6vDQbPf5W8+OKLePnll7FhwwZ0794d8+bNQ15eHo4fP95qwBpNnARNFFmsxlYGZhupGXNEOSQfhL700kuYNm0apkyZAgBYs2YNdu3ahXXr1mHu3Lle32MwGGC1Wr2+JggCli9fjt/+9rd4+OGHAQAbN25EcnIy3nnnHYwdO1aaAwkBJ0ETRZ63pfi4BJ+8mG2kdlzSUxkkHYRev34d5eXlKCwsdD9nNBqRm5uLsrIyn++7cuUKbrvtNjidTtx111144YUXcMcddwAAKisrYbPZkJub697eYrEgOzsbZWVlXgeh9fX1qK+vd3/tcDgicXgBcWk7IumxQlt+zDbSKuaJvCSdE3rx4kU0NjYiOTnZ4/nk5GTYbDav7+nduzfWrVuHd999F5s2bYLT6cS9996Ls2fPAoD7fcF8ZlFRESwWi/uRlpYW7qGJwknQRNJihXZ0MNtIi5gn8lNcdXxOTg4mTpyIjIwMPPDAA3j77beRlJSEP/zhDyF/ZmFhIex2u/vx9ddfR3CP/eMkaCJpsEI7uphtpCXMk+iQ9HZ8586dERMTg+rqao/nq6urfc75bKlt27YYOHAgTp48CQDu91VXVyMl5WbIVVdXIyMjw+tnmEwmmEymEI4gMjgJmijyuARf9DHbSCuYJ9Eh6ZXQ2NhYZGZmorS01P2c0+lEaWkpcnJyRH1GY2MjvvjiC/eAs3v37rBarR6f6XA4cODAAdGfGQ2uSdAPZ3RFTs9bGNJEYWKFtjIw20gLmCfRIXl1fEFBASZNmoRBgwYhKysLy5cvR11dnbtafuLEiejatSuKiooAAM8++yzuuece9OrVC7W1tVi6dClOnz6Nxx57DEBT5fxTTz2F5557Drfffru7RVNqaipGjx4t9eGIxuo6ImmxQls5mHekdsyT6JB8EDpmzBh88803mD9/Pmw2GzIyMlBSUuIuLDpz5gyMxpsXZL/99ltMmzYNNpsNCQkJyMzMxCeffIL09HT3Ns888wzq6uowffp01NbWYvDgwSgpKVFMj1BW1xFJjxXaysC8Iy1gnkSHQRAE3c2ydTgcsFgssNvtMJvNEf1sV3Vdy2+q65oAJ+wTRY7r9w2Ax++cHL9vUuZINAVzXMw70pJo5onWiM0RxVXHqxmr64jkxQrt6GHekdYwT+Qn+e14PWF1HZH8WKEdHcw70iLmibw4CI0gVtcRRYe/JfhIGsw70irmiXw4CI0gVtcRRR8rteXBvCM9Yr5EFgehEcTqOqLoYqW2fJh3pDfMl8hjYVIEcT1loujhus/yYt6RnjBfpMFBaISxuo5IfqzUjg7mHekB80U6vB0vAVbXEcmLldrRw7wjrWO+SIeDUImwuo5IPqzUji7mHWkZ80U6HIRGEKvmiKKDldrRx/wjrWK+SIeD0Ahh1RxR9LBSO7qYf6RlzBfpsDApAlg1RxRdrNSOHuYfaR3zRTochIaJVXNEysBKbfkx/0gvmC/S4O34MLFqjkg5WKktL+Yf6QnzJfI4CA0Tq+aIlIWV2vJh/pHeMF8ii4PQMLFqjki5WLEtLeYf0U3Mm+BxEBomVs0RKRMrtqXH/CNqwrwJDQuTwsSqOSLlYcW2PJh/RMybcHAQGgGsmiNSDlZsy4v5R3rGvAkPb8dHCKvmiJSBFdvyY/6RXjFvwsNBaASxao4o+lixHR3MP9Ij5k14OAgNE6vhiJSFFdvRwzwkvWHehIeD0DCwGo5IeVixHR3MQ9Ij5k14WJgUIlbDESkTK7blxzwkvWLehIeD0BCwGo5I2VixLR/mIekd8yZ0vB0fAlbDESkfK7blwTwkYt6EioPQELAajkgdWLEtPeYhURPmTfB4Oz4ErIYjImrCPCSiUPFKaAhYDUekPmwfJA3mIZF/zB7fOAgNgasabsamwzAAHsHLajgi5WH7IOkwD4l8Y/b4x9vxIWI1HJE6sH2Q9JiHRK0xewLjldAwsBqOSNkCtQ8yoKl90PB0q2J/b1etWoWlS5fCZrNhwIABWLFiBbKysrxue+zYMcyfPx/l5eU4ffo0/ud//gdPPfWULPvJPCS6SQvZIwdeCQ2Tqxru4YyuyOl5i67/MhEpTTDtg5Ro27ZtKCgowIIFC3D48GEMGDAAeXl5uHDhgtftr169ih49emDJkiWwWq0y7y3zkMhF7dkjFw5CiUiz1N4+6KWXXsK0adMwZcoUpKenY82aNWjfvj3WrVvndfu7774bS5cuxdixY2EymWTeWyJyUXv2yIW340PEajci5VNz+6Dr16+jvLwchYWF7ueMRiNyc3NRVlYWsT+nvr4e9fX17q8dDkfYn8l8JL1Tc/bIiYPQELDajUgd1Nw+6OLFi2hsbERycrLH88nJyfjyyy8j9ucUFRVh0aJFEfs85iORurNHTrwdHyRWuxGph6t9EHCzXZAL2wc1KSwshN1udz++/vrrkD+L+UjUhNkjDgehQQhU7QY0Vbs1Or1tQUTRoNb2QZ07d0ZMTAyqq6s9nq+uro5o0ZHJZILZbPZ4hIL5SORJrdkjJ96OD0Iw1W5cP5ZIOdTYPig2NhaZmZkoLS3F6NGjAQBOpxOlpaWYOXNmdHfOC+YjUWtqzB45cRAaBFa7EamXq32QmhQUFGDSpEkYNGgQsrKysHz5ctTV1WHKlCkAgIkTJ6Jr164oKioC0FTMdPz4cff/nzt3DkeOHEHHjh3Rq1cvSfeV+UjknRqzRy4chAaB1W5EJKcxY8bgm2++wfz582Gz2ZCRkYGSkhJ3sdKZM2dgNN6cVXX+/HkMHDjQ/fWyZcuwbNkyPPDAA9i3b5+k+8p8JKJgcRAaBFa7Eamf2toHzZw50+ft95YDy27dukEQojPnkvlIFBy1ZZEUZClMWrVqFbp164a4uDhkZ2fj4MGDPrddu3YtfvCDHyAhIQEJCQnIzc1ttf3kyZNhMBg8Hvn5+VIfBqvdiFSu5GgVBhfvwbi1n+JXW49g3NpPMbh4D6u2I4D5SCQes6iJ5IPQYJed27dvH8aNG4e9e/eirKwMaWlpGDFiBM6dO+exXX5+PqqqqtyPN954Q+pDafpzWe1GpEpsHyQ95iNRYMyimwyCxPdusrOzcffdd2PlypUAmqo709LSMGvWLMydOzfg+xsbG5GQkICVK1di4sSJAJquhNbW1uKdd94JaZ8cDgcsFgvsdntY7Uj0fhmdSC0anQIGF+/xWb3tulX80Zyhon+PI5EjSsR8JJKOFFmkRGJzRNI5oZFYdu7q1atoaGhAYqLnPKJ9+/ahS5cuSEhIwNChQ/Hcc8/hllu8V59JsSwdq92I1IPtg+TFfCTyjlnkSdLb8f6WnbPZbKI+Y86cOUhNTUVubq77ufz8fGzcuBGlpaUoLi7Ghx9+iAcffBCNjY1eP6OoqAgWi8X9SEtLC/2giEh12D6IiJSAWeRJ0dXxS5YswdatW7Fv3z7Exd2cYzR27Fj3//fv3x933nknevbsiX379mHYsGGtPqewsBAFBQXurx0OR8gDUd5mIlIftg+KHmYm0U3MIk+SDkLDWXZu2bJlWLJkCT744APceeedfrft0aMHOnfujJMnT3odhJpMJphMpuAPoIWSo1VYtPO4x6X0FEscFvwonRPuiRSM7YOig5lJ5IlZ5EnS2/HNl51zcS07l5OT4/N9L774IhYvXoySkhIMGjQo4J9z9uxZXLp0CSkp0oUaq9mI1Ivtg+THzCRqjVnkSfIWTQUFBVi7di02bNiAiooKzJgxo9Wyc80Ll4qLizFv3jysW7cO3bp1g81mg81mw5UrVwAAV65cwa9//Wt8+umn+Oqrr1BaWoqHH34YvXr1Ql5eniTH0OgUsGjnca//anE9t2jncTQ6o9MkmogCY/sg+TAziXxjFt0k+ZzQYJedW716Na5fv46f/OQnHp+zYMECLFy4EDExMfj888+xYcMG1NbWIjU1FSNGjMDixYsjcsvdG1azEWlDfr8UDE+3co6ixJiZRP4xi5rIUpgUzLJzX331ld/PateuHf7yl79EaM/EYTUbkXawfZD0mJlEgTGLZFq2U+1YzUZEJB4zk4jEUHSLJqVgNRuR9rB1kHSYmUSh01M2cRAqgquabcamwzAAHqGqx2o2IrVj6yBpMTOJQqO3bOLteJFYzUakDWwdJA9mJlFw9JhNvBIaBFazEalboNZBBjS1DhqebuXvdQQwM4nE0Ws2cRAaJFazEakXWwfJj5lJFJhes4m344lIN9g6iIiUSK/ZxCuhIumpWo1Iq9g6KLqYo0Te6TWbOAgVQW/VakRaxdZB0cMcJfJNr9nE2/EB6LFajUirXK2DgJutglzYOkg6zFEi//SaTRyE+hGoWg1oqlZrdHrbgoiUiK2D5MUcJRJHj9nE2/F+6LVajUjr2DpIPsxRIvH0lk0chPqh12o1Ij1g6yB5MEeJgqOnbOLteD/0Wq1GRBQpzFEi8oVXQv3Qa7UakV6wZZD0mKNE4dNqVnEQ6oerWm3GpsMwAB4BquVqNSI9YMsgeTBHicKj5azi7fgA9FitRqR1bBkkL+YoUWi0nlW8EiqC3qrViLQsUMsgA5paBg1Pt/J3PIKYo0TB0UNWcRAqkp6q1Yi0jC2Dooc5SiSeHrKKt+OJSFfYMoiI1EAPWcUroSJotSqNSI/YMkg5mK1EvukhqzgIDUDLVWlEesSWQcrAbCXyTw9Zxdvxfmi9Ko1Ij1wtg4CbLYJc2DJIHsxWosD0kFUchPoQqCoNaKpKa3R624KIlIwtg6KH2UokntazirfjfdBDVRqRnrFlUHQwW4mCo+Ws4iDUBz1UpRHpHVsGyY/ZShQ8rWYVb8f7oIeqNCIiuTFbiciFV0J90ENVGhGxTZDcmK1EkaGF7OIg1AdXVdqMTYdhADzCUitVaUR6xzZB8mO2EoVPK9nF2/F+aL0qjUjP2CYoepitRKHTUnbxSmgAWq5KI9KrQG2CDGhqEzQ83crfdYkwW4mCp7Xs4iBUBK1WpRHpFdsEKQOzlSg4Wssu3o4nIt1hmyAiUiOtZRcHoUSkO2wTRERqpLXs4u14EbTQBoGIbmKbIOVhzhIFprXs4iA0AK20QSCim9gmSFmYs0TiaC27eDveDy21QSAiT2wTpAzMWaLgaCm7eCXUB621QSCi1tgmKLqYs0Sh0Up2cRDqg9baIBCRd2wTFD3MWaLQaSG7eDveB621QSAiUhrmLJG+8UqoD1prg0BE3l2/4cTrZV/hdM1V3JbYHj/L6YbYNvz3uRyYs0SRpbY8k2XPVq1ahW7duiEuLg7Z2dk4ePCg3+3feust9OnTB3Fxcejfvz/ef/99j9cFQcD8+fORkpKCdu3aITc3F//4xz8ius+uNgi+ZlcY0FS9qZY2CETUWtH7x9Fn3p+xeFcFNpadxuJdFegz788oev94tHfNQ6QzVCmYs0SRo5Y8a07yQei2bdtQUFCABQsW4PDhwxgwYADy8vJw4cIFr9t/8sknGDduHKZOnYrPPvsMo0ePxujRo3H06FH3Ni+++CJefvllrFmzBgcOHECHDh2Ql5eHa9cid8vG1QYBQKuAVGMbBCLyVPT+cfxhfyWcLapinALwh/2VigluKTJUKZizRJGhljxrySAIgrfCxIjJzs7G3XffjZUrVwIAnE4n0tLSMGvWLMydO7fV9mPGjEFdXR3ee+8993P33HMPMjIysGbNGgiCgNTUVMyePRv/+Z//CQCw2+1ITk7G+vXrMXbs2ID75HA4YLFYYLfbYTab/W7L/nVE2nP9hhN95v25VWA3ZzQAXy5+0OetrGByJByRztBA5Dqu5pizRKGLRJ5FmtgckXRO6PXr11FeXo7CwkL3c0ajEbm5uSgrK/P6nrKyMhQUFHg8l5eXh3feeQcAUFlZCZvNhtzcXPfrFosF2dnZKCsr8zoIra+vR319vftrh8Mh+hi00gaBiG56vewrv4ENNF1BeL3sK0z9QQ95dsoLKTK0pXDyMVKYs0ShU0ueeSPpIPTixYtobGxEcnKyx/PJycn48ssvvb7HZrN53d5ms7lfdz3na5uWioqKsGjRopCOAdBGGwQiuul0zdWIbicVKTK0pXDzMVKYs0ShUUueeaPckqkIKiwshN1udz++/vrraO8SEUXRbYntI7qdmjEfidRNzXkm6SC0c+fOiImJQXV1tcfz1dXVsFqtXt9jtVr9bu/6bzCfaTKZYDabPR7BaHQKKDt1Ce8eOYeyU5fQGOi6NxEp2s9yuiHQnV6joWm7aJIiQ1sKNx+lxOwlCkwteeaNpIPQ2NhYZGZmorS01P2c0+lEaWkpcnJyvL4nJyfHY3sA2L17t3v77t27w2q1emzjcDhw4MABn58ZjpKjVRhcvAfj1n6KX209gnFrP8Xg4j1cz5hIxWLbGDHtB939bjPtB92j3l9PigxVC2YvkThqyTNvJN+jgoICrF27Fhs2bEBFRQVmzJiBuro6TJkyBQAwceJEj0n3v/rVr1BSUoL//u//xpdffomFCxfib3/7G2bOnAkAMBgMeOqpp/Dcc89hx44d+OKLLzBx4kSkpqZi9OjREd33kqNVmLHpcKtl5Wz2a5ix6TDDkEjFCkem4xf3d291BcFoAH5xf3cUjkyPzo61EOkMVQNmL1Fw1JJnLUm+YtKYMWPwzTffYP78+bDZbMjIyEBJSYl74vyZM2dgNN4cC997773YsmULfvvb3+I3v/kNbr/9drzzzjvo16+fe5tnnnkGdXV1mD59OmprazF48GCUlJQgLi5yq2o0OgUs2nkc3m7+CGjqYbdo53EMT7eygpNIpQpHpmP2iD6KXmFEigxVMmYvUWjUkGctSd4nVInE9K8qO3UJ49Z+GvCz3ph2Dys6iXQoGv005RDt42L2Eqmf2BxR7vA4yi5cFrf6ktjtiIgoMGYvkX5wEOpDl07ibu2L3Y6IiAJj9hLph+RzQtUqq3siUixxsNmveZ2bZABgtTSt6kFE6tXoFLhSj4Iwe4kiR+n5xkGoDzFGAxb8KB0zNh2GAfAIQ9ePb8GP0hX1wySi4HDNcuVh9hJFhhryjbfj/cjvl4LVE+6C1eJ528dqicPqCXcp5odIRMFjGyDlYvYShUct+cYroQHk90vB8HSroi9nE1Fw2AZI+Zi9RKFRU75xECpCjNHAViBEGnKwsqbVFYLmBABV9ms4WFnD3/0oYvYSBU9N+cbb8USkO2wDRERapaZ845VQEZReXUZEwWEbIPViHhP5p6Z84yA0ADVUlxFRcNgGSJ2Yx0SBqSnfeDveD7VUlxFRcFxtgICbbX9c2AZImZjHROKoKd84CPUhUHUZ0FRd1uj0tgURKR3bAKkH85goOGrJN96O90FN1WVEFBq2AVIH5jFR8NSQbxyE+qCm6jIiCh3bACkf85goNErPN96O90FN1WVERFrGPCbSJl4J9UFN1WVEFBq2+1EH5jFR+JSYdxyE+uCqLpux6TAMgEfwKa26jIiCx3Y/6sE8JgqPUvOOt+P9UEt1GREFh+1+1Id5TBQaJecdr4QGoIbqMiISL1C7HwOa2v0MT7fy91xhmMdEwVF63nEQKoLSq8uISDy2+1E35jGReErPO96OJyJdYbsfItILpecdr4SKpMSqMiIKHtv9aA/zmcg7pecdB6EiKLWqjIiCx3Y/2sJ8JvJN6XnH2/EBKLmqjIiC52r3A9xs7+PCdj/qwnwm8k/pecdBqB+BqsqApqqyRqe3LYhIqdjuR/2Yz0TiKDnveDveD6VXlRFR6NjuR92Yz0TiKTXvOAj1Q+lVZUQUHrb7US/mM1FwlJh3vB3vh9KryoiI9Ir5TKR+vBLqh9KryogodGzro27MZ6LwKCEDOQj1w1VVNmPTYRgAj6BTQlUZEYWGbX3Uj/lMFDqlZCBvxweg5KoyIgoe2/poB/OZKHhKykBeCRVBqVVlRBScQG19DGhq6zM83crfb5VgPhOJp7QM5CBUJCVWlRFRcNjWR5uYz0TiKC0DeTueiHSDbX2ISM+UloEchBKRbrCtDxHpmdIykLfjg6SElgZEFBq29dEPZjVRa0rLQA5Cg6CUlgZEFBq29dEHZjWRd0rLQN6OF0lJLQ2IKHRs66NtzGoi/5SUgbwSKoLSWhoQUXjY1kebmNVE4iglAzkIFUFpLQ2IKHxs66M9zGoi8ZSQgbwdL4LSWhoQEVFrzGoideGVUBGU1tKAiELHqmntYlYThSZauSjpldCamhqMHz8eZrMZ8fHxmDp1Kq5cueJ3+1mzZqF3795o164dbr31Vjz55JOw2+0e2xkMhlaPrVu3SnYcrpYGvn4cBjRVXrKtC5GylRytwuDiPRi39lP8ausRjFv7KQYX72GxikYwq4mCF81clHQQOn78eBw7dgy7d+/Ge++9h/3792P69Ok+tz9//jzOnz+PZcuW4ejRo1i/fj1KSkowderUVtu+9tprqKqqcj9Gjx4t2XG4WhoAaBVubOtCpA6smtY+ZjVRcKKdiwZBELwVEoatoqIC6enpOHToEAYNGgQAKCkpwciRI3H27FmkpqaK+py33noLEyZMQF1dHdq0aZo9YDAYsH379pAHng6HAxaLBXa7HWazWfT72HuOSJ0anQIGF+/xWbTiatD80ZyhogcooeaI0mnhuJjVRIFJkYsuYnNEsjmhZWVliI+Pdw9AASA3NxdGoxEHDhzAj3/8Y1Gf4zoA1wDU5YknnsBjjz2GHj164PHHH8eUKVNgMHj/JtXX16O+vt79tcPhCOGIlNPSgIiCw6ppfWFWEwWmhFyUbBBqs9nQpUsXzz+sTRskJibCZrOJ+oyLFy9i8eLFrW7hP/vssxg6dCjat2+Pv/71r/jlL3+JK1eu4Mknn/T6OUVFRVi0aFFoB9KCEloaEFFwWDWtP8xqIv+UkItBD0Lnzp2L4uJiv9tUVFSEvEMuDocDo0aNQnp6OhYuXOjx2rx589z/P3DgQNTV1WHp0qU+B6GFhYUoKCjw+Oy0tLSw9o8VtkTqwappcmF2EzVRQi4GPQidPXs2Jk+e7HebHj16wGq14sKFCx7P37hxAzU1NbBarX7ff/nyZeTn56NTp07Yvn072rZt63f77OxsLF68GPX19TCZTK1eN5lMXp8PFecbEamLq2raZr/mdTUd19wnVk1rG7Ob6CYl5GLQ1fFJSUno06eP30dsbCxycnJQW1uL8vJy93v37NkDp9OJ7Oxsn5/vcDgwYsQIxMbGYseOHYiLCzwCP3LkCBISEiI60PQl2pVkRBQ8Vk0Ts5vIkxJyUbIWTX379kV+fj6mTZuGgwcP4uOPP8bMmTMxduxYd2X8uXPn0KdPHxw8eBDAzQFoXV0dXn31VTgcDthsNthsNjQ2NgIAdu7ciT/+8Y84evQoTp48idWrV+OFF17ArFmzpDoUt0DrEgNN6xI3OiVpOEBEYcjvl4LVE+6C1eL5D1urJQ6rJ9zFK2Eaxuwm8i7auSjpikmbN2/GzJkzMWzYMBiNRjzyyCN4+eWX3a83NDTgxIkTuHr1KgDg8OHDOHDgAACgV69eHp9VWVmJbt26oW3btli1ahWefvppCIKAXr164aWXXsK0adOkPBQAyqgkI6LQsWpan5jdRL5FMxclHYQmJiZiy5YtPl/v1q0bmrcpHTJkCAK1Lc3Pz0d+fn7E9jEYSqgkI6LwsGpaf5jdRP5FKxclXTFJa5RQSUZE+hDssscA8Morr2DIkCEwm80wGAyora2VZ2cVjtlNpEwchAaB6xITqVejU0DZqUt498g5lJ26pPj5f8EuewwAV69eRX5+Pn7zm9/ItJfqwOwmEkfunJT0drzWuCrJZmw6DAPgMcmdFbZEyqW21jwVFRUoKSnxWPZ4xYoVGDlyJJYtW+Zz2eOnnnoKALBv3z7Rf1akVpRTMmY3UWDRyEleCQ1StCvJiCg4amzNE2jZ40gqKiqCxWJxP8JdyEOpmN1EvkUrJ3klNASssCVSh0CteQxoas0zPN2qqN/fSCx7LJYUK8opFbObqLVo5iQHoSFihS2R8imtNY9cyx4HI9Iryikds5vIUzRzkoNQItIspbXmkWPZYyKiYEQzJzkIDVOjU+CtHSKFUlprnqSkJCQlJQXcrvmyx5mZmQDELXtMoWOWk15FMyc5CA2D2ipuifTG1ZrHZr/mdb6TAU2FKUprzdN82eM1a9agoaHB67LHw4YNw8aNG5GVlQUA7mWOT548CQD44osv0KlTJ9x6661ITFTWMSoJs5z0LJo5yer4EKmx4pZIb1yteQC06hGp9NY8mzdvRp8+fTBs2DCMHDkSgwcPxiuvvOJ+veWyxwCwZs0aDBw40L2M8f3334+BAwdix44dsu+/WjDLSe+imZMGIdA6mRrkcDhgsVhgt9thNpuDfn+jU8Dg4j0+J/K6/tXw0Zyhijy5EemNFFe6ws0RpdLqcXnDLCe6KZI5KTZHeDs+BEqruCUi/9iah7xhlhPdFI2c5CA0BEqruCWiwNiah1pilhN5kjsnOQgNgdIqbonIO1Y8kz/MciLv5MpODkJDoNaKWyI9YcUzBcIsJ2pNzuxkdXwI1FxxS6QHrHgmMZjlRJ7kzk4OQkOU3y8FqyfcBavF8zaN1RKH1RPu4pUWoigJtA4y0LQOcqNTd41ByAtmOVGTaGQnb8eHgRW3RMrDimcKFrOcKDrZyUFomFhxS6QsrHimUDDLSe+ikZ0chEYQK3GJoo8VzxRpzHbSg2hkJwehEcJKXCJlYMUzRRKznfQiGtnJwqQIYCUukXKw4pkihdlOehKN7OQgNEysxCVSHlY8U7iY7aRHcmcnb8eHiZW4RMrEimcKB7Od9ErO7OQgNEysxCVSLlY8U6iY7aRncmUnb8eHiZW4RETaw2wnkh6vhIaJlbhEysJ2OhQJzHbSOzmylIPQMLmqyWZsOgwD4BFWrMQlkhfb6VCkMNtJz+TKUt6OjwBW4hJFH9vpUKQx20mP5MxSXgmNEFbiEkVPoHY6BjS10xmebuXvJAWF2U56IneWchAaQazEJYoOttMhKTHbSS/kzlLejici1WM7HSKi8MmdpbwSKhFW6BLJh+10KFqY9aQlcmcpB6ESYIUukbzYToeigVlPWiN3lvJ2fISxQpdIfq52OsDN9jkubKdDUmDWkxbJnaUchEZQoKoyoKmqrNHpbQsiCgfb6ZBcmPWkZXJmKW/HRxArdImii+10SA7MetI6ubKUg9AIYoUuUXSwOITkxKwnPYgxGpDVPdGdrQcrayKerRyERhArdInkx+IQkhuznvRAjmzlnNAIclWV+fo3ggFNP0BW6BJFBotDKBqY9aR1cmUrB6ERxApdIvmwOISihVlPWiZntnIQGmGs0CWSRzDFIUSRxqwnrZIzWyWdE1pTU4NZs2Zh586dMBqNeOSRR/C73/0OHTt29PmeIUOG4MMPP/R47he/+AXWrFnj/vrMmTOYMWMG9u7di44dO2LSpEkoKipCmzbKmOLKCl0i6bE4hKKNWU9aJGe2SjpqGz9+PKqqqrB79240NDRgypQpmD59OrZs2eL3fdOmTcOzzz7r/rp9+/bu/29sbMSoUaNgtVrxySefoKqqChMnTkTbtm3xwgsvSHYswYoxGjxaczQ6BZSdusSgIooQFoeQErTM+ubYtYHUSM5slWwQWlFRgZKSEhw6dAiDBg0CAKxYsQIjR47EsmXLkJqa6vO97du3h9Vq9fraX//6Vxw/fhwffPABkpOTkZGRgcWLF2POnDlYuHAhYmNjW72nvr4e9fX17q8dDkeYRxccVu8SRR6X6iQlY+6TWsmZrZLNCS0rK0N8fLx7AAoAubm5MBqNOHDggN/3bt68GZ07d0a/fv1QWFiIq1evenxu//79kZyc7H4uLy8PDocDx44d8/p5RUVFsFgs7kdaWlqYRyceq3eJpMHiEFIq5j6pmZzZKtkg1GazoUuXLh7PtWnTBomJibDZbD7f9+ijj2LTpk3Yu3cvCgsL8frrr2PChAken9t8AArA/bWvzy0sLITdbnc/vv7661APKyis3iWSFotDSGmY+6QFcmVr0Lfj586di+LiYr/bVFRUhLxD06dPd/9///79kZKSgmHDhuHUqVPo2bNnSJ9pMplgMplC3qdQcWk3IumxOISUhLlPWiFHtgY9CJ09ezYmT57sd5sePXrAarXiwoULHs/fuHEDNTU1Pud7epOdnQ0AOHnyJHr27Amr1YqDBw96bFNdXQ0AQX2uHFi9SyQPf8UhRHJi7pOWSJ2tQQ9Ck5KSkJSUFHC7nJwc1NbWory8HJmZmQCAPXv2wOl0ugeWYhw5cgQAkJKS4v7c559/HhcuXHDf7t+9ezfMZjPS09ODPBppsXqXKPJYcUxKxtwnLZE6byWrju/bty/y8/Mxbdo0rFmzBg0NDZg5cybGjh3rrow/d+4chg0bho0bNyIrKwunTp3Cli1bMHLkSNxyyy34/PPP8fTTT+P+++/HnXfeCQAYMWIE0tPT8bOf/QwvvvgibDYbfvvb3+KJJ56Iyi13f1i9SxRZrDgmpWPuk1aofu34zZs3o0+fPhg2bBhGjhyJwYMH45VXXnG/3tDQgBMnTrir32NjY/HBBx9gxIgR6NOnD2bPno1HHnkEO3fudL8nJiYG7733HmJiYpCTk4MJEyZg4sSJHn1FlYLVu0SRw4pjUgPmPmmBXHlrEARBdyV6DocDFosFdrsdZrNZ8j+PV2+IwtPoFDC4eI/Pgg/X1aWP5gyV7eQud47IRavHJTfmPqlVJPJWbI4oY51LjWP1LlF4WHFMasPcJ7WSM285CJUJl/EkCh0rjkmNuKQnqZFm1o4n73ibhig4rDgmLeE5gJRMzryVtDCJWmNxBVHwXBXHvq4TGdB0EmfFMSkdzwGkdHLmLQehMuJybkShYcUxaQHPAaQGmlg7nloLZrIvEXniOvGkdjwHkFoodu14Ch2LK4jCw4pjUjOeA0hNFLl2PIWOxRVEwfNWRcw2TKRGPAcQeeIgVEZczo0oOKwiJi3hOYDURPXLdpInFlcQiccqYtIangNILeTKXw5CZcbiCqLAWEXcpKamBuPHj4fZbEZ8fDymTp2KK1eu+N1+1qxZ6N27N9q1a4dbb70VTz75JOx2u4x7Tf7wHEBKJ2f+8nZ8FLC4gsg/LtPZZPz48aiqqsLu3bvR0NCAKVOmYPr06diyZYvX7c+fP4/z589j2bJlSE9Px+nTp/H444/j/Pnz+NOf/iTz3pMvPAeQknHZTh3wtpwbl3EjasIqYqCiogIlJSU4dOgQBg0aBABYsWIFRo4ciWXLliE1NbXVe/r164f//d//dX/ds2dPPP/885gwYQJu3LiBNm1aR359fT3q6+vdXzscDgmOhlrikp6kVFy2U4dYgEF0E6uIgbKyMsTHx7sHoACQm5sLo9GIAwcO4Mc//rGoz7Hb7TCbzV4HoABQVFSERYsWRWSfKXw8F1C0cdlOnWEBBpEnLtMJ2Gw2dOnSxeO5Nm3aIDExETabTdRnXLx4EYsXL8b06dN9blNYWAi73e5+fP3112HtN4WO5wJSAi7bqSMswCBqTctVxHPnzoXBYPD7+PLLL8P+cxwOB0aNGoX09HQsXLjQ53Ymkwlms9njQfLjuYCUgst26giXcSPyTqtVxLNnz0ZFRYXfR48ePWC1WnHhwgWP9964cQM1NTWwWq1+/4zLly8jPz8fnTp1wvbt29G2bVspD4kigOcCUhIu26kTLMAguqllQcbwdKvmqoiTkpKQlJQUcLucnBzU1taivLwcmZmZAIA9e/bA6XQiOzvb5/scDgfy8vJgMpmwY8cOxMVpd96slvBcQEqT3y8FQ/sk4/Wyr3C65ipuS2yPn+V0Q2ybyF2/5CA0yliAQdSEBRme+vbti/z8fEybNg1r1qxBQ0MDZs6cibFjx7or48+dO4dhw4Zh48aNyMrKgsPhwIgRI3D16lVs2rQJDofDXe2elJSEmJiYaB4S+cFzASmNt0z+40eVXDFJS1iAQcSCDF82b96MPn36YNiwYRg5ciQGDx6MV155xf16Q0MDTpw4gatXrwIADh8+jAMHDuCLL75Ar169kJKS4n6w4EjZeC4gJeGKSTqh5QIMIjFYkOFbYmIitmzZgsuXL8Nut2PdunXo2LGj+/Vu3bpBEAQMGTIEADBkyBAIguD10a1bt+gcBInCcwEphZyZzEGoAmi1AINIDBZkEDXhuYCUQM5M5pxQhfC1jBsAlJ26pJmiDKKWWJBBdJOYJT25ohJJiSsm6VTLZdxYqEF6wIIMIk/+lvTkeYGkxhWTiIUapBssyCASh+cFkgNXTNI5FmqQnrAggygwnhdILlwxSedYqEF6w4IMIv94XiA5ccUkHWOhBumBHlZHIooUnhdIbsPTrehkaouyf14E0DRP+Z4et0Q0kzkIVSAWapDWsbiCKDg8L5CcvGX0/x4+G/GM5u14BWKhBmkZiyuIgsfzAslFzozmIFSBWKhBWsXiCqLQ8LxAcpA7ozkIVSgWapAWsbiCKHQ8L5DU5M5ozglVMH8rZ3DFDFIjFlcQhYcrKpGU5M5oDkIVztvKGSzqILVicQVR+LiiEklF7ozm7XiVYVEHqRmLK4ikw/MDhUvujOYgVEVY1EFqx+IKImnw/ECRIHdGcxCqIizqIC1gcQVR5PH8QJEiZ0ZzTqiKsKiD1MhbkYSY4goiEo/nB4qURqcAS7tYPJPXGzV115HY0QSrWZqM5iBURVjUQWoTqEjCV3EFEQWH5weKBH+ZLcVFAt6OVxEWdZCasEiCSD48P1C4opHZHISqCIs6SC1YJEEkL54fKBzRymxJB6E1NTUYP348zGYz4uPjMXXqVFy5csXn9l999RUMBoPXx1tvveXeztvrW7dulfJQFMPfhOFVjw6EpV0s3j1yDmWnLvEET1HDIgki+YkpKGl0Cig7dYnnCfIQrcyWdE7o+PHjUVVVhd27d6OhoQFTpkzB9OnTsWXLFq/bp6WloarK83LvK6+8gqVLl+LBBx/0eP61115Dfn6+++v4+PiI779SeSvq+LbuOhbvYoNiUgYWSRBFh7+iPzayJ1+ildmSDUIrKipQUlKCQ4cOYdCgQQCAFStWYOTIkVi2bBlSU1NbvScmJgZWq9Xjue3bt+OnP/0pOnbs6PF8fHx8q231pPmKGSVHq/DElsOtLqO75nGw7Q3JjUUSRNHja6W9GZt4niDvopXZkt2OLysrQ3x8vHsACgC5ubkwGo04cOCAqM8oLy/HkSNHMHXq1FavPfHEE+jcuTOysrKwbt06CILvWwr19fVwOBweD63g3DtSIhZJECkHzxMUSLQyW7JBqM1mQ5cuXTyea9OmDRITE2Gz2UR9xquvvoq+ffvi3nvv9Xj+2WefxZtvvondu3fjkUcewS9/+UusWLHC5+cUFRXBYrG4H2lpacEfkEJx7h0pRfO5ZgcrazBvFIskiJSA5wkKJMZowLxR6V7/oSJlZgd9O37u3LkoLi72u01FRUXIO+Ty3XffYcuWLZg3b16r15o/N3DgQNTV1WHp0qV48sknvX5WYWEhCgoK3F87HA7NDEQ5946UwNdcs+n3d8eOv1d5PG/lHDQiWfE8QYGUHK3C4l3Hvb4mZWYHPQidPXs2Jk+e7HebHj16wGq14sKFCx7P37hxAzU1NaLmcv7pT3/C1atXMXHixIDbZmdnY/Hixaivr4fJZGr1uslk8vq8FnDuHUWbv7lmr+yvxKpHByKhg4krIxFFCc8T5I+vDHeZN6qvZBcNgh6EJiUlISkpKeB2OTk5qK2tRXl5OTIzMwEAe/bsgdPpRHZ2dsD3v/rqq3jooYdE/VlHjhxBQkKCZgea/rjmcdjs13xeRrdy7h1JJNBcMwOAxbsq8NGcoRx4EkUJzxPki78MB25meF6/FHWtmNS3b1/k5+dj2rRpOHjwID7++GPMnDkTY8eOdVfGnzt3Dn369MHBgwc93nvy5Ens378fjz32WKvP3blzJ/74xz/i6NGjOHnyJFavXo0XXngBs2bNkupQFE1sg2IA7A1HEce5ZkTKF0wje/YR1ZdoZ7ikfUI3b96MmTNnYtiwYTAajXjkkUfw8ssvu19vaGjAiRMncPXqVY/3rVu3Dt/73vcwYsSIVp/Ztm1brFq1Ck8//TQEQUCvXr3w0ksvYdq0aVIeiqK5GhS3nJPnmscBAIOL97A3HEUc55oRqUOg80R+vxT2EdWhaGe4QfDX20ijHA4HLBYL7HY7zGZztHcnYhqdQqsGxbuP27zO9XD965e94SgcZacuYdzaTwNu98a0e1r1LVQ7reaIVo+Lmng7T7ga2fNcoT9SZbjYHJH0SijJq2WDYjHz9RbtPI7h6VbO16OQcK4Zkbp4a2TPc4V+RTvDJV07nqIr2nM9SJvYD5RIW3iu0Lexd98qe39QF14J1bBoz/Ug7WE/UCLt4blCn7zleXNyZDgHoRrG3nAUSewHSqRNPFfoT6DeoE/n3o6ZQ2+XPMN5O17DuH43RYqYtacX76pAVvdEPJzRFTk9b+EAlEgleK7QFzG9Qbce+lqWfeEgVMPE9IabN6ovDlbWsCcc+cU5Y0TaxT6i+qKkPOfteI3z1xvuoQEpWLyrgj3hKCDOGSPSNvYR1Q8l5TkHoTqQ3y8Fw9OtHr3hvq27jie2eJ/fN2PTYfaEIw+cM0akfd7OFYH6iPKcoT5KynMOQnWieW+4RqeAwcV72BOORIt2Lzkikgf7iGqfkvKcc0J1SEnzQUiZWs77AiB6zhgRaQvPGdoTzd6gzfFKqA4paT4IKY+/eV+B5owRkfbwnKEdSugN2hwHoTqkpPkgpCxi5n19NGeo1zljRKRNPGdog1J6gzbH2/E6JLYnXOZtCWzFoSNieoEu2nkcAJDT8xb2AyXSiWD6iLKFkzIpqTdoc7wSqkOunnAzNh2GAfD4S+kKmYcGpOCBpXvZikNHgpn31bJwgYi0S8w5Y8GP0rH7uI0tnBRKqfnOK6E65eoJZ7V43j6x/msd8Ff2V7b6C+u6JVtytErOXSWZcN4XEfni75yxesJdAIAZmw7zvKFQSs13XgnVMW894TJvS8ADS/eyFYdONDoF98//4uV6Ue/hvC8iffLVRxQA2/4pXOeOJlHbyZ3vHITqXMuecGWnLinykj1FnrcqSaMB8DWFi71AichbH1GeN5St5GgVFu445nebaOU7B6HkQamX7CmyfFVJ+huAAuwFSkSt8byhXIEq4oHo5jsHoeRB7KX4zh1MKDt1iW16VChQlSTQ+oooe4ESkS/BtHBqPgWI5w5picl6ILr5zkEoeRCznJelfVvMfuvvsDlYAalGgaokgaYB6LxRfdG5k4knCiLyS+wykN/WXcfg4j2snpeJmKwHgGU/GYD7bu8swx61xup48uBqxQF4X55RAFB7tcFjAAqwAlJNxN4S69zJxF6gRBRQoPMG0NT274ktrJ6Xk9isv1gnrihVChyEUiu+WnEkm02Ib9/W63uaNzNnc2Ll8NY4mqufqEdNTQ3Gjx8Ps9mM+Ph4TJ06FVeuXPH7nl/84hfo2bMn2rVrh6SkJDz88MP48ssvZdpj0it/LZxWPToQO/5eFXAhDJ47IqfRKaii4wlvx5NX3lpxOAUB4/94wOd7WAGpLL7WgJ83qq+oW2esgo++8ePHo6qqCrt370ZDQwOmTJmC6dOnY8uWLT7fk5mZifHjx+PWW29FTU0NFi5ciBEjRqCyshIxMTEy7j3pja8WTkptlK5VgdaHd1FC1nMQSj61bMXx7pFzot7HCsjo87cG/BNbPnMvSOBv9RPego+uiooKlJSU4NChQxg0aBAAYMWKFRg5ciSWLVuG1NRUr++bPn26+/+7deuG5557DgMGDMBXX32Fnj17ttq+vr4e9fU3r5g4HI4IHwnpibcWTqyel4+YanhAOVnP2/EkWrCV81w7ODrErAG/4+9VWPWo79VPWCQQfWVlZYiPj3cPQAEgNzcXRqMRBw74viPRXF1dHV577TV0794daWlpXrcpKiqCxWJxP3xtRxSqYKvnef4IjdhqeEA5Wc8roSQaK+fVQeytr4QOsfhozlC2S1Eom82GLl26eDzXpk0bJCYmwmaz+X3v73//ezzzzDOoq6tD7969sXv3bsTGxnrdtrCwEAUFBe6vHQ4HB6IUUayel4fYavh5o/pi8n3dFZH1vBJKorFyXplaXjlo+f335cLla+5bZ6yCl8/cuXNhMBj8PsItJBo/fjw+++wzfPjhh/j+97+Pn/70p7h2zfvfC5PJBLPZ7PEgiiRWz0vLdQ74s8jvUedOJsVkPa+EUlBcFZAtJz0nm024dsOJ2qsNrd7DtYOl420CemIH7x0MWmL1e3TMnj0bkydP9rtNjx49YLVaceHCBY/nb9y4gZqaGlitVr/vd91av/3223HPPfcgISEB27dvx7hx48LdfaKQ+Dp3WP9VLLl4VwXXng+B2CKk5pSU/RyEUtBYOa8Mviag19S1/odAc0qoiNSzpKQkJCUlBdwuJycHtbW1KC8vR2ZmJgBgz549cDqdyM7OFv3nCYIAQRA8io+IooHV85EltgjJRYnZz0EohSTUynmb/Tsu9xkBYiegs/pdvfr27Yv8/HxMmzYNa9asQUNDA2bOnImxY8e6K+PPnTuHYcOGYePGjcjKysI///lPbNu2DSNGjEBSUhLOnj2LJUuWoF27dhg5cmSUj4govOp5nj9uCqYICVBu9nMQShEh9vL+4l0VqKm77v6ak85DI3YCekKHWI/vN9eAV5fNmzdj5syZGDZsGIxGIx555BG8/PLL7tcbGhpw4sQJXL16FQAQFxeH//u//8Py5cvx7bffIjk5Gffffz8++eSTVkVORErB80fwxJ4DXJSa/RyEUkQEqn50aR4gwM1J50poFaFkjU7B4xaW2OKjeaP6wmppxysHKpWYmOi3MX23bt0gCDd/41JTU/H+++/LsWtEEcPzR3AanQI+PnlR1LYTc27Dg/1SFJv9HIRSRLiqH2dsOtzqFrA/nHQeWDjFR1ZLO86hIiJF4/lDvGALkR7sl6LocwBbNFHE+Fo7ONCAqfmkczYq9uSaeN4ycMQUH6UobAI6EZEvPH8E5ut84I1azgG8EkoR5a360ea4hqe3HQn43t3HbSh484guGxW3vN3uCg4WHxGRXvD80Zrr3GCzf+ezjVVLajoHcBBKEdey+rHs1CVR71v38VetntPDnB9vt1dSLHEYe3cai4+ISFd4/rgplB6ggLrOARyEkuTETDo3GgBvd06az/kZ2icZ5ae/1VSRja8+bzb7NfzPB/8Q9RksPiIirdLr+SPYHqAuM3/YE08P763oY2uOg1CSnL9J566v/U3dcc35uaeoVLXtOYK93R5M8LD4iIi0Si/nj+bniM4dTFi4Q3wP0Obu65WkmgEowEEoycTfkm0j+1nxqpdbKS0Fas/hbaCnhF/GcG+3+6LE1S+IiCJN6+ePUG+7N6fW8wEHoSQbf0u2iQmRlprfanE6gcW7Wg/05PyXrrcQ233cFvbtdoDFR0Skb1o5f7Q8T3xbdx1PbAn+tntzaj4fGITmnY51wuFwwGKxwG63w2w2R3t3dK/RKWBw8Z6AjYqD5fpVFPsv3XBe9/YvWavZhGs3nKi96r+dUiBP534fWw+d0VzVp9ppNUe0elykTXKcP7wNfmOMBr/nBF8XJVqeJ3zNZw2GEs8HYnOEg1CGrCK4JmEDwc2HDMR1i2LeqHS//9L1dctczOsAQppALnbfP5ozFAAUOdVAz7SaI1o9LtIuKc8flvZtEdcmxmOVuhRLHB4akIIdf6/yeU5oeb6Ib9827AsSzSV2aIt5/3YHrGZlng+iPgh9/vnnsWvXLhw5cgSxsbGora0N+B5BELBgwQKsXbsWtbW1uO+++7B69Wrcfvvt7m1qamowa9Ys7Ny5072W8u9+9zt07NhR9L4xZJXJ18pAgRqzh8r1Kzv9/u54ZX9lq/AS87qAyISLr9vtam0togdazRGtHhdpm9znD2+CWe0pnD8DUP65IeqD0AULFiA+Ph5nz57Fq6++KmoQWlxcjKKiImzYsAHdu3fHvHnz8MUXX+D48eOIi2taReHBBx9EVVUV/vCHP6ChoQFTpkzB3Xff7Xd95ZYYssrV8hZG5m0JeGDp3ojfamku0O2QSNwu8Ye329VJqzmi1eMi7YvG+UNuajk3RH0Q6rJ+/Xo89dRTAQehgiAgNTUVs2fPxn/+538CAOx2O5KTk7F+/XqMHTsWFRUVSE9Px6FDhzBo0CAAQElJCUaOHImzZ88iNTXV62fX19ejvr7e/bXD4UBaWhpDViV83WqR41+dUuLtdnXT6mBNq8dF+qTm84cBQLLZhP/+aQYuXqlX1blBbI4oZu34yspK2Gw25Obmup+zWCzIzs5GWVkZAKCsrAzx8fHuASgA5Obmwmg04sCBAz4/u6ioCBaLxf1IS0uT7kAo4nytKWy1xOH3jw5EiiUOyv+V9NSymtG1SsjDGV2R0/MWVYQMEZHSqfX84dqnhQ/dgft6ddbsuUExLZpsNhsAIDk52eP55ORk92s2mw1dunTxeL1NmzZITEx0b+NNYWEhCgoK3F+7roSSevhqzxFjNMBoNPhtZBxNvia2q2lZNSIiNQvl/CG3ltO+9HKOCGoQOnfuXBQXF/vdpqKiAn369AlrpyLNZDLBZDJFezcoTC3XFHbx18h43qi+WLyrIuCSb4LgO4D8ve4aZNr/VZjkrbhoyb/39xmAREQkvaDPH/9qsWe/2iDp4NR1Flg5biASOph0d44IahA6e/ZsTJ482e82PXr0CGlHrFYrAKC6uhopKTdH/tXV1cjIyHBvc+HCBY/33bhxAzU1Ne73kz6FeqUUAKb9oKn6PdTXl/x7fwCtW3K0/Jcsl9YkIlIeX+cP12IjYq+SNt/O1525lt1U9HLF05egBqFJSUlISkqSZEe6d+8Oq9WK0tJS96DT4XDgwIEDmDFjBgAgJycHtbW1KC8vR2ZmJgBgz549cDqdyM7OlmS/SD1CuVLq+uUfeGtCWK8D4NVOIiKV8nb+8HXu8NUn1OqnT6jrNZ4nPElWHX/mzBnU1NRgx44dWLp0Kf7v//4PANCrVy93T88+ffqgqKgIP/7xjwE0tWhasmSJR4umzz//vFWLpurqaqxZs8bdomnQoEFs0UQBSbliEumPVnNEq8dFFCpf2R/sikl6Ol+IzRHJCpPmz5+PDRs2uL8eOHAgAGDv3r0YMmQIAODEiROw2+3ubZ555hnU1dVh+vTpqK2txeDBg1FSUuIegALA5s2bMXPmTAwbNszdrP7ll1+W6jBIQ3xdKY3U60REpD2+st/fOYHnC3G4bCf/pU9EIdBqjmj1uIhIPqrrE0pERERE+sFBKBERERHJjoNQIiIiIpIdB6FEREREJDsOQomIiIhIdhyEEhEREZHsOAglIiIiItlxEEpEREREsuMglIiIiIhkx0EoEREREcmOg1AiIiIikh0HoUREREQkuzbR3oFoEAQBAOBwOKK8J0SkVq78cOWJVjAfiShcYvNRl4PQy5cvAwDS0tKivCdEpHaXL1+GxWKJ9m5EDPORiCIlUD4aBK39M14Ep9OJ8+fPo1OnTjAYDAG3dzgcSEtLw9dffw2z2SzDHiqHXo+dx83jDkQQBFy+fBmpqakwGrUzsynYfFQLvf7dlgK/l5Gj1e+l2HzU5ZVQo9GI733ve0G/z2w2a+ovSTD0euw8bn0J9ri1dAXUJdR8VAu9/t2WAr+XkaPF76WYfNTOP9+JiIiISDU4CCUiIiIi2XEQKoLJZMKCBQtgMpmivSuy0+ux87h53KQt/BlHDr+XkaP376UuC5OIiIiIKLp4JZSIiIiIZMdBKBERERHJjoNQIiIiIpIdB6FEREREJDsOQomIiIhIdhyE+vD888/j3nvvRfv27REfHy/qPYIgYP78+UhJSUG7du2Qm5uLf/zjH9LuaITV1NRg/PjxMJvNiI+Px9SpU3HlyhW/7xkyZAgMBoPH4/HHH5dpj0O3atUqdOvWDXFxccjOzsbBgwf9bv/WW2+hT58+iIuLQ//+/fH+++/LtKeRFcxxr1+/vtXPNi4uTsa9jYz9+/fjRz/6EVJTU2EwGPDOO+8EfM++fftw1113wWQyoVevXli/fr3k+0mRFWye1dTUYNasWejduzfatWuHW2+9FU8++STsdruMe60Mes1HKQTzvVy7di1+8IMfICEhAQkJCcjNzQ34vVczDkJ9uH79Ov7jP/4DM2bMEP2eF198ES+//DLWrFmDAwcOoEOHDsjLy8O1a9ck3NPIGj9+PI4dO4bdu3fjvffew/79+zF9+vSA75s2bRqqqqrcjxdffFGGvQ3dtm3bUFBQgAULFuDw4cMYMGAA8vLycOHCBa/bf/LJJxg3bhymTp2Kzz77DKNHj8bo0aNx9OhRmfc8PMEeN9C0nFzzn+3p06dl3OPIqKurw4ABA7Bq1SpR21dWVmLUqFH44Q9/iCNHjuCpp57CY489hr/85S8S7ylFUrB5dv78eZw/fx7Lli3D0aNHsX79epSUlGDq1Kky7nX06TUfpRDs93Lfvn0YN24c9u7di7KyMqSlpWHEiBE4d+6czHsuE4H8eu211wSLxRJwO6fTKVitVmHp0qXu52prawWTySS88cYbEu5h5Bw/flwAIBw6dMj93J///GfBYDAI586d8/m+Bx54QPjVr34lwx5GTlZWlvDEE0+4v25sbBRSU1OFoqIir9v/9Kc/FUaNGuXxXHZ2tvCLX/xC0v2MtGCPW+zffzUBIGzfvt3vNs8884xwxx13eDw3ZswYIS8vT8I9o0gKNc9aevPNN4XY2FihoaFBit1UJL3moxSC/V62dOPGDaFTp07Chg0bpNrFqOKV0AiprKyEzWZDbm6u+zmLxYLs7GyUlZVFcc/EKysrQ3x8PAYNGuR+Ljc3F0ajEQcOHPD73s2bN6Nz587o168fCgsLcfXqVal3N2TXr19HeXm5x8/KaDQiNzfX58+qrKzMY3sAyMvLU83PFgjtuAHgypUruO2225CWloaHH34Yx44dk2N3o0oLP2+9CyfPmrPb7TCbzWjTpo0Uu6k4es1HKYSauc1dvXoVDQ0NSExMlGo3o4qD0Aix2WwAgOTkZI/nk5OT3a8pnc1mQ5cuXTyea9OmDRITE/0ew6OPPopNmzZh7969KCwsxOuvv44JEyZIvbshu3jxIhobG4P6WdlsNlX/bIHQjrt3795Yt24d3n33XWzatAlOpxP33nsvzp49K8cuR42vn7fD4cB3330Xpb2iYISaZ81dvHgRixcvFjUlSSv0mo9SCOV72dKcOXOQmpraapCvFboahM6dO7dVkUXLx5dffhnt3Yw4qY97+vTpyMvLQ//+/TF+/Hhs3LgR27dvx6lTpyJ4FBQNOTk5mDhxIjIyMvDAAw/g7bffRlJSEv7whz9Ee9dIp+TKcYfDgVGjRiE9PR0LFy4Mf8eJgrRkyRJs3boV27dvV2VBqBj6uL/wL7Nnz8bkyZP9btOjR4+QPttqtQIAqqurkZKS4n6+uroaGRkZIX1mpIg9bqvV2mqy9I0bN1BTU+M+PjGys7MBACdPnkTPnj2D3l+pde7cGTExMaiurvZ4vrq62udxWq3WoLZXolCOu6W2bdti4MCBOHnypBS7qBi+ft5msxnt2rWL0l4RIE+eXb58Gfn5+ejUqRO2b9+Otm3bhrvbqqHXfJRCOJm7bNkyLFmyBB988AHuvPNOKXczqnR1JTQpKQl9+vTx+4iNjQ3ps7t37w6r1YrS0lL3cw6HAwcOHEBOTk6kDiEkYo87JycHtbW1KC8vd793z549cDqd7oGlGEeOHAEAj8G4ksTGxiIzM9PjZ+V0OlFaWurzZ5WTk+OxPQDs3r076j/bYIRy3C01Njbiiy++UOzPNlK08PPWKqnzzOFwYMSIEYiNjcWOHTs0ewXKF73moxRCzdwXX3wRixcvRklJicecZk2KdmWUUp0+fVr47LPPhEWLFgkdO3YUPvvsM+Gzzz4TLl++7N6md+/ewttvv+3+esmSJUJ8fLzw7rvvCp9//rnw8MMPC927dxe+++67aBxCSPLz84WBAwcKBw4cED766CPh9ttvF8aNG+d+/ezZs0Lv3r2FAwcOCIIgCCdPnhSeffZZ4W9/+5tQWVkpvPvuu0KPHj2E+++/P1qHIMrWrVsFk8kkrF+/Xjh+/Lgwffp0IT4+XrDZbIIgCMLPfvYzYe7cue7tP/74Y6FNmzbCsmXLhIqKCmHBggVC27ZthS+++CJahxCSYI970aJFwl/+8hfh1KlTQnl5uTB27FghLi5OOHbsWLQOISSXL192/w4DEF566SXhs88+E06fPi0IgiDMnTtX+NnPfube/p///KfQvn174de//rVQUVEhrFq1SoiJiRFKSkqidQgUgmDzzG63C9nZ2UL//v2FkydPClVVVe7HjRs3onUYstNrPkoh2O/lkiVLhNjYWOFPf/qTx9+/5mMPLeEg1IdJkyYJAFo99u7d694GgPDaa6+5v3Y6ncK8efOE5ORkwWQyCcOGDRNOnDgh/86H4dKlS8K4ceOEjh07CmazWZgyZYrHX/7KykqP78OZM2eE+++/X0hMTBRMJpPQq1cv4de//rVgt9ujdATirVixQrj11luF2NhYISsrS/j000/drz3wwAPCpEmTPLZ/8803he9///tCbGyscMcddwi7du2SeY8jI5jjfuqpp9zbJicnCyNHjhQOHz4chb0Oz969e73+PruOddKkScIDDzzQ6j0ZGRlCbGys0KNHD4/fdVKHYPPM198TAEJlZWV0DiJK9JqPUgjme3nbbbd5/fu3YMEC+XdcBgZBEATZLrsSEREREUFnc0KJiIiISBk4CCUiIiIi2XEQSkRERESy4yCUiIiIiGTHQSgRERERyY6DUCIiIiKSHQehRERERCQ7DkKJiIiISHYchBIRERGR7DgIJSIiIiLZcRBKRERERLL7/ybFxxobt9LmAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1200x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.path import Path\n",
    "import matplotlib.patches as patches\n",
    "\n",
    "fig = plt.figure(figsize=(12, 6))\n",
    "\n",
    "# 円の描画\n",
    "ax = fig.add_subplot(131)\n",
    "theta = np.linspace(0, 2 * np.pi, 100)\n",
    "x = np.cos(theta)\n",
    "y = np.sin(theta)\n",
    "plt.scatter(x, y)\n",
    "\n",
    "\n",
    "# 楕円の描画\n",
    "ax = fig.add_subplot(132)\n",
    "L = np.array([[4, 0],\n",
    "              [2, 3]])\n",
    "L_T = np.linalg.inv(L).T\n",
    "x2, y2 = [], []\n",
    "for _x, _y in zip(x, y):\n",
    "    v = np.dot(L_T, np.array([_x, _y]))\n",
    "    x2.append(v[0])\n",
    "    y2.append(v[1])\n",
    "plt.scatter(x2, y2)\n",
    "\n",
    "\n",
    "# 楕円の条件のチェック\n",
    "# x.T H x == 1.0\n",
    "H = L @ L.T\n",
    "for _x, _y in zip(x2, y2):\n",
    "    np.testing.assert_almost_equal(np.array([_x, _y]).T @ H @ np.array([_x, _y]), 1.0)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "また、このような楕円を考えると、楕円の中の線形最適化問題はすぐに解くことができます。\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\min_x c^T x\\\\\n",
    "& (x - \\bar{x})^T H (x - \\bar{x}) \\leq n\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "については、KKT条件から最適解$x$は$c+2\\lambda H (x - \\bar{x})=0$を$\\lambda \\geq 0$で満たします。\n",
    "よって、\n",
    "\n",
    "* 最適解：$x = \\bar{x} - \\sqrt{\\frac{n}{c^T H^{-1}c}}H^{-1}c$\n",
    "* 最小値：$c^T \\bar{x} - \\sqrt{nc^TH^{-1}c}$\n",
    "\n",
    "です。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "さらに、\n",
    "\n",
    "* $\\mathcal{E}_*(X)$を「$X$を内包する最小の楕円」\n",
    "* 楕円の中心からの相似スケール：$\\alpha \\mathcal{E}:= \\{\\bar{x} + \\alpha z : \\bar{x} + z \\in \\mathcal{E}\\}$ \n",
    "\n",
    "とします。このとき、$\\mathbb{R}^n$内の凸体（内部が空ではないコンパクトな凸集合です）について、次が成立します。\n",
    "\n",
    "1. 相似スケール$\\frac{1}{n}\\mathcal{E}_*(X)$は$X$の内部にある\n",
    "2. $-X=X$なら、$\\frac{1}{\\sqrt{n}}\\mathcal{E}_*(X)$は$X$の内部にある\n",
    "\n",
    "つまり、最小楕円は凸体を中から抑える場合でも$n$しか小さくしなくて良いため、凸体を表現するのに向いているわけですね。\n",
    "また、このスケールはtightになっています。\n",
    "詳細は教科書2.4章を参照してください。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 最小楕円の見つけ方\n",
    "\n",
    "表記\n",
    "* $S^k$：実数の$k\\times k$の対称行列の空間\n",
    "* $S^k_+$：$S_k$のうち、positive semidefiniteな行列によるCone\n",
    "* $S^k_+$：$S_k$のうち、positive definiteな行列によるCone\n",
    "* $Y \\prec Z$：$Y - X$がPositive definite\n",
    "\n",
    "上で述べたように、「楕円はCholesky行列$\\sqrt{n}L^{-T}$を使って変形された円」とみなすことができます。\n",
    "よって、楕円の体積は\n",
    "\n",
    "$$\\text{円の体積 } \\times |\\det (\\sqrt{n}L^{-T})|$$\n",
    "\n",
    "であることがわかります（$M$でアフィン変換されると集合の体積は$|\\det M|$倍されます。詳しくは[Determinants and Volumes](https://textbooks.math.gatech.edu/ila/determinants-volumes.html#:~:text=The%20theorem%20on%20determinants%20and,a%20number%20A%20a%20B%20.)を参照してください）。\n",
    "\n",
    "このことから、\n",
    "\n",
    "$$\n",
    "\\operatorname{vol} (\\mathcal{E}(H, \\bar{x}))=\\frac{n^{n/2}\\Omega_n}{\\sqrt{\\det H}}\n",
    "$$\n",
    "\n",
    "です。ここで$\\Omega_n$は$\\mathbb{R}^n$での半径１の球の体積です。\n",
    "よってEllipsoidの体積を小さくするには、$-\\ln \\det H$を小さくすれば良いことがわかります（$\\ln \\det H$は$H$が正則出ないならば$-\\infty$とします。）\n",
    "\n",
    "つまり、$x_1, x_2, \\dots, x_m$を含む最小楕円を見つけるためには、次の最適化問題をとけば良いことがわかります。\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\text{ 主問題 }\\;  &\\min_{H\\in S^n} f(H) := -\\ln \\det (H)\\\\\n",
    "&x_i^T H x_i \\leq n\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "ここで、データは中心が$0$であることを仮定しましたが、これは一般性を失いません（2.3章を参照してください）。\n",
    "また、$x_i$は$\\mathbb{R}^n$を張ることを仮定します。つまり、$X$はfull row rankとします（こうしないと楕円の体積が０になります）。\n",
    "\n",
    "このとき、ラグランジュの未定定数法を考えてみましょう。それぞれの制約についての未定定数を$u_i$とすると、ラグランジアンは\n",
    "\n",
    "$$\n",
    "L(H, u):= -\\ln \\det (H) + \\sum_i u_i(x_i^T H x_i - n)\n",
    "$$\n",
    "です。ここで、$u\\in \\mathbb{R}^m$、$U:= Diag(u) \\in S^m$とすると、ラグランジアンは\n",
    "\n",
    "$$\n",
    "L(H, u):= -\\ln \\det (H) + \\operatorname{Tr}[H^T(XU X^T)]- ne^T u\n",
    "$$\n",
    "\n",
    "であり、$\\nabla_H L(H, u)=-H^{-1} + XUX^T$です。よって最適な$H$は\n",
    "\n",
    "$$\\bar{H}=(XUX^T)^{-1}$$\n",
    "\n",
    "になります。これは$XUX^T$がpositive definiteなときだけしか適用できません。semidefiniteなときは解が無限に飛んでいきます（教科書p.12参照）。\n",
    "\n",
    "また、最小値は\n",
    "\n",
    "$$\n",
    "\\min_H L(H, u)=\\ln \\det (XUX^T) + n - ne^T u\n",
    "$$\n",
    "\n",
    "になります。\n",
    "よって、双対問題はこの右辺を最大化させる問題なので、\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\text{ 双対問題 }\\; &\\max_{u\\in \\mathbb{R}^m} \\ln \\det (XUX^T) \\\\\n",
    "& e^T u = 1, \\\\\n",
    "& u \\geq 0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "を考えても同じ解が得られます。（$e^Tu$の制約はこっそり追加しました。詳細は教科書p.12を参考にしてください。）\n",
    "また、最適解の$H$や$u$が満たす条件は教科書2.2賞を参照してください。\n",
    "\n",
    "この双対問題はD-最適計画問題と呼ばれます（線形回帰の分散を最小にするような実験計画問題と等価です。教科書1.3章参照）。\n",
    "\n",
    "教科書2.3章（中心がずれた場合の最小楕円問題）を参考に、最小楕円を描画してみましょう。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fbbbf546e20>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiUAAAGdCAYAAADNHANuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVAElEQVR4nO3de5hU1Zk/+m91c7GjTQeVix66CWhGCN5AUcCZqEcmXiAzaqJkRsZEGSajGE3QnKBzEpNJFDwyjr94vxJD4qAZdYwgOEqixtDGVtDRX+SiTLoRRkRDugH5tUDX+WP3bqqrq/Zea++19rrs7+d5eLSbqq7dVUXtd7/rfd9VKBaLRRAREREZVmP6AIiIiIgABiVERERkCQYlREREZAUGJURERGQFBiVERERkBQYlREREZAUGJURERGQFBiVERERkhX6mDyBKV1cXtmzZgvr6ehQKBdOHQ0RERAKKxSJ27NiBww8/HDU14vkPq4OSLVu2oLGx0fRhEBERUQKbNm3CiBEjhG9vdVBSX18PIPilBg0aZPhoiIiISERHRwcaGxt7zuOirA5KwiWbQYMGMSghIiJyjGzpBQtdiYiIyAoMSoiIiMgKDEqIiIjICgxKiIiIyAoMSoiIiMgKDEqIiIjICgxKiIiIyAoMSoiIiMgKDEqIiIjICgxKiIjS2LMb2PlB8F8iSsXqMfNEROXWtG3Hf3+4C6MOPRDjmwabO5DWZqD5DmDdMqDYBRRqgKOmAVOuAJommTsuIocxKCEiI5IEFwuWv427X9jY8/U/njoa884eq+3xqt6v5X5g2TVATW0QkADBf9cvB9YuBab9CzBxlvBjEFGAQQkRJZb0RJ8kuFjTtr3XfQDg7hc24sxxw2MfO2kwU/F+n/tTEJCgCHTt7X2H8OtlVwPDxjFjQiSJNSVElMiC5W/jvDtXYe6jb+C8O1dhwfK3he5XLbhY07Y98n7//eEuqe+nfbxq99v+q38NMiRRamqxfeW/4vHV78U+DhHtx6CEiLCmbbvUCTTpiR5IHlyMOvRAqe+nfbxKfz8Qn+DTrc/2zZCU69qLQX94Btc+2iIVsBHlHYMSopxLkvFIeqIHkgcX45sG4x9PHd3re5edOjp26Sbp41X6+3rsRgFdkfcL1RaKqEfQkSMasAHyASKRT1hTQpRjSes0kp7ogf3BRenjigQXADDv7LE4c9xwqTqWpI9X6X5/8xdjgVdr9he3RthXLGAH6nq+/u8Pd2mrfSHyBYMSIo/IFp5GZTyi7p8msACSBReljy3bCpz08Sreb8e0oMsmYglnb7EG/9l1AjoxoOd7cQFbmkJeIl8wKCHyRJKr7DQZjzSBBZAsuEgj6eP1ud/kOUHbb4TaQhce2HtOz9ciAVvSABGwaHYLUUoMSog8kPQqO23GI+vAwgojJwdzSJZdHXThlGZMavoBXftQmPYv+H+HnS8VKCQNELnkQz5hUEJkoayWYYD0GY9cmjgrmEPSfEeQNemZ6HpOkElpmoTxgPSylGyAyCUf8g2DEiLLZL0MA+Q045FW06Tgz57dQOcOYGA90L8u/n4RZAPENMEokY3YEkxkkaTzP5K2y5IC/euAg4amDkhC45sG4/wJI7Qu+YTYfky2YaaEyCJchiEZaWqCWItCNmJQQqSRbG0Il2FIVpJglLUoZCsGJUSaJLkSTdsNQ/kkG4yyFoVsxaCESIM0V6JchiHdkmbkOA+FdGNQQqRB2itRLsOQTkkycqxBoSwwKCESJHOVmLY2hEg3mYwca1AoKwxKiATIXiWyNoRcIJqRYw0KZYVBCVGMpFeJrA0hX6TJ/LEOhWQwKCGKkeYqkbUh5IOkmT/WoZAsBiVEMVgfQiSf+WMdCiXBMfOUSzLjtTnCnSggMwI/KsNIVA0zJZQ7SVLKrA8hksMMIyXBTAnlStIN7wC5q0SivGOGkZJgpoRyha2NRNlJUofCbGS+MSihXGFKmShboh1o7NQhgMs3lDNMKRPZJ82yKvmFmRJynmzKl0WrRHbhsiqFGJSQ05KmfDnUjMgeXFalEJdvyFlM+RL5gcuqFGKmhJzFlC+RP5Isq7Jbxz8MSshZTPkS+UVmWZXdOn7i8g05iylfonzi0q2/mCkh68ikZNlJQ+SYPbuBzh3AwHqgf12iH8GlW38xKCGrJEnJspOGyAGtzUDzHcC6ZUCxCyjUAEdNA6ZcATRNkvpRXLr1F5dvyBpMyRJ5quV+YNHZwPrlQUACBP9dvxx48Cyg5QGpH8elW38xU0LWYEqWyEOtzcCyawAUga69vf8u/HrZ1cCwcVIZE5mlW3bpuINBCVmDKVnKJQU1FlZrvgOoqe0bkJSqqQ1uJ7mMI7J0yy4dt3D5hqwhnJLdsxvY+UHwXyJXvfs88LMLgBsOAxZ+FrjxcGDJTKDtZdNHps6e3UENSVRAAgR/v3ap8n/TXBJ2DzMlpJ2ybhqFhXJExrQ2A09/G9j6Zu/vhzUWa5cC0/4FmDjLzPGp1Lljfw1JnGJXcHuF2SIuCbuHQQlppaybpuX+YF26prZvoZxPH+Iu8n35QaWW+4P6iWpS1FhYaWB9cPEgEpgUaoLbK8QlYfdw+Ya0UZY6jS2UKwYf4j6lvV3Q2hwsN9x4uL/LDyr1vI8FhDUWrutfF2Qza2Kuf2v6AWOmKw9q2aXjHmZKSBtlqVONhXKUEDNX8prvAGpqgK598bctrbFwPfs0eU7wu0Tp2hfcTgMOWHQLMyWkjZLUqeFCOaqAmSt5Pe9jgYAkFNZYuG7k5CBIRaFvxqSmX/D9af+i9YJifNNgnD9hRGxAsqZtOx5f/R4LYQ3SGpRs3rwZM2fOxCGHHIK6ujocc8wxePXVV3U+JFlESeo0SaEc6RVmrqL4svygisz7OKShxsKYibOAS1cAR50T/F5Ad6H6OcH3LciqLVj+Ns67cxXmPvoGzrtzFRYsf9v0IeWStuWb7du345RTTsHpp5+O5cuXY8iQIdiwYQMGD2bqLE9Sp04NF8pRmfCKP+718Gn5QQWZ9zEAFGqBMdP8eu6aJgV/LCyMrlb/dua44VzuyZi2oOSmm25CY2MjFi1a1PO9UaNG6Xo4sliqvWnCQrn1y2NqSvoFV12WfMh5y3CLp7NE38ehYpe2Ggvj+tdZ955g67A9tC3f/PKXv8SJJ56ICy64AEOHDsX48eNx3333Rd6ns7MTHR0dvf6QnTJde508J34tXmOhXKS8DXILr/hFMHPVm8j7OKS5xoJ6Y+uwPbQFJRs3bsRdd92Fz372s3jmmWdw2WWX4corr8RDDz1U9T7z589HQ0NDz5/GxkZdh0cpZL72akGhXB95bYc13OLptKj3cWj4scClz1hRY5EnbB22R6FYLBZ1/OABAwbgxBNPxKpVq3q+d+WVV6KlpQXNzc0V79PZ2YnOzs6erzs6OtDY2Ij29nYMGjRIx2GSpDVt23Henav6fP+Jy6fo/wfc9nJQPLl26f6JrmOmB1egWQYkpe2wpan4mn7BlbDv7bCtzcGOr4j66CgEBYy82u+r0vv4yL8M3sejTzV9dLnGjfvU6ejoQENDg/T5W1tNyWGHHYbPfe5zvb43duxYPPbYY1XvM3DgQAwcOFDXIZECRtdebSiU07TjqVPCK/5lV0cHZr7+/mnZ8D6milLVv5ES2oKSU045BevWrev1vfXr12PkyJG6HpIyYMXaq8lCOQ5yC0ycFQRe5Vf8R52TfebKVRYWfJIYZlT00RaUfOtb38KUKVNw44034sILL8Qrr7yCe++9F/fee6+uh6QMhGuvpe1zuVl7ZTtsb7zipxxKsp8XidNWUwIAS5cuxbXXXosNGzZg1KhRmDt3LmbPni18/6RrUqRfLq8Udn4QFLWKumYDcNBQfcdDRJkyWlPnGOtqSgBg+vTpmD59us6HIIVkAo1crr1ykBtRrnGeiX7ckI8AMCUphIPciHLNipo6z3FDPqo6YpmbUlVg8yA3ItKK80z0Y6aEmJKUwXZYolxLvZ8XRWJQQkxJymI7LFGuidTU5bIZQAEGJeRem68NLahsh7UbXxdxfK6UY41eclpbgtNiS3C2rI/sW5uD7EQ4K6RQExSeTrlCLjvBD2F/qXqP5AGfKy3YNhywsiWY3GJ1m2/pfjNhS26xK+iEWbtUbL8ZFz+EGUCJU/EeyQs+V9qwRi8ddt+Q/WL3mykGhadRO/S23B9sIrd+ed8P4QfPAloe0HX0yeR1F+KkVLxH8oLPlVas0UuHQUkOrGnbjsdXv+dui2+430yUcL+ZSlz7EHYtgLJB2vdInvC50optw+lw+cZzzhdcqdhvxqVN9LgLsTzuSSSOz1Um2DacHDMlHvNiKFrnDrGx7kBwu84dvb8XfghHBSRA7w9hk3gVKy/teyRP+FxlZnzTYJw/YQQDEkkMSjwWVXDljHC/GRGV9ptx6UPYtQDKFmnfI3nC54osx6DEY14UXIX7zdTErDTW9APGTO+banbpQ9ilAMomad8jeSL8XNXyudLI+To/jRiUeMybgqs0+824dMJyKYCyjY49ifbsBnZ+4F9GSvS52vWRPcXfHlmw/G2cd+cqzH30DZx35yosWP626UOyCoMSz807eyyeuHwKbrnwODxx+RR8x6Ui11C43wwKfYOLmn7B96P2m3FlEz2XAijbpH2PlPK9HTvquSr13u/Y7aWYF3V+mjEoyQEvCq4mzgIuXRHsLxNmE8L9Zi5dET3oSeUJSzeZAMrXK/mk0rxHQnlpxw6fq8aTq9/GxnZ5x3lR56cZW4IdZ/1oeJXS7DfjyiZ6IrsQT74CWHW7W5Nps5LmPZK3duymSUDdwUChFihGBMK2tMt7wIs6P80YlDjM+RkkSfWvS7Z04comelEB1KebgObbk40Ht/33VinJe8SleTYqcGZJ5pzb/NQABiWOqrY2eea44XyDx0ka1GSpUgC15fVgaUH2St7FPX+ylscTdJJuL9d/ZwtwsFo01pQ4imuTOdG/DjhoaPDfJIPV8lIjkVYe27HZ7WWMF3V+mjAocRTXJnMmyWA11/b8MSmPJ2h2e5GFGJQ4ypsZJCQmyZU8R9aL2/I6cODQ+Nv5doJ2pV0+Z/I8XI01JQ7j2mSOhFfyIoFJoSY4eeatRiKplvuDjJJIpsS3E7RIt5ct7fI5kdsGhm7MlDiOa5M5IZtq79qbvxqJJEqXuKLaYgu1sGqejUoq5ruQEhyuxkwJkTsmzwmyGlHCK3nZzIoPNRJJiLQBA0Gx8QU/8S8gCbnSLq+DRb9zVANDXi48GZRYLFeD0SiebKr9qGlBl03k3I1+wRVxXk5ApUTbgAFg51bgsOP0H5NpLrTLq2JhqzwbGLh8Yy1u2uSgLMa+y6TaWcQYLY9twBSwtFWeDQxAoVgsFk0fRDUdHR1oaGhAe3s7Bg0aZPpwMrOmbTvOu3NVn+8/cfmUXL05nWHqiksk7dzyQHxmJa81A3t2B5vtiS5xXbclP1kEn7U27x9CWFUhCPINZUx8yJInPX8zU2IhDkbLgKqshskrrtLBatWwiLE6zunIJwda5fPcwMCaEgtxXVGjJFmNahkJVzZwy3MRYxyZ4uEk+JzbJY/bCTiGmRILcV1RE9msRmszsGRmkOJf+Nngv0tm7p+A6sAVVy8imZW8CYuHUeibManph8RtwHHvHTLD8TqiPAxVY02JxXxYV7SG7DpyOFCrWi3GWfOBZ65jPYIvmYC2l/vuyjxmepAhkQ1I4t47ea7jMc3hOiLXhqolPX9z+cZi45sGMxhRRWZb+mIxfllmxTzxx/Zxh1UL2ylTiVrikgm8XFnSy6uwjsixVvk87QrPoIT8J7uO3LU3PoAp1EZPAO11Ww+Gk5WemF//+f5MQPky2NqlbmcCSud0JAm8ZIJfBiVm6K4j0iBPQ9UYlBjE5ZmMyK4jr18Rf/swIKmpjZ4FYtkVl7RKJ+bwufE5E1C6BCMaeLGI0g0O7veTp+YHBiWGuLY+6DSZkesoiAcwANAVdwKy64pLSrUTcxzXMwFJl2CSFFEyKDFj4qzg9SuvIzrqnGR1RJqFzQ+l5wxfmx8YlBiQp/VBK8isI3/2zN7dOVEKNcCZ84P6EkeuuIRFnZjjuJ4JSLoEk6f9hnwocHasVT4vu8IzKDEgT+uD1hBdRz7lyuCEIVoIN+kfgcOPd+aKS5joRnXVuJoJSLME42gRpRTfCpwBp/b7yUPzA4MSA/K0PmgNmXXkYlGuEM6xK65YMhvVVeNqJiDtEoyDRZTCktTZEEni8DQDOBzNENGR60kHavkynEzmxFyJy2PZwyUYEZUCL13D2EyLrbMpBgE/h8MZ4dNQNWZKDMnL+qB1RLMajhXCKSVVGFyBq5kAQM0SjI/vHbY6W8u3pglOdCWK48uyjIwlM+NPzOV8mViqchdZH947Dk9B9Z3NO8pzl2CiUqp2AQbEl2VUPqZpk+dEz1/pUej+j0c7D6tcgvFhSc/x/WJ85uOO8ly+Ib+Y6A7wsSNBtDD4+L+1PxOQJFvh4xJMUnlqdXaMj00TXL7RjFNbM2RiIzTfN19TuVFd1lQFiz4swaQlspwX1tnMWJzdcVGfmpLLTh2N71hQU5L0/M2gRCPfCpCsprIOwObHNMW1E7PvwWLW8vRed5CNF7+sKbFMtamtPrRsWSnsDogSdge4/JimuFQbwfZV9XxtdfbE+KbBOH/CCGsCkjQYlGjiYwGStcJhX3GdIuEUzo//mP1j+lD86oo8BYtZEp3zkwWfisqpFxa6auJjAZK1ZLsD/r/RQV1EmkJUbr5mJ+7Uq1f5nJ+afsFzmVVxq49F5dQLMyWacGprhmSmcAIAikHR3oNnAS0P6H9MmzoSfL/CZPtqNFWv/5bXgaVzgZuPABZ+NphjsmSm3iWxlvuDupbSDTPDMfdp/i17yOUJr8yUaMSprRkRncJZKmoLepWPacvma3m5wpSaRluwJ1jUTeXr33J/8O8GBfQUvha7gszT2qeAabeoX8qJrRNC8n/LnnG9wYKZEs18KkCymvCwrzJpagtEHtOGkes2XGFmlaEJg8VCTE0JAKAI/M8beo/HBipf/9bm7oAE6NuJ0/31srnqMyasExLiQ4MFgxLyQ1R3QJQ0hagqOxJET9qyJ3fTnSitzUFa/8bDs0vzT54DFAUC1EIOTmKqX/9f/bPY7Vb+UOYoo7GoXJgPDRZcviF/VJrCKSJNIWrayZ+iafWk6XeTG6mZ2ur+8OPRa2mhmuI++4td086HUfn679kNtPbdZ6Wi1pfUPa8sKhfmQ4MFgxJFbBxek0thd8DHfwy6bOJOTED6QlTRnYfLiZ60k57cTXaimKwB6NwBodcdsPckpqIGRPXrv3Ob+PGHtx/cJHefSjjmXljYYFE+4dWlcxKDEgVcLyzy0qcODtp+syhELQ1GDhoqdh/Rk3axCDyd8ORu8grTZIbG9ZOYqgyTLxkGk0Xlrk0yhvsNFqwpScmHwiJv6S5ETVMvIVq495ubkxf4mWpb1lUDIFpPE57E4mqLavoFgatNJxuVNSCqX/+Dhoj9rKS3j5J1UbmJWiiFXG6wYFCSkg+FRd7SORo7TUeDzEl7x/vJT+6mTs6qZ4UkOUG40hlVTmWXierXv38d0DQl/nFDj81WdxLPcsy9Dd1qOcagJCUfCou8pmo0dulVetqrWZmTtqhqJ3cTJ2eVV+hJTxA27NUi2ymlI8Ok+vU/43titwPUn8SzGHNvuluNWFOSlg+FRd5LWogKVC44PHBod81CxId9VL2E1IAvQdVO7uHJednV0Tvmqjw5h1fo656OeY5iagDSFsum7YxKKmmRqo4aENWv/8jJwXC0ZXPjb6ujoDnNv2URJmuhNHKpEYNBiQKuFxblRv86uQ+wagWHO9+Pv29UR4NM4d6BhwK7PkxX4Gfi5Dx4ZPyskLgrdBUnCN0nsXJpilR1Feiqfv1Lf97bv4y/vY6TuOy/ZRGe7pvkWiMGgxJFxjcNZjDik6irdFFRV7OT5wQfbFG69gF/8e3u7puY28Wl37M8OYeZgjiTIzIHqk8QOk5i5d59fv+00ySZHZ1dJqpf/6ZJwGHHATc8hdj2a1dO4r50K5Wo1ohx5rjh1p6vWFNCVIlIwWGcqKtZ0ZqHk/5ebW1E/7qgbVnnh6nIc1eoBf7UWv3vXdpYLyzEXfzX8beNK1LVXQOk6vVvbQYeuRjS82Bs5uommxFcbMRgUEJUTrTgMIpIR4No4V4WBX6qiD53pdNUK3HlBNFTiPu02O3jilRFt0tYdbu5Ysvwd964Uvw+LpzEXW4lr8LFRgwu30hyqWAo95KmqlV0x1S7mi0/JtG0eta1EUmpSoG7sAtzryU+ic0g41L/vWo2Ki2PFPWP6q8mye9sy07ZIkSXVW1rJa/CxUaMzIKSBQsW4Nprr8VVV12FW2+9NauHVcq1gqHcSjuiO013TLWOhrhjEq15yKI2Ig2VxZpCJ4i9wMS/lztGVUQKcSsRyRo0TQqm+VYrJNU9qr+aJL+zQydxI91qmrnWiJHJ8k1LSwvuueceHHvssVk8nBac3OoIFYOPRNO4hVqg/rD4JZU8DWNSmQIXXcpYfG720zaTLvHJpP6b7xB4HjPc6Vj2dy7UIpN5MKq5tFwqyKUJr9ozJTt37sRFF12E++67Dz/60Y90P5w2UQVDLrzQuaByEziRq/RiF3DBT4IuhGpLKqLHdMiRwNCx9i7LyFCZAq+683PJLsBZ7DxcLukSX9desd/bxvZU2d/5yDOAv7jarYAk5MpyqYe0Z0rmzJmDadOmYerUqbG37ezsREdHR68/tnCxYCh3VI7olpkIGtXRINTFUwR++ldO7rFRkeppqk2TgBmLgb97IrgvgD51FllP25QpxE3Cxu4j2eLjC3/qZkBSKotuNepFa1CyZMkSrF69GvPnzxe6/fz589HQ0NDzp7GxUefhSQkLhkrZXjCUK3t2B1eMKkd0p03jJknx+7Ksk+S5ixvL/sr96oLOtESXqcrV9BM7Phu7jzzsTsmjNW3b8fjq96wtPdC2fLNp0yZcddVVePbZZ3HAAQcI3efaa6/F3Ln7xxd3dHRYFZi4VjCUK7+7G9IzE0Q+NNOkcdOk+AH9RYy6U9Oiz51IYbKNyxkiy1TlZIa92dh95Fl3St640KyhLVPy2muv4YMPPsCECRPQr18/9OvXDy+88AJ+/OMfo1+/fti3r2872cCBAzFo0KBef2zjUsGQE2Q3LauktRl47gfit09yZZkkjZs2xa/rqj/rbdmjnjvRImAblzNKl6kKEoP2RI/Pxp2ObdjokBJxpVlDW1Byxhln4M0338Trr7/e8+fEE0/ERRddhNdffx21tSmnZZLbVJ4YZaavZplaTpriD8ksNYmyqRNIZkdWG5czgP3LVH92lvh9RI/P1gDAw+6UPHBluqu25Zv6+nocffTRvb534IEH4pBDDunzfcqZNJuWlRNN64dEux9USZLiL6Vyjw2V3UkqyGy4N2OxncsZwP5lqn/7G2D9M+l2Ri5naqfjOOxOcY4rzRqc6BqDE1wVU31ilK3bmPqDbD/Io4YxiVB51W/TtuxJakRsr2eYciWwbnn0bZIcn80BgO3D/KiHK9NdMw1Knn/++SwfLjUXioKco/rEKDtB9OSvix+rKiKzNipRedVvW6FokhoR26dt6j4+BgCUkgvNGsyUVOHils/W03FiTNqlkPVVZ/nV7ge/B356bvR9VF7127Yte9Jx9LYuZ4SyPD4bMydkvfFNg60+hzEoqYITXDXQdWKUSeun3RcnrfBq96Ch2V71q9yTRoU0La82L2cA+o/P9HuYSKNM9r5xkStFQU7R1UEh2qWw9S17Ok+AbLsYbBx8lbbl1fZpm+HxAenb3kM2dU8RacBMSRWuFAU5RedAqLi0ebEYfJjb0nkSyvKq37ZCUdtrRNKKy2jIvua2dU8RaVAoFouCYzCz19HRgYaGBrS3txsbpMbuG8Vam/cHB1UVgkxB0g/WSh/2S2aKB0MzFid7XBe0PBAfBGQ9Z6Lt5b7B5JjpQXAUtdmhzUrb3vs8z3uBYccAH/xvueUXvodJId3ntqTnbwYllL2sT4x7dgfD2UTrKa7b4tYJUFZUEGDyCrs0mNzyurt1E0KBd5m4934W72Fba3RIuSw6S5Oev7l8Q9nLuoPCts4T02wtFA2LgFUO1zNBpO29XNzyi873MAtnc8X2zlIGJWRGlidG2zpPbGHj3AtddRNZBWCyE4bLVZvRo+s97HoASNJs7yxl9w2ZlUUHhY2dJ1SZyD5GMhsVZr35YNKdoUPhjJ7tbb27dXrew3HPjcR7WGbvIfKG7Z2lDErKrGnbjsdXv2fdzomUko07rlJvYZYhbtlDdKPCau2z657W1z6bdmdoIDjG/3VM7wCqtRn4+CO172HVASA5IewsLWVTZymXb0pwrLzHfG8/9YHKuomoLEC4Yd6yuerbZ0Xb3kX0LKM8FXwdlemTfQ/btu0AZcrmcfPMlHSrVvzDjIlHZIaV7dmtbuAViVE5XE8kCwAAT39b7PFkiGTlRJUGNlFBTuMkuYF7SQJA8sr4psE4f8IIqwISgJmSHrYX/5AicQW27EQwR9VwPZli0/f/C9j4AjD61GTHXEnanaFlFWqBTx0s9/40VfxtW8cXWYeZkm62F/+QYpUKbDnC2zwVtT+yxaY6aiaqZeWGHxv8f1zRtYziPrEam1JZF39nXXBMzmJQ0s324h/SzIZOBC4Zie9jFJUVGFgf3E7UO8/qec6bJgWTVa/bAlyzIfjvP/4GuPSZ3sGKCkmWWLIq/mawby0bGzu4fFPC5uIf0kxk4FW1GRJpccmot7TD9frXAUf+JfDOf4o9nu6BeeXzYMIlxIf/Bli/AkCKFuJQkiWWLIq/uV+PtWxt7GBQUmZ802AGI3ljshOBw6sqSztcb/Ic8aDExMC8PbuBDYoCkiQbWIZ0T1c2GexTVTZPdWVQQmRqDD2vIuMlnTp7xGnBpndb34y+XZoTehpph6yVSrvEomu6MtuOrWVzYwdrSkgNl+shVLaiyvB9eJXp98Q5N8ffxtTAvCRD1pLW2IhSPV2ZbcfWsrmxg0EJpeNDVb2JMfSqp5faxJb3xMjJwLRbEBTNlgV/5Sf0rAMo0fccAIw8pW9xbLX5OjYxFexTLJsbOwrFYlFif+1sJd36WNaatu0sbk2itB6iWpGcrR+Y5YS2my8EJwEVV6U7PwhO2KKu2RBcxdrOxvdE28t9aybGTA8yJMWiuSLjJO+5LOd8qHisJTPF587MWJzsMSgxnee+pOfv3AcltlYgWy/rk3gpXR/MLQ/EdyKoOqHu2R1kEESHV123xcx6u8xzbfI9IaL8d7EhgMryPSdiz27g3V8DqxcHhbhpAzXb3xOkTdLzd64LXW2uQLaeiap63a2zujsRSqmaXqpLkufa9k6L0qJZW4qMs3zPRQlf77VL0SeASNMNxj2nSFKugxKbK5CtZqKqXrZ1Nmk2RVcnQiWT53SfBCKYKMRM0qbsWqeFTQFUlu+5SsLXu1CDqhmNNIGaLYEXOSHXQYnNFchWy7qFVuaqVlWNQNJWVBk2XkUmzSCYaqtOwtYAKnzPhUW3WQQnpa93uHtylKSBmunAi5yR66AkrEAuXcKxpQLZallv5iV6Vfv0NcD7b7k1iMy2q8gkGYQ9u4E9/8fMBm9JqAygVJ5kTUz2FXm9S6UN1LII9kmKbY0euQ5KAI6WTyTLeog9uyuvc5fr2gu8/+b+/y//O8DeQWS2XEXKZhDefT5I/fe6TwGRr5WpGplSKoJq1QGEicm+MrsplzKd6SJlbGz04JwSBBmT8yeMYEAiI6vNvH53N2IDElG2DyJTPbxKlmwGYfFf995kLfiL6PuZGlZWKu1cGtUbzJnaDDLpVFnTmS5Solqjh+nN+RiUUDIqdnON09oMPPeDNEfZm6pBZKYnleqSZMqoaNpf9fTRtJIG1ToCCFOTfZNOlVU1QJCMimr0MIlBCSU3cVYwX0DXpEmRD2tZacZZ2zKpVBfhKaMFuZ9r4/TRpEG16gDC5GRfmamyPcdhQaaLlLC10SPXNSW2Ffg4yfRmXrKSpp7zspuvSJuy8HJaAbjqv4CDhthzZV36PpUtMtbRtWO6a0no9QZQ6H7f25LpotRsbfTIbVBiY4GP01RX1cuudx86BvjjO2oKb8sDLFsGbekWFm9WE7YpCwclRaD/AXoDEtFgOKowdcZisZ+jI4DIupOt3MjJwFkLgBXzgJqaKktaBWDMNM4U8ZCNjR65DEo4ydUBsh/WZ98ELD43+nZxqedqJ66P/2jPoK1KVGSpSjNBFYOOQhDQTfz74Hk23for0/0imuUyEUCYnOxb+hyi2DsgKdQAR50NjP87YPRp9mS6SLnxTYOtOu/lsqbE1gIfKiHbIXHEaekKb6t2VDwNtP3Wzt18VdW4RGWCehSDgG70qdnvqFxOpvtFZWGqrt2ks+pkK1XpOQSCZRoAOPNG4CsPB4EJAxLKUC6DElsLfKiM7Id10sLbyBOXwJTLUJoiWlkq21KFijf77V/aMXESDckGGaoLU3X87ll0spWKeg7Dqa4rrvWngJucksugJCzwKWVDgQ+VSfJh3TQpqBG4bgtwzYbgvzMWR3+gq+ryyWp+g8qr/yTdH1mfREvJBBk6Olt0/e66O9lKmWpBJus80tKG6x7/LzzS0mb6UHrksqYEsLPAhypIOoZdtPBWVZdPlpNKVW4ml7R408R4fNnul53b9HS26Prds5jsm9W+P6anE4tw4Rg1OveOl/D6pnYAwMOvbMK/vdKG/5jz54aPKsdBCWBfgQ9VofPDOulUy3JZzW9QfVJJU7yp6nURvb9sABUes46iXJ3vSZ37w+huQTaxf48sF45Rs0da2noCktDrm9rxSEsbZkxsMnRUgVwu35CjdIxhl51qWShLe6dJ2SeZDJvkpBJFRfFm0tdFtlBX5rUq1ATzUXQX5ZreGkCW7HMoE6ipHr+vgwvHmIE3Nv1J6vtZYlBC+SZzUh7558G8hrRr/mm6ZnScVFQUb8oGWC/f1X1yeFr85BC+VuWBYbnSIMNkUa6NdHUQmdq/R4YLx5iR4xo/LfX9LDEoIRI9cZ3xXfki2nJpr9R0nFRGTgaOPj/6Nkd/qfLvKRtgtTYDD54dDOsqn40BRJ8cWpuB3R/t7xCppjTIMFmUaysdgZoLxbMuHGNGZkxswvGNDb2+N76xwfjSDZDjoGRN23Y8vvo94zsikgVkT1xplitUXKkJnVT2BkOvRI/rrcejb/PWY32PSzbACm/ftir+mMpPDuF9N70ScZ8qQUaWnS2iTG7qqDpQM7l/jygXjjFj/zHnz3HTl47B357UiJu+dAyesKDIFchpoStHzOdQXEFiFt0kIl0zALDyh8Aly6r/fXhSWXY1UChUrjEp1AR/D8SfdJN088iO3i+9vYjSk8OW1wWGuwFoPBk443vAYccFJ/zS1zqLzhYRYZHl2qUInotCkNHKushS5fvd9P49Ilw4RgNmTGyyIjtSKndBCUfM54xMpb3OE5dw63ERaH0pqLmYdFn1m02cBRSLwNNXV/kx3Y8TtydP0m4e2UBGNCAr/x06d4g/FgCsuj36tdbZ2RKn5f79gWKPIrD2qeDPtFuyzdqoer+b3r8HiP8dbDhGEpK75RuOmM+RpPUbOjoqZFuPV8yLX8bZ+LxAbUnMGnmSK0jZVPjHfxS7fblCTfD7CT3WPqD1t3KFs1lqba4QkJRZNtdMkWXa97uu4lkRojVNJo+RpOQuKOGI+ZywrdJ+YD2AgvjtCzHBhKo18iTdPLKBzI735WfBhCeHrr1y95UpnM3Syn9OdjuTtScybNm/JyoQZSeWE3IXlHDEfE7orrSXPVm8/nMI11MAQYdJVDChal5JkitI2UCmfrjcLBhg/8lBdo5MNSa7KvbsFivuBYJsz57d6jZbzIpN+/dUC0TZidXD5kaP3NWUABwx7z2do7STTIPs+QCVFFVwp3KNfPKc7sLLCKVXkGEgs355TJ1H9+j9Tx0sdvvgYIP/lJ4chO8bdfwpx6aXkq3B2LlN7uevug349Y1BIFWeAVi7NHhuTHQMxcly64GkWy2Y2B7BMrY3euQyKAE4Yt5ruirtW+4PggvZk0WSIk8gOpiQDQyifr/Sbp7y46zpFwQk5VeQsoGMyO0BYOQpwTwY2ccSkbarIqvx5L++IfivSFeTbVzYv8eWTiwDXGj0yN3yDeWAjqmnSWtURGs/yokU3KlcI5ed5SGbCo+6fTih9awFQSt0+clW5L4i0nRVpBl6d9AQuceKnVgrsBRluhZF5/h9lUuXLm0RoIALjR65zZSQx1RmEUJJ08VJN/wTCSaSZDiiyF5ByqbCq91+zLT41HnUfT/+CNj0O3WvdTnZmSzl+tcBI6cArYJ1JbETayOWovKw2RzbexNzodEjt0HJmrbtrClxiWyqVXZ5Ie6xk6aLZT5AgeAqudglHkzoWCOXmeUhG8ikSZ1Xu29rc5DFiJKmqyJpQFrq//4esOisZI9fSaWlqKTLi67RcdGRE2GjR+kSjm2NHrkMSmwv9KESSa/8VGYR0tSoiH6AhkSyBuVsWCOXHUqWZohZ+X1VZ4xKqSqaHjk5GI62bC6CYt7STqzur8+6CXjm2mQZgLTZnJArdRYqLzpyxvZGj9zVlFQr9LGxNSr30m5ep2rPk7Q1KiK1HygAF/9SfoO/UjlcI++hen+bsCZj5zY19Qs9x/gMMPaLvY9x7BeD70/6x+QDvtK2wLMFOVfGNw3G+RNGWBeQADnMlEQV+tj4AuWWqis/FVmE/nXA0HHA1jejb1ctXSx6JT/6VLnjot5UvNblmTmpgXcC9Qtxx5gkA5A2m+Pqsg/be72Uu6DEhUIfgpp1/FJplgtam+MDEiA41mrpYn6AZifpa13p5Cw68K48II0LjKodY5KlqDTLi6qCf1NsWLokpXIXlLhQ6JN7OoefJdF8R/cJIaYmZPix+rM2pEfUyVlEmL1Q0f0iG8Cm6UZRHfybYnKjRVIqd0EJYH+hT+7ZtM248O6+ALa+JRYglX+AMkgxT3jAXVmRamn2YutbyZZBKr3+MgFs0m4U24J/0s6FrtNcBiUAJ7pazaY5BDoDpDzMlBBhOiiTCTxR3P/eLM1eFIvdbckSyyAir79oBiBJLYpNwT9p50rXaW6DErKYTXMIdAVIJosLTQcBIVuCMtkBd1e+AfQ/oPfzt2Sm3DKI6tc/SS2KTcE/aeXCePlQ7lqCyRG2bDOeZBfdOElH1qdlU9tn2nZvlWRbvg8a0rv1WnQrgXAZ5N3n9bz+sm3ROt7bZCUXxsuHGJSQnWyaQ6A6QEo7UyIJm4IAU0FZNWlPzrLLIDpf/6ZJwayb67YA12wI/hvOvqm0H44twT9p5VLXaW6DkjVt2/H46vc4NM1mqgdiJaUyQJK9qlaxoZptQYCJoCxOmpOzbKblnWf1v/6lg/SqZcg2vgAcckSwEaINwT9pE3adlrK16zSXNSWuFPwQ7GmjLW3TfPspBB0YBfk5IyaKC21q+7S14yPNqHqZGqgjzgA2PCN2TCpe/2q1K2uXAmufCr4u1ABNU4BCAWhbxRk6nnKl61RrpmT+/PmYOHEi6uvrMXToUJx77rlYt26dzoeMxTHzjrJhhHqxGPwpdE/5LBSCr2WkHVkvy0RmJoqqbed1SJOZk8m0ZPX6R85eKXnfFruA934HtP4WOPPGvss+5A2bx8uHtAYlL7zwAubMmYOXX34Zzz77LPbs2YMvfOEL2LXLXHGNSwU/pFmlNfZqVNVkZF1caFsQkHVQJiuqJiOK6BLf6FOze/1FlslC4TLeimuBP25kUSsZo3X5ZsWKFb2+/slPfoKhQ4fitddew+c//3mdD12VSwU/VIGKZRzZVlTVo7iz3OHUtrZPm9q9oySZECo6iTWL119q9koJ26a3ml62pcxlWlPS3t4OADj44IMr/n1nZyc6Ozt7vu7o6FB+DBwz7yhVMy2SzIdQXZORpn5BVk8Q8HT08kKWQYDP286L1EBl8frLzl4J2TK91ZYZNh5wYYprqUKxKLsonkxXVxf+6q/+Cn/605/w0ksvVbzN97//ffzgBz/o8/329nYMGjRI6fG49kLlWmkgUe0DXKQTp7V5/9TNqgpB/UD4wbdnd9CtIJppuG6L+Id528t9r6rHTFdbXNjaDKz856CAMVLZ761bywPxJ2Ubd6ZVSefrL/O+reSaDUENVxbKAzhV/97JaFNHR0cHGhoapM/fmQUll112GZYvX46XXnoJI0aMqHibSpmSxsZGLUEJOSJJIFHNkpniywYzFgdf7/wgaKMUleTDXFeKutqHeymTH/RZBGU2Kn+9db3+Iu/3SmSD66QqZUMaJ9kXQDtqTdt2nHdn3+fyicunZHIhnjQoyWT55oorrsDSpUvx4osvVg1IAGDgwIEYOHBgFodErlC1dJK0FXVgPfpswlZN0poMHTuciu562zgJOOO7wGHHBQFY6YlR93q+Le3eWcl6SUJkmaxcVst41ZZRYwMS2Ff3Yqmopg6bVwe0BiXFYhHf+MY38MQTT+D555/HqFGjdD4c+UblTIuk80G2vA6xgKQWGDPNnpOqSDBXqAVQBFbd3vdqtVAA2pq5nq+Kib2OompXqsmilkc0YK7GlroXy7na1KE1KJkzZw4efvhhPPnkk6ivr8f7778PAGhoaEBdnR1vJtaWWEzloLGkXSjNdwQn72LMDIqiRYWZosFccV8wm6IQc7Wq6+SZl2JG1d1bMip1BFWiusA6ikjAHIe7FsdytalDa1By1113AQBOO+20Xt9ftGgRvva1r+l8aCGc7Go5le2sSVpRpdoqC8ESiA1LEbKdF3EBF6D+5Glyl+SsmZ6oW75M9sHbwfMf1basS9JW5XLctViIK1NcS2lfvrGVS1s5W8HEyVb1TAvZVlSpk3sReORi4N3nzF/1ywRzslScPE1mDrJm01j9sHbpoKHBEDcT/6aTtiqXMj3DxjHjmwY7dU7L7YZ8nOwqyPR29yp3MZXdWE9m+igAbFxpfgdeQHxqbBIqxtHbuCGfLrJLkL+7R+/xlDKxdYPsv6lKXJ1hQ0JyG5S4WgSUKRu2u1e5Qy8gt7+J7Mm9PHgysQNvSCSYSyrNOHrb9uLRTfYk/Nz12b9XsiT1b6rQ+0vuWpwLuQ1KXNrK2QibtrtPs1FaJTL7m6g4uZu46o8M5gT3Q6kmzXq+bXvx6CYb2Nb08yNDFEX039TIU9T8eyenZDpm3jYuFgFlxmRxXqW1bh0zLUTmg0S1VYp05QDmWhir7sUyDdj1UbAzrGwHRNr1fNv24smCzLwQ3e8VGwqxRcfsT5xlx/E6xIdu0lwHJYB7RUCZMFWcJ9IiqmPQWJyKbZUFYPRpwLsrxX6GqRbGasFcz6RcSWnX813ZkE+lkZOBqd8PlmZE6Hiv2NZ+Lbp5oYl/747ypZs090EJVaByPogo21tEmyYBxe6lrPUrgmPb+Gvx+5u+6i//cJcdrKVyjoXPG/JVc/LXgee+D62Tgaux9d9W3ib6auRTN2lua0oogkxxnooPUJvqV6oJi343/GfvD/byYrxKamqDPV1s+8CtVqsz8hRg5J/rW89XXbzsgv51wXsgrrakpp/a94oL/7ZMdAF5xqduUmZKqK+sU+wm6ldkrs4ix2ILXPl27QtqONpetu9EG3W1qvMKVjR97xMTGSLTg9soEz51kzJT0m1N23Y8vvo9rGnbbvpQ7KByPkiUrFtEk8xdEZmrEee935mZWSKq0tWq7ivYYnH/n9KvfZV1hihv7dc55lM3KTMl8KdASCnRCvm0H6BZ1q8kWVtXNRbbt0mlabXcHzwXvXZgLgavw9qngGm3+Nn6mWWGyERtGBnjSzdp7oMSnwqElIv7AK203b2srFpEk442VzEWuxRT5d2vxdXdX5RnRrq/XjbX3+AtqwLPPLZf55wP3aS5D0qiCoRcf3GVqPQBuuX1vtvdJ20vzKp+Jenauup9ZLjtOvCrfxa73cofAjP/3f3OjGrBh+521zy2X5Pzch+U+FQgpFX4AaqjvVB3AWCauSuiH+wyXEqVq76a37MbaF0ldtvWl4KaHxvmaiRhw2yQPLZfk9NyX+jqU4GQdrraC3UXAKYdba56HxkXUuW6NmLcuU3u9jZscJiEDftGAflsv84JX5szcp8pAfwpENJOZ3uhzgLAtGvrsoPGoriQKrd12JYrxcJJ65d0yWP7ted8bs5gUNLNhwIhrbIYPa+rAFDF2nrFUfMJ2J4q131CPWhI6kOsGvjaMhnUxtkgnJ7qDd+bMxiUkJgs2wt1FACqWFsv/WDfuQ34X8dCaHga0L15X5f9qXLdJ9T+dUDTFKBNsK6kkvLA14bajZCpfaNEcS8Z5/nenJH7mhISlPXoedVUrq33rwMGN4mNDQ+NmWb/tutZDds643vJ7lcqDHxN1W7s2R20w5c/B2nrl4hi+N6cwaCkijVt2/EfLe/gzXXrOekQ2L8EkvXeHSpV2+sl6b4uQgWwBeDiXwIzFtudIQGyO6GOnBwMR0ujUAN88Pvs93WJKwBWEbxXC3iI4H9zBpdvKvjZo/+GQ998AF+seRW1hSK6UIOaMY61I+rgQ3uhyrV10am3o09Vc+y6ZTlsq7RG5+2nECyDFYD6YcHSWDEi2Atrf165P9vaDdEC4KT1SzYtQ5HVfG7OKBSL9m420dHRgYaGBrS3t2PQoEGZPOamZ36M/2vVd9GFGvQr7P9wLhZqUQhrAmxOwevW8kD8Sdj356c8oGl7uW9nw5jpbnY2LJkpfkKdsVjNY5YP5lt0NqJrdQrAxU8Ci88VD6Cu25Iue9faLHZcl64I9u8RvW34/igNeLL4d8WCV9Is6fmbmZJSrc0Y0fw9FApADXp/2BXCKzfb2xF109Ve6MKHZNSV7IzFbvwOcUxkw0qLL0WyT5OvCCYKZ7mvS/Md3VmkqAxOd1ZmxmK5faOybCFmNoYsx0xJqSUzUVz39P4ApBLVV4kuU3ESduVDMusrWZNsyIZVyz59ukksQCiVNlOyZzdww2EQ6rQqfSzRDFpW2ak8vYfJOGZK0uruPCjY2spno7TthbYO6Sr37vP7N5CzYRiWbjYM26q251K4LCIakKgYVvfuryHc+l2alRGpX8qqhdi2gW4kbE3bdi9rR6phUBLiNt/ZcuFDMszirH0q/ra+7f5ry7Ct0sC3+Q6gUAhqNkR17QMm/n263axXS2QnKhUARwXvWX3u2DjQjWL5PLm1GgYlIW7znS3bPyR7sjiC7Z2+ZtBsGba1Z3d3rYvMsLp9wLCj9xfEJlka3LMb2LBC/DiPOlvu+cric0dFNsZ0cJpDvk9urYZzSkKuzOHwYYZBVkO6kuqVxZHYiI/DsPTp3AHhgATYv7zzwe/TDVWTyWQAwPi/E78tkM3nTpr5M7o2ZqRYUZNbfcagpJTIMCxTczh8+nCwfeplmMWRxQyaPqKTc8uV157IDlWTGYaGAjD6NMkDhP7PnaQD3WzZ6TinfJ/cWg2DklIRo8iLhVoUUUDblB9mv5zg24eDzSPrRbM45Uxn0HyXZmfmSsKlwTiimYxCLTD2i8lef5VbIFSSJBsTW/OlYVou9eL75NZqGJSUqzCKvAs1WL53Ar7c+T18/lejsWD529kdj48fDjYvlcmm60O2T7J1TfkypVTGQoDM0qBIJqPYle71V70FQjnZbIxItlA0sKPE5p09Fk9cPgW3XHgcnrh8Cr7jeZErwELXyko6D97cuAlfXvS/0YkBPX+dabGR7QWhSdk6sl6m8BBwZ/dfWaYKG6Pm1hw1TaxgU5RoN4vodgJpX3+dHU8yv4PtOx3nzPimwd5nR0oxUxKlfx027PpUr4AklEmxke0FoWnoTlknJZrFCbmw+68Mk7VLccuUgz8j1w4cR2ZpUHcmo1T/OuCgoepP9KK/g+01X+Q1ZkpiGC028n12ig1DuioRyeKE+6+Em+350DJpcpidyNya5tuDEfPNt8dnD+MkGapmy+yWNER+B45HIIMYlMQIi41K+8UzKzbKw4eDjR/0Mrv/ujImP47pYXaiy5R/ag2u6ksD2STSLA3aMrsljajfIcwWJtnpmBLL2+TWarj3jSBjbxgTu7ZSIG7vEp/2EjH5PtuzO1gmkt3xd8/uYAT8kr8RfywXXxsTZHZFVhGk2nRRYoCPk1u5941mxoqNbC0IzYOoLI7pzIJKpgsbky5T9q8DxpwDTLulclYrnOqKAoCiHUuDrsiquNeXTGMKeZ3cWg2DEttl9eFA1VVKdfvUFWWydqm1Gfjtj8VvX2mZslpt0phpQQBy2HFmrsJdv/rXXfPlyoacmkVNbmVQQlIyW9KxtSA0r0xnFlQzVbtUelISEVXDEFebZEtbs2v/VnXVfPmUaUwpr5Nbq2FQklDma4A2FoTmlW9dUSYKG6NOStWILFOaLkL19epf9fPqU6YxJaPNFBZiUJKA0TVA0x+65GdXVNa1SyInpZAry5S8+hfjW6ZRgXlnj8WZ44az+wYcnpZIXndvpG42j8lPKsthdrL7C/3ZmW4MqONodjEczlbR+KbBOH/CiFwHJACDkkS4BuiY8n1UVLB5R+mksppaKru/0PRb7c8s+Dx9WTWbN+Qk47h8kwDXAB2hs+DQ166oLGqXfFz+8q3OSKecDmfjcDQxDEoS4hqg5bIoOPS5K0pn7ZKPJyUfAy2dcjZ/ycfhaLowKEkhb7s3OiPLgkN2RSXj20nJx0BLJ18zjRVwOJoc1pSQf0wUHOra2dVXtu4SnYaPdUYqldd2ZbnzskFsjJDDTAn5he2G7vBt+StHV/9S4mq7PM80sjFCDoMSDVjQZBALDuWZPCH4dlLyLdBKS7S2y+P5S2yMkMOgRDEWNBnGgkNxNo1D9+mk5FuglRSHyfVgY4Q41pQoVK2gaU3bdkNHlEM+DjZLIm42S8v9wdb065f3vYJ98Cyg5YHsjtVXea8z4jC5XjgcTQyDEoWMFDSlGQymY6iYLjLHmqeCw/LnpbUZWDITuPFwYOFng/8umQm0vbz/PrFXsMXgCrb0PkQyOEyOEuLyjUKZFjSlSb3blLaPk+RY81BwWOl5GToO2Ppm8DtGrd9zMzTSLYe1XawlVKNQLBaLpg+imo6ODjQ0NKC9vR2DBg0yfThCymtKLjt1NL6juqaktHis2gm3WjtdmvtmLe2xtr3ct+BwzHT3Cw6rPS+xCsDFTwKLzxWvublui/MnCzJgz+4gS5eT9xlrCftKev5mUKKB1oi5tTmoBUDUy1YI+vzLT7xp7ps1lcfqU8Gh0PNSRU0/4IgzgA3PiN/nmg1BXQSRrCUzxYfJzVic3XEptqZtO867c1Wf7z9x+ZRcZ0ySnr9ZU6KB1oKmNMVjLhWeqTxWnwoORZ6Xarr2Au88C6Agdvu03Uku1SyRejmp7eJwNLUYlLgkTfGYS4VnLh1rlkSflyjFLuDIv9TbnSRSbEtiXA7sfJzaWwGHo6nFoCRja9q24/HV7yVrE05SPKbivllz6VizJPO8VFXovoKNC/gSXsGy1VgNXwK7HIySD4ejleJwtOTYfZOh1MVQaQaDuTRUzKVjzZLM81JVsXsJJ0aSK1gOy1Ijix2us5SDYXIcjqYOMyUZUTJYLc1gMJeGirl0rFkSfV4i1QDNt8ffbNg4+R/tUs2SrXyeIeNTbVcFHI6mBoOSjCgrhkpTPOZS4ZlLx5olkeclkkCWpaaffODAOiA1GNhRzjEoyYiyYqg0xWMuFZ65dKxZinpeVEkSOLAOKD0GdtZKVQtIUlhTkhGlO0Wm2YnUpV1MXTrWLFV8XgqAypFDslM2WQeUXg6noFqrpP5lwXN/4GC0DDEoyZDSYqg0xWMuFZ65dKxZKn1efnc38Nz3EcwfURSYyAYOYb2L6LAsvoZ9MbAzr2z7hmKhBsfvPQEnFM7Ba8WjAAS1gGeOG87aEU24fJMx5cVQaYrHXCo8c+lYs7TldeC5H3R/IRiQFOJqFhIWELMOKB0WeJtVoZ29UOzC1JrX8IsBP8BFtc/13JSD0fRhUGIZrl2SFKkJrwVg8hXxV+Jde4Hj/1a+ZoF1QOkxsDMjouupX6ELNQXgh/0exAmFdQA4GE0n7UHJHXfcgc985jM44IADcPLJJ+OVV17R/ZDOWrD8bZx35yrMffQNnHfnKixY/rbpQyKbyU54vfhJ4MwbqgcOPRmUAvBvX0k2sCsHw7K0YmBnhkBw34UazOr3NAejaaZ1Q75HHnkEF198Me6++26cfPLJuPXWW/GLX/wC69atw9Ch8Zt8ubohXxLc1Imk7fwgmPYpqnRzvfIdlMN6lEItUCy5Uk+zezTrgJLzdYdrG0nsaFws1KDg+I7GWUl6/tZa6HrLLbdg9uzZuOSSSwAAd999N5YtW4YHH3wQ8+bN0/nQzomaY8KghCpKUxhZWij77q+BJX8bfL9YtnSQZhJr/zp+eCfFAu/sSHQ9Fdj1pJ225ZtPPvkEr732GqZOnbr/wWpqMHXqVDQ3N1e8T2dnJzo6Onr9yQtu6kTSVBRG9q8DXn+YA7tsxQJvrda0bcd/vN2BYkHwVMiuJ+20BSUffvgh9u3bh2HDhvX6/rBhw/D+++9XvM/8+fPR0NDQ86exsVHX4VmHmzpRImkLIzmwy28u7zKsWVjD983H1uGZvSdgHzR1pZEUq+aUXHvttZg7d27P1x0dHbkKTLipE0kLCyOXXR1kM0qDi9J6kGrLLr4O7Mr7skfZvI2g2HgaMOUK1qSg715k9+89G1+oaQlKq6ph11MmtAUlhx56KGpra7F169Ze39+6dSuGDx9e8T4DBw7EwIEDdR2SE8Y3DWYwQnLSTL71bWAXT8b+7TKsQXkN36vFMfju3kvxo/6LUEgS3JMy2pZvBgwYgBNOOAErV67s+V5XVxdWrlyJyZMn63pYygOmpPtqmgTMWAxctyXosrluS/B13IeoTwO7Kgy/6jkZP3gW0PKA2ePLgs+7DCtUqVbv5/umYsO0X7Cd3TCtyzdz587FV7/6VZx44ok46aSTcOutt2LXrl093Tgkb03b9vwu7/AqOF6SjpfJc4Ir6Ci2p65jT8ZI1kHkmnDeRuSo/+6iZZ+fhxjV9iL7s4ljgYl/yeU/g7QGJTNmzMC2bdvwve99D++//z6OP/54rFixok/xK4lZsPzt/G4MxZS0PmnrUmyQ15Nx6ckT2B+wRyktWs7xCTeyho/t7MZoHZ6WVp6Gp8XJ9XC11uYgLR+5t0shSLH6dMLJmqsDuySGX6FQEyxtuX7CqZQ1PGIq8M5/iv+M0mF6RIpZOTyN1Mn1cLW8XgVnzZWBXeXH52sHUTXVsoYbV0bfr5QLRcsp5HqZ23EMShyR2+Fq4RwNpqSzY2vqulpN0Umz/eogihJZOxMzryZU0y8o3rTxNVYg18vcHuAuwY7I7XC1JFfB5J+ozpqf/jUwdJwfHURxpHaFrsL2ouUUyuePAMDdL2zkrusOYabEIbkcrubbHA2SJ9JZs/XN+J/j+slYNGsYitpc0dNlzlwvc3uCmRLHjG8ajPMnjIj8B7ambTseX/2eH1cHPs3RoGREsgM1/YDhxwIo9H2v1PQLvu/6yVgmawgAR56Ru3kbuV3m9ggzJZ7xcj3VhzkacWwvLjVFpqZo61vAxU8GSz2yk21dIJs1vPCnwf/n6H1Vbf4IsyTuYFDikWrrqWeOG+72P0of5mhUw4Fw0WRrioaODSbZ+hjkhVnD9ctjOtHKCll9+f0F5XKZ2yNcvvFI1Hqq8ybOClLPPo2ATjMW3YdR+yK/Q5gdEFFaU9S/LpjB4dsJOe2u0DkhssxNdmKmxCPer6e6MkdDRNKx6D5kVmR+h6TZAV/5nDWMwdkj+cBMiUdy0zbsw1WwUPFm90C4kA8bziX5HZgd6M3HrGGMBcvfxnl3rsLcR9/AeXeuwoLlb5s+JNKEY+Y9xCsKyyUZi77ldfdH7afZLqDlgfjsgIcn41g+ZA1j5HqLDYclPX8zU+Ih0fVUr1qHXZJkIFySzIpt0vwOOcwOCPEhaxjD61o56oM1JTnlZeuwK2RbO2v6uT9qX8V2AT7VFJEw72vlqBdmSnKIo5gNkx0I17XX/VH7KrcLyEF2gPbLTa0cAWCmJJc4itkCMgPhfBi178PvQMZw9kh+MFOSQ0yHWiBs7RQZi+7DqH0ffgdSSramjbNH8oFBSQ4xHWoJmeJNH9piffgdSAm2+FI1bAnOMbYOW0SkeNOHtlgffgdKhS2++ZD0/M2akhwb3zRY6EOAwUsG+tfFL1lMnBVMeG2+w90N53z4HSgV1rRRFAYlFImtw5bxoS3Wh9+BEmNNG0VhTQlVxdZhi/nQFuvD70DSWNNGUZgpoaqYZqXcYzZHmMwyL1t8qRoGJVRVbtKsH/8R2PE+UD8c+NTBpo9GDZ5M03F9N+aMX/8ky7yiNW2ULwxKqKowzVr6YeNVmvWV+4AXFwI739//vYOGA6f+P+52gLh+MrVBy/3AsmuCDqHynYzXLrW7Q8jA619tmffMccP9+aygzLAlmGLJpGWd6dT590uBtx6r/vdHfxn48gPZHY8KpSdTttsmk2YnY9MMvf6Pr34Pcx99o8/3b7nwOJw/YYTyxyM3sCWYtBFNszrTqfPKfdEBCQC89e/AyCnunMRbm4MTEoq9T0jA/q+XXR2049p2MrVJuJNx+XNYKtzJ2Kbn0eDrn5tlXsoEu29ICac6dV5cKHi7m/Ueh0rhyTRKeDKlysKdjKMCEqD3Tsa2MPj6s5uGVGKmhJRwplPn4z/2riGJsuN/gtvbXvwankzjNrsrPZmy+LWvJDsZ2/A8anr92U1DJjAoISWcSeHuEAxISm9ve1Di6snUNq7uZKzh9Wc3DZnC5RtSwpkUbv1wvbc3ITyZirDpZGobV3cyVvz6O7UUS95hpoSUSZLCzbxb51MHB22/Iks49YfpzZKomiURnkzXL48p0OwX7DFjy8nURpPnBEscUWzbyVjx6+/MUix5iUEJKSWTwjXWrfP5a4CnrxG43bf1PL6OWRIunkxtNHJy0Dobt5OxTZ03gNLX35mlWPISl2/ICKMp4pNmB3NIohz9ZT3twC33B3Mw1i/vO5jrwbOAloSzUcKTKQp9lx9q+gXft/FkaqOJs4I5JEeds39ZJNzJ+NIVdraJK3z9nVmKJS9xeBoZkXrgkoqlj5YHgrbfHf+z/3v1hwUZEh0nniwGc7W9HGRh1i7dn4UZMz24QmZAIs+1cf1VXv/1oy/GW7VjpZZJnRmESFbi8DRySuIUscqlj4mzgj9Z7X2TxWCupknBH9dOprbqX+fW81fh9V/w3B9w92MbAQQXAaLLpOymIRO4fENGJEoRt9yP4qKz0bXuabVLH586GBj2Of1FrVkO5upfBxw01K0TKqnT/fqv+Z//w04acgozJWSMVLdOazOKy65GAUChuK/337kwRp2zRMgAdtKQaxiUkFGiKeLtv/pX1Bdr0K8QcWK3cU+SkKuDucg6MrUe7KQh1zAoIfvt2Y1Ptz6LQlRAAtg9Rp2zREgB2Tb6cJm09D7spCGbMSgh+3XuQAHJlj6s6iDgLBFKoVob/Znjhke+t7kvDbmEQQnZL+HSh7HhbNW4OpiLrJCmPoSdNOQKdt+Q/RLsSWLt/h0uDuYiLda0bcfjq98Tfk+yPoTygJkScoPk0ofVXQecJZJ7SXfhZX0I+Y5BCblBcukjzVVlZnUorg3mIiWS1oYArA8h/zEoIXdMnBXMISkfo33UOX3GqCe9qrSuDoW8kzaLx/oQ8hmDEnKLxNKH7FVlmitYyi/ZzBprQ4iqY1BCbhJc+pC5qkxzBWtV6zFlhrUhRGoxKCHqlvQKlks++cTaECL12BJM1C3JJoHWth6TNNkW3ajMmojxTYNx/oQRDEiISjBTQlRC9go2bdEil33skCTbxdoQIvUYlBCVkalDSXNi4rKPHrKBXtJlGNaGEKnHoIQohaQnprSdPnnIsCT5HZMEemmyXawNIVKLQQlRSklOTGlOhGkyLEmDmazvl+R3TBropV2G4dwQInUYlBApIHtiSnoiTJNhSRrMZH2/pL9j0kCPyzBE9mD3DZEBSTp9gOQdH0m7hLK+H5D8d0yT8Zh39lg8cfkU3HLhcXji8in4Dmt7iIxgpoTIkCTLPklPvEmzCFnfD0j+O6bNeHAZhsg8BiVEBsmeCJOeeJOe6LO+H5AuuGDhKZHbCsVisWj6IKrp6OhAQ0MD2tvbMWjQINOHQ2QNFZ0pl506WmiZIuv7hfLQYUTkq6TnbwYlRDniSvcNEbmNQQkRERFZIen5m903REREZAUGJURERGQFBiVERERkBQYlREREZAUtQckf/vAHzJo1C6NGjUJdXR2OOOIIXH/99fjkk090PBwRERF5QMvwtLVr16Krqwv33HMPjjzySLz11luYPXs2du3ahYULF+p4SCIiInJcZi3BN998M+666y5s3Lgx/sbd2BJMRETknqTn78zGzLe3t+Pggw+OvE1nZyc6Ozt7vu7o6NB9WERERGSJTApd33nnHdx22234+te/Hnm7+fPno6GhoedPY2NjFodHREREFpAKSubNm4dCoRD5Z+3atb3us3nzZpx11lm44IILMHv27Miff+2116K9vb3nz6ZNm+R/IyIiInKSVE3Jtm3b8NFHH0XeZvTo0RgwYAAAYMuWLTjttNMwadIk/OQnP0FNjVxihjUlRERE7smkpmTIkCEYMmSI0G03b96M008/HSeccAIWLVokHZAAQBgvsbaEiIjIHeF5W7aXRkuh6+bNm3Haaadh5MiRWLhwIbZt29bzd8OHDxf+OTt27AAA1pYQERE56KOPPkJDQ4Pw7bUEJc8++yzeeecdvPPOOxgxYkSvv5OJmg4//HBs2rQJ9fX1KBQKyo6vo6MDjY2N2LRpE5eFMsbn3gw+7+bwuTeHz7057e3taGpqiu26LZfZnBKbsFbFHD73ZvB5N4fPvTl87s1J+txz7xsiIiKyAoMSIiIiskIug5KBAwfi+uuvx8CBA00fSu7wuTeDz7s5fO7N4XNvTtLnPpc1JURERGSfXGZKiIiIyD4MSoiIiMgKDEqIiIjICgxKiIiIyAq5Dkr+8Ic/YNasWRg1ahTq6upwxBFH4Prrr8cnn3xi+tC8dMcdd+Azn/kMDjjgAJx88sl45ZVXTB+S9+bPn4+JEyeivr4eQ4cOxbnnnot169aZPqzcWbBgAQqFAr75zW+aPpRc2Lx5M2bOnIlDDjkEdXV1OOaYY/Dqq6+aPizv7du3D9/97nd7nVN/+MMfSk1y1zJm3hVr165FV1cX7rnnHhx55JF46623MHv2bOzatQsLFy40fXheeeSRRzB37lzcfffdOPnkk3HrrbfizDPPxLp16zB06FDTh+etF154AXPmzMHEiROxd+9eXHfddfjCF76A3//+9zjwwANNH14utLS04J577sGxxx5r+lByYfv27TjllFNw+umnY/ny5RgyZAg2bNiAwYMHmz407910002466678NBDD2HcuHF49dVXcckll6ChoQFXXnml0M9gS3CZm2++GXfddRc2btxo+lC8cvLJJ2PixIm4/fbbAQBdXV1obGzEN77xDcybN8/w0eXHtm3bMHToULzwwgv4/Oc/b/pwvLdz505MmDABd955J370ox/h+OOPx6233mr6sLw2b948/Pa3v8VvfvMb04eSO9OnT8ewYcPwwAMP9HzvS1/6Eurq6vCzn/1M6Gfkevmmkvb2dukNhCjaJ598gtdeew1Tp07t+V5NTQ2mTp2K5uZmg0eWP+3t7QDA93hG5syZg2nTpvV675Nev/zlL3HiiSfiggsuwNChQzF+/Hjcd999pg8rF6ZMmYKVK1di/fr1AIA33ngDL730Es4++2zhn5Hr5Zty77zzDm677TYu3Sj24YcfYt++fRg2bFiv7w8bNgxr1641dFT509XVhW9+85s45ZRTcPTRR5s+HO8tWbIEq1evRktLi+lDyZWNGzfirrvuwty5c3HdddehpaUFV155JQYMGICvfvWrpg/Pa/PmzUNHRwfGjBmD2tpa7Nu3DzfccAMuuugi4Z/hZaZk3rx5KBQKkX/KT4abN2/GWWedhQsuuACzZ882dORE+syZMwdvvfUWlixZYvpQvLdp0yZcddVV+PnPf44DDjjA9OHkSldXFyZMmIAbb7wR48ePxz/8wz9g9uzZuPvuu00fmvceffRR/PznP8fDDz+M1atX46GHHsLChQvx0EMPCf8MLzMlV199Nb72ta9F3mb06NE9/79lyxacfvrpmDJlCu69917NR5c/hx56KGpra7F169Ze39+6dSuGDx9u6Kjy5YorrsDSpUvx4osvYsSIEaYPx3uvvfYaPvjgA0yYMKHne/v27cOLL76I22+/HZ2dnaitrTV4hP467LDD8LnPfa7X98aOHYvHHnvM0BHlx7e//W3MmzcPX/nKVwAAxxxzDFpbWzF//nzhLJWXQcmQIUMwZMgQodtu3rwZp59+Ok444QQsWrQINTVeJo+MGjBgAE444QSsXLkS5557LoDgamblypW44oorzB6c54rFIr7xjW/giSeewPPPP49Ro0aZPqRcOOOMM/Dmm2/2+t4ll1yCMWPG4Dvf+Q4DEo1OOeWUPm3v69evx8iRIw0dUX58/PHHfc6htbW16OrqEv4ZXgYlojZv3ozTTjsNI0eOxMKFC7Ft27aev+MVvFpz587FV7/6VZx44ok46aSTcOutt2LXrl245JJLTB+a1+bMmYOHH34YTz75JOrr6/H+++8DABoaGlBXV2f46PxVX1/fp27nwAMPxCGHHMJ6Hs2+9a1vYcqUKbjxxhtx4YUX4pVXXsG9997LLHgGvvjFL+KGG25AU1MTxo0bhzVr1uCWW27BpZdeKv5Dijm2aNGiIoCKf0i92267rdjU1FQcMGBA8aSTTiq+/PLLpg/Je9Xe34sWLTJ9aLlz6qmnFq+66irTh5ELTz31VPHoo48uDhw4sDhmzJjivffea/qQcqGjo6N41VVXFZuamooHHHBAcfTo0cV/+qd/KnZ2dgr/DM4pISIiIiuwgIKIiIiswKCEiIiIrMCghIiIiKzAoISIiIiswKCEiIiIrMCghIiIiKzAoISIiIiswKCEiIiIrMCghIiIiKzAoISIiIiswKCEiIiIrMCghIiIiKzw/wMbjPUpPpkIegAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "m = 200\n",
    "Y = np.random.randn(2, m) + np.random.rand(2, m) * 5\n",
    "X = np.vstack((Y, np.array([1,] * m)))  # 中心を計算する用\n",
    "n = 2 + 1\n",
    "\n",
    "# minimum volume ellipsoidの計算\n",
    "# 双対問題を解きます\n",
    "def dual_fn(u):\n",
    "    x = X @ np.diag(u) @ X.T\n",
    "    return -np.log(np.linalg.det(x))\n",
    "\n",
    "u0 = np.ones(m) / m\n",
    "cons = (\n",
    "    {\"type\": \"eq\", \"fun\": lambda u: u.sum() - 1}\n",
    ")\n",
    "bnds = [(0, None)] * m  # u >= 0\n",
    "res = minimize(dual_fn, x0=u0, constraints=cons, bounds=bnds, method=\"SLSQP\")\n",
    "u = res.x\n",
    "H = np.linalg.inv(X @ np.diag(u) @ X.T)\n",
    "H = H[:-1, :-1]\n",
    "L = np.linalg.cholesky(H)\n",
    "y_ = Y @ u\n",
    "\n",
    "# 楕円の描画\n",
    "theta = np.linspace(0, 2 * np.pi, 100)\n",
    "x = np.cos(theta)\n",
    "y = np.sin(theta)\n",
    "L_T = np.linalg.inv(L).T\n",
    "x2, y2 = [], []\n",
    "for _x, _y in zip(x, y):\n",
    "    v = y_ + np.sqrt(2) * np.dot(L_T, np.array([_x, _y]))\n",
    "    x2.append(v[0])\n",
    "    y2.append(v[1])\n",
    "plt.scatter(x2, y2, s=10)\n",
    "\n",
    "# データの描画\n",
    "x, y = Y\n",
    "plt.scatter(x, y, s=50)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "いい感じですね！\n",
    "今回は最小楕円問題をscipy.optimizeで解きましたが、もっと効率がいい方法を見てみましょう。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 最適解の性質\n",
    "\n",
    "主問題と双対問題の関係を踏まえると、最適解が満たすべき性質が見えてきます。\n",
    "まずは主問題と双対問題の関係を確認してみましょう。\n",
    "\n",
    "---\n",
    "\n",
    "**主問題と双対問題の関係**\n",
    "\n",
    "$H$と$u$がそれぞれ主問題と双対問題について実行可能であるとき、\n",
    "\n",
    "$$f(H)\\geq g(u)$$\n",
    "\n",
    "が成立します。証明しましょう。\n",
    "まず、変形して\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "f(H) - g(u) \n",
    "&= -\\ln \\det(H) - \\ln\\det(XUX^T)\\\\\n",
    "&= -\\ln \\det(HXUX^T)\\\\\n",
    "&= -\\ln (\\prod_j \\lambda_j)\\\\\n",
    "&= -n\\ln \\left(\\prod_j \\lambda_j\\right)^{1/n}\\\\\n",
    "&\\geq -n\\ln \\left(\\sum_j \\lambda_j / n\\right)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "です。ここで、$\\lambda_j$は$HXUX^T$の固有値であり、最後の不等式では算術平均と幾何平均の関係を使いました。この固有値について、\n",
    "\n",
    "$$\n",
    "\\operatorname{Trace}\\left[HX\\left(\\sum_i u_i e_i e_i^T X^T\\right)\\right]\n",
    "= \\sum_i u_i x_i^T H x_i\\leq n e^T u=n\n",
    "$$\n",
    "\n",
    "が成り立ちます。上の不等式ではラグランジアンが最小化されるとき、$H=(XUX^T)^{-1}$であることを使っています。\n",
    "\n",
    "このことから、$HXUX^T$の固有値の和は最大でも$n$であることがわかります。よって、\n",
    "$$\n",
    "f(H) - g(u) \n",
    "\\geq -n\\ln \\left(\\sum_j \\lambda_j / n\\right)\n",
    "\\geq -n \\ln (n / n) = 0\n",
    "$$\n",
    "が成立します。\n",
    "\n",
    "証明は省略しますが、$H$と$u$はそれぞれ唯一の最適解を持ち、それは$f(H)=g(u)$を満たします。\n",
    "\n",
    "---\n",
    "\n",
    "続いて最適解が満たす条件を確認します。\n",
    "\n",
    "---\n",
    "\n",
    "**最適解の必要十分条件**\n",
    "\n",
    "主問題の解$H$と双対問題の解$u$がそれぞれ最適である必要十分条件は次である。\n",
    "\n",
    "1. 実行可能性：$i=1,\\dots, m$について、$e^T u=1\\quad, u \\geq 0, \\quad x_i^T H x_i\\leq n$\n",
    "2. $H=(XUX^T)^{-1}$\n",
    "3. $i=1,\\dots, m$について$u_i>0$ならば$x_i^T H x_i=n$\n",
    "\n",
    "これは最適解が$f(H)=g(u)$を満たすことや、主問題と双対問題の関係の不等式を確認すれば出てきます。\n",
    "\n",
    "---\n",
    "\n",
    "最後に近似最適な解について述べておきます。\n",
    "\n",
    "---\n",
    "\n",
    "**近似最適解**\n",
    "\n",
    "* 実行可能解$u$は次を満たすとき$\\varepsilon$-主実行可能であると言う。\n",
    "    * $H(u):=(XU^TX)^{-1}$が $x_i^T H(u)x_i \\leq (1 + \\varepsilon) n, \\quad i=1,\\dots, m$ を満たす。\n",
    "* 実行可能解$u$が$\\varepsilon$-主実行可能かつ次を満たすとき$\\varepsilon$-最適であると言う。\n",
    "    * $x_i^T H(u)x_i \\geq (1 - \\varepsilon) n, \\quad \\text{ if } u_i > 0, \\quad i=1,\\dots, m$ を満たす。\n",
    "    * つまり、$(1 - \\varepsilon) n \\leq x_i^T H(u)x_i \\leq (1 + \\varepsilon) n$ です。\n",
    "* $u$が実行可能で、$u$の目的関数が最適なものから$\\delta$以内であるとき、その$u$は$\\delta$-以内最適解であると言う。\n",
    "\n",
    "以上の定義を踏まえて、近似最適解について以下のことが言えます。\n",
    "\n",
    "\n",
    "**近似最適解の性質**\n",
    "\n",
    "* $u$が$\\varepsilon$-主実行可能解（もしくは$\\varepsilon$-最適解）であるとき、$u$は双対問題について$n\\ln (1+\\varepsilon)$-以内最適な解であり、$(1+\\varepsilon)^{-1}H(u)$は主問題について$n\\ln (1+\\varepsilon)$-以内最適な解である。\n",
    "* また、$\\mathcal{E}((1+\\varepsilon)^{-1}H(u))$は全ての$x_i$を含み、かつ最小楕円の$(1+\\varepsilon)^{n/2}$倍以内の楕円である。\n",
    "\n",
    "証明してみましょう。まず主実行可能解の定義より、$u$と$(1+\\varepsilon)^{-1}H(u)$は実行可能解です。また、$H(u)$の定義より、\n",
    "\n",
    "$$\n",
    "f((1+\\varepsilon)^{-1}H(u)) - g(u) = n\\ln (1+\\varepsilon) + f(H(u)) - g(u) = n\\ln (1+\\varepsilon) \\leq n\\varepsilon\n",
    "$$\n",
    "\n",
    "になります。最適解は$f(H^*) - g(u^*)=0$であることから、１つ目が成立します。\n",
    "２つ目は主実行可能の定義と[Determinants and Volumes](https://textbooks.math.gatech.edu/ila/determinants-volumes.html#:~:text=The%20theorem%20on%20determinants%20and,a%20number%20A%20a%20B%20.)から直ちに成立します。\n",
    "\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 最小二乗法との関連\n",
    "\n",
    "MVEE問題は最小二乗法を考えると自然に出てきます。\n",
    "例として、$i=1, \\dots, m$について\n",
    "\n",
    "$$y_i = x_i^T \\theta$$\n",
    "\n",
    "なる線形な関係を持つ入力$x_1, x_2, \\dots, x_m \\in \\mathcal{R}^d$と出力$y_1, y_2, \\dots, y_m \\in \\mathcal{R}$を考えましょう。\n",
    "また、実際の観測値として、それにノイズが乗った\n",
    "$$\\hat{y}_i = x_i^T \\theta + \\epsilon$$\n",
    "を考えます。\n",
    "\n",
    "また、行列形式でまとめるため$X = [x_1, \\dots, x_m]$, $Y=[\\hat{y}_1, \\dots, \\hat{y}_m]=x^T\\theta + \\varepsilon$とします。ここで$\\varepsilon$は$m$次元のノイズであり、$N(0, \\sigma^2I)$に従うとします。\n",
    "\n",
    "このとき、最小二乗法の解は\n",
    "\n",
    "$$\\hat{\\theta}:=(XX^T)^{-1}XY$$\n",
    "\n",
    "になります。\n",
    "書き直すと、\n",
    "$$\\hat{\\theta}=(XX^T)^{-1}X(X^T\\theta + \\varepsilon) = \\theta + (XX^T)^{-1}X\\varepsilon$$\n",
    "なので、$\\hat{\\theta}$は普遍推定になっています。\n",
    "ここで、推定値の共分散行列を見てみましょう。\n",
    "\n",
    "$$\\mathbb{E}[(\\hat{\\theta} - \\theta)(\\hat{\\theta} - \\theta)^T] = (XX^T)^{-1}X\\mathbb{E}[\\varepsilon \\varepsilon^T]X^T(XX^T)^{-1} = \\sigma^2 (XX^T)^{-1}$$\n",
    "\n",
    "推定値の共分散行列はデータ点$X$によって変化するわけですね。\n",
    "ここで、それぞれのデータ点$x_1, \\dots, x_m$について重み$u_1, \\dots, u_m$を割り当てることを考えましょう。\n",
    "このとき、$U = \\operatorname{Diag}(u_1, \\dots, u_m)$とすると、推定値は\n",
    "\n",
    "$$\\hat{\\theta}:=(XUX^T)^{-1}XUY$$\n",
    "\n",
    "になり、その共分散行列は\n",
    "$$\\sigma^2 (XUX^T)^{-1}$$\n",
    "になります。この共分散行列が「小さく」なるように$U$を選択できないでしょうか？選択肢の一つは、「共分散行列によって与えられる体積が小さくなるような$U$の設計」であり、これは**D-optimal design**と呼ばれます。\n",
    "\n",
    "つまり、次の問題を解きます\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\max_{u\\in \\mathcal{R}^m} \\ln \\det (XUX^T)\\\\\n",
    "&e^T u = 1 \\\\\n",
    "&u \\geq 0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "また別の選択肢として、「一番大きな分散が小さくなるような$U$の設計」が考えられます。これは**G-optimal design**と呼ばれます。\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\min_{u\\in \\mathcal{R}^m}\\max_i x_i^T(XUX^T)x_i\\\\\n",
    "&e^T u = 1 \\\\\n",
    "&u \\geq 0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "実はD-optimalな$u$はG-optimalでもあります。（Kiefer-Wolfowitzの定理と呼ばれます。[BANDIT_Kiefer_Wolfowitz.ipynb](BANDIT_Kiefer_Wolfowitz.ipynb)参照）\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "shumi-VTLwuKSy-py3.9",
   "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"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "d6b7cac5e0d2ff733f340da4d53ae5ecfef7f7ad39623f5982b029a09306b36b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}