Appendix_D_sensitivity_of_Algorithm_1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import tensorly as tl\n",
    "import time\n",
    "import os\n",
    "from tensorly.decomposition import tucker, parafac, non_negative_tucker\n",
    "from tensorly import tucker_to_tensor\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn.utils.extmath import randomized_svd\n",
    "from itertools import combinations\n",
    "import multiprocessing as mp\n",
    "from collections import defaultdict\n",
    "from copy import deepcopy\n",
    "\n",
    "from util import std_logistic_function, get_square_set, square_unfolding, normalized_error, tenips_general, \\\n",
    "tenips_general_paper1, generate_orthogonal_mats, one_bit_MC_fully_observed, one_bit_TC_fully_observed_gd\n",
    "\n",
    "import apgpy\n",
    "from scipy.sparse.linalg import svds\n",
    "import mkl\n",
    "mkl.set_num_threads(1)\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "from matplotlib import rc\n",
    "rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'weight': 'bold'})\n",
    "rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "generate data tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "get core and factors\n",
      "get real_X\n"
     ]
    }
   ],
   "source": [
    "# Low rank + noise\n",
    "r = [5, 5, 5, 5]\n",
    "I = 100\n",
    "N = len(r)\n",
    "noise_para = 0.1 # noise parameter\n",
    "\n",
    "np.random.seed(1)\n",
    "# real_core = tl.tensor(np.random.uniform(-1, 1, r))\n",
    "real_core = 100 * tl.tensor(np.random.randn(*r))\n",
    "factor_matrices = [generate_orthogonal_mats((I, r[0])) for i in range(N)]\n",
    "print(\"get core and factors\")\n",
    "real_X = tucker_to_tensor((real_core, factor_matrices))\n",
    "print(\"get real_X\")\n",
    "fro_norm = np.linalg.norm(real_X)\n",
    "noise_coeff = noise_para * fro_norm / I ** (len(r)*1./2)\n",
    "noise_ten = tl.tensor(np.random.normal(0, 1, [I]*N))\n",
    "\n",
    "X = real_X + noise_coeff*noise_ten"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "generate propensity tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "get core and factors\n",
      "get real_A\n"
     ]
    }
   ],
   "source": [
    "A_noise_para = 0.1 # noise parameter\n",
    "np.random.seed(2)\n",
    "A_real_core = 100 * tl.tensor(np.random.randn(*r))\n",
    "# A_real_core = np.random.uniform(-1,1,r)\n",
    "A_factor_matrices = [generate_orthogonal_mats((I,r[0])) for i in range(N)]\n",
    "# A_factor_matrices = [np.random.uniform(-1, 1, (I, r[n])) for n in range(N)]\n",
    "print(\"get core and factors\")\n",
    "real_A = tucker_to_tensor((A_real_core, A_factor_matrices))\n",
    "print(\"get real_A\")\n",
    "A_fro_norm = np.linalg.norm(real_A)\n",
    "A_noise_coeff = A_noise_para * A_fro_norm / I ** (len(r)*1./2)\n",
    "A_noise_ten = tl.tensor(np.random.normal(0,1,[I]*N))\n",
    "\n",
    "A = real_A + A_noise_coeff * A_noise_ten\n",
    "A_true = A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = std_logistic_function(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([       0., 10000000., 20000000., 30000000., 40000000., 50000000.,\n",
       "        60000000., 70000000., 80000000.]),\n",
       " <a list of 9 Text yticklabel objects>)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAExCAYAAACnAX83AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dT1ccR7om8OfxnW1bCLW2HrnwftSAPkC3C83i7nzBmg9ggWbWNkjzASyBez1uYL6ALLl3vfAF+e5bgGaWczyUfNeyoNRn9u8sIhKSIv9HVlZm1fM7pw5SVVZGkJXUmxEZ8QbNDCIiIk34aNIVEBGR2aGgIyIijVHQERGRxijoiIhIY/7DpCvQJSRXATwAcAbg1Mx2JlwlEZFOoUavXSK5aGYnsf+vAhgCWDSzHZI9MxuQ7AM4MrPhxCorItJB6l7zfCB5Efv/IgCY2SGAoQ9Ig+h1BRwRkfIUdDwfXAaxpx7AtXLgn+8DF62fs2ZrJyIyHWYm6JBcJNmL/b9Pci7jLXO4Glxu+Z89XA1OIiJS0MwEHX+vpk+y57vSKnWRmdmOutZERKqZmaADAGa2B2ALQM93p2UZApj3/54D8H6cdRMRmQUzFXTigwXiXW0pnsN1pcH/zAtSIiKSY2aCTnw0mm/xXLmn4wcILPufUXdcFKiG8aHUIiJSjebpiIhIY2ampSMiIpM3M2lwfv/739udO3cmXQ0RkU45Pj7+zcxu17W/mQk6d+7cwdHR0aSrISLSKST/vc79qXtNREQao6AjIiKNqRx0SN4deXwce+1rkr+QfE/yOck7dVRWRES6LaSlswLgBMATAMvRkySfAdgG8AouSeY+gF0FHhERCRlIMASwYmavoidI/gHAJoBdM/uvsedfA3gMF6BERGRGhbR0bsQDjvcAgMG1dC6Y2QdoOQARkZkXEnQ+JDwXpYz5NeE1pT4QEZlxIUHnShAheQPAIoC0yTC3Up4XEZEZERJ0RhdA6/ufL0Y39AFJ3WsiIjMuZCDBW5JfmNlf/XDpbQDnZvY/4xuRvAtgIz6wQESuuvP4bxMp99dn/zyRcmV2VW7pmNmPAO6TPIMbydaDb+2QvEHyG5JHcMOq10l+VUeFRUSku4Jyr5nZI5KPAXxqZm9GXj7xj0jp7jWSp3ABrcj7d83sZdkyRESkOcEJP81sSNJ8S+bMzP7qh0i/irrfAnaft7pn3EZAOSIi0oCg3Gs+/c1rAMcA9nD9i/8Nyb/4+zpl970IYAfATTNj2gMuM8KGmQ1CfhcRERm/kNxrNwD8DOAcwJcAlgAcxrcxs7dm9gjASoU0OMtwXWbDtA38ctMbfvlpERFpuZCWzjMAD83svpn96O/pJE4ANbPvUL77a65A62UbwMOS+xURkQkJykjgR7AVxTI7N7OdzJ2RqwCOs1pCIiLSLiEDCX4ruf2NgLKu8N1qD8xsra59iojI+IW0dMqmtakzDc42gK0a9yciIg0ICTok+cfR51I2/B7A3wPKiu9rEQCKjFYjuU7yiOTRu3fv6iheREQChGQkeAxgn+S3JH8XPR3fxg+p/gnAspn9OaCecdsAdgvWcc/Mls1s+fbt2zUVLyIiVYVODr0P4F8BbJEcAjgj+cC/1oNLCnqCy2SgQUj24ALYSe7GIiLSOkGTQ81sYGafwa0IOgSwADdfZwmuq+2Rmd3zGQrqsIGRuUAiItIdQUEnYmY7ZrZgZh/BBZ6bZjZvZvt17D9mHYAyD4iIdFQtQSfOZyG41rIh+aeQ/foBBHMAXofsR0REJqf2oJOh0M3/DNF9IU0GFRHpqNyBBCS/BnBiZj+PPP8U11cPzVImY3SSFf9TK5CKiHRUZtAh+TlcpudzXJ/cuYZygSQxL1sJUVlq6YiIdFRe99oRgJcAHie8NoDL8PxRkQfCg8V84PtFRGTCMls6fkDAlykvDwAclCjrbYltk5TpyhMRkRaqPDnUr5NTZvvlqmV5hwB6WqxNRKS7gperboqZreRvJSIibRaycui/kPwldP6NiIjMjpB5Ohtw2QdCh0KLiMiMCOleOzCz+7XVREREpl5QRgKSH5fY9qeQskREpPtC1tP5DsBOiXs66oYTEZlxlbvXSH4FN3l0g+Qu3Lo5AwDvEzbXvR8REQm6p7MD4AYul6heyNk+NA2OiIh0XMg9nTO49Dg3C6TAmUdgGhySPZK7JI9JHvjHZsg+RUSkWSFBZwDgRZFVQc1siIA0OCTX4VLu7JrZkpmt+MmiJwo8IiLdEZIGp9Rw6appcHzA2QbwqQ9ecdEaPTtV9i0iIs1qchG30kj24ALLw4SAA9TQbSciIs0ZS+41P39nHsCZmf0jYFfbAIZm9jLpRTO7GbBvERFpWMiQ6adm9iTl5QUAy367BQCWsW3a/ucArALYq1pHERFpl5CWTh9AYiAxszcA3kT/J3kjJ0glidbxOa5eRRERaZOQoMP8TRwz+0AW3jyy5H8OfKtnPfbaLbjcb4dldyoiIpMTEnQKT/b093hSW0Yp4hkMnpjZ1sg+D0iujD4/ss06fLD65JNPShQtIiLjkBt0SH4OYAXXg0yP5LfIb/H04O7NpAaHjPcBwFrKe9cAnJN8nTHQYA/+ntDy8rIyIoiITFhu0DGzVyQHcEFgDa7lYHDB5nGBMgYAHpnZftVKJg2XNrMhyRMA+wASg46IiLRLoe41M3sLl1HgFckduESfp3BdZlnvy81WkCEKNAcZ2xwBWCS5aGYnAWWJiEgDSk8ONbMB/H0SM/uQ9Qis25n/mTX5M3qtUrYDERFpVtWBBAcY/xf9wP88y9gmWkYhL8O1iIi0QKU0OL4lU+R+zoUSi71Fovk5Wevw3PI/k9bwERGRlmky99p2ye2jOTj3MraZG9lWRERaLDj3Gsm7cK2R+YzN5gAsltmvmQ386LSs9/UBDDSIQESkG4KCDskjAH9AsewEVebJbAE4INnzAxjiZffhgt1S4jtFRKR1QhJ+/gXuJv+XuLzpn+YWgJ/KlmFmhySjwLMUzdfxSx68ALChVo6ISHcEtXTKLORG8k3+Voll7ESTQH0OtmiY9OcKOCIi3RISdE5Lbr9WtSCf2FODBUREOq6x0Ws+q4GIiMywkKAz8CPXCiH5dUBZIiIyBSoHHTP7EcAKyS8KvmWjalkiIjIdgpar9v/sk9yHS76ZNoptHtmZBUREZAaEDCTYAHADl3N0VnK213o2IiIzLuSezhncPJmPijyQnS1aRERmQNBAApQbxqzRayIiMy5kIMF9M/u1xPall0Ig2SN5QHLdZyGIP79O8kX8eRERabfghJ/ARdLPZQBnZvbX2PNfxP9fwRxcUs++31/8tQGAldGcbCIi0l6hCT/vAtiHywRNuMXd4kHmjc/R9hcz+18Vi3kJd/8oymQ9APDczF5WrriIiExEyJDpGwB+hhsqHSX97Me38VkIHpH8huSwTHdczIGZ7VWtp4iItEfIQIJnAB76ezs/mtkbpAyLNrPvoMmhIiIzLyTofPBZCYoqsuaOiIhMsZCg81vJ7W8ElCUiIlMgZCDBrTFvf8EPi14d2deBX/JAREQ6IiTokOQfzezf4s+lbPg9gL9XLGcFwMDMdkb2eUzy+ejzI9usA1gHgE8++aRi8SIiUpeQyaGP4Vbz/Jbk76Kn49uQvEvyJwDLZvbnCsUMAOymtGi2AGyT7Ce8FtVxz8yWzWz59u3bFYoXEZE6hS7idh9uuPSQ5HsAGyRf+8d7AMdwc2tSA0MWMxumdaH55wcAdqtVXUREmhYUdMxsYGafAXgCl9BzAcCSfxDAIzO7Z2YfgmuabACgp1Q4IiLdUMty1Wa2Y2YLPpv0AoCbZjZvZvt17D9DlAJHQUdEpANqCTpxZva2rpZNgRZMtFyCgo6ISAfUHnTqQnIXwCnJ1YzN5vzPswaqJCIigVobdOAGIAyRvfhb1MLRfB0RkQ5oc9B5DWApZwJoH8ChmWlVUhGRDmht0PGTPrfSXie57f+51kyNREQkVGuDjrebtDpoLNPAklo5IiLdUcvKoeNiZickH8JlHpjH5cCBEwCfKuCIiHRLbtAh+Scz+7mJyiTxgUVr8YiITIHM7jW/1PQBya8aqo+IiEyxvHs6QwDfIWFIsk/kKSIiUlhe0DEze2xmvya8Vmp9HJKvy2wvIiLTJy/o3CT5x5TXLOX5NHP5m4iIyDTLG0iwDeCY5CmAI1zNDtAj+T+QsnDbiHkoP5qIyMzLDDpm9tbPkdmHW8FzNHA8KlFW2ZaRiIhMmdwh037I8rVZ/ySPzGy5aEF+UTcREZlhTWYkeFv3Dkke171PEREZn5Cgk5oXLUmZVlERJDcBLNa5TxERGa/KQcfMXtVZkTL8faYnkypfRESqqaV7jeTHJL8m+ZzkLyR/IvktyTt17D/BBrSGjohI5wQHHZJfAzgHsAM34GABbqTbY7iVP78NLWOkvFUAz6HVQkVEOicoyzTJZwBW4YZOH+JqIOgB+C8AHpGcM7P/FlKWL28OwD0ze0kWmR4kIiJtUjnokPwDgEUz+yxlkzf+sUXyX2vKVv3EzEoNYBARkfYI6V5bN7P7RTb02wWt8EmyD7eEtYiIdFRI0Cnbv/UhoCwAWDGzl2XeQHKd5BHJo3fv3gUWLyIioUKCzv8tuX3lNDh+Ts5u2feZ2Z6ZLZvZ8u3bt6sWLyIiNQkJOgu11SKDn5MzNLNBE+WJiMj4hASdk6IripJ8CqBq7rUtM9ur+F4REWmRyqPXzGzfj0rrAdhLWuiN5F24iZzLZnavbBkk11GhW01ERNopaJ4O3Ii0V3DDogdw6+2c4XL9nDkAJwC+LLtjPydnQa0cEZHpERR0zOwDgGXfIlkHsBR7eQDgmZl9V3H3fQD9jEzSPeBqpmkzW0rZVkREWiC0pQPAjRIDsAcAJD81s+BlDPzw6NQh0iQPAPQVaEREuqP29XTqCDgiIjKdmlzETUREZlyXg848cDGPR0REOqBTQYfkKsljkobLVUNPSZ6S3J5k3UREJF8tAwmakje4QERE2q1TLR0REek2BR0REWmMgo6IiDRGQUdERBqjoCMiIo2pHHRI/gvJX0j+qc4KiYjI9App6WzALeSmyZkiIlJIyDydAzO7X1tNcvglq6M1eeb8z20zO2yqDiIiEibong7Jj0ts+1PFMuZIvgBwaGZr/rECYBvAgX9NREQ6oHLQ8evk7JS4p1O1G24fwEMzOxkp/xDADoBVv56PiIi0XOXuNZJfATgCsEFyF26F0AGA9wmbV7r341cPXYXrTltJ2OQAwCbcCqZaYVREpOVC7unsALgBgP7/CznbW4Uy5v3P5ZTXz/xPDWYQEemAkHs6ZwAeA7hpZh9lPeCCx7BsAWY2gAtmn6ZsEgWbk5TXRUSkRUJaOgMAL8zsQ96GZjYkWWlFUR940kRdbrtV9i0iIs2qHHTKDpc2s7Quskr8/Z51AC81bFpEpBtqSYND8i7Jr0h+MfL8F2nvCSxvDsArADtmtpax3TrJI5JH7969G0dVRESkhNB5OndJvgZwDDd6bGNkkzck/0Lybkg5vqxFktskDwC8BfDUzLay3mNme2a2bGbLt2/fDq2CiIgECsm9dgPAzwDOAXwJYAnAlW4uM3trZo8ArJC8U72agJmdmNmWma2Y2U0A836Z6n7IfkVEpDkhLZ1ncJM275vZj2b2BinDov1E0tFWUBAz28NlVgJNDhUR6YCQoPPBzH4ssT3zNynHB54hgF2SmqsjItJyIUHnt5Lb3wgoK8uR/7k6pv2LiEhNQoLOrTFvD5IvSJ6TXMzYLJrHcy9jGxERaYGQoEOSfxx9LmXD7wH8vUIZUd61rDk+Ubfa6wr7FxGRBoVkmX4MYJ/ktyR/Fz0d38YPqf4JwLKZ/blCMYcANvy9mzRRQFIqHBGRlgudHHofbrj0kOR7uIzTr/3jPdz8nXkAVYc1byMjkSjJqCWkrAQiIh0QFHTMbGBmnwF4AjeKbAFuvs4SXFfbIzO7VyQ/W8r+DwGc+ns7V0an+fk5+/CLu4X8HiIi0oyQhJ8XzGwHbqkDkPwUwFnVQJOw7z2SPwDYJhktdRAtV/3QzF7WUY6IiIxfLUEnzswqZZPO2ecQNU8uFRGR5tUadEh+DDeabGhmv9a5bxER6b66skw/9QMHzuEGD5ySfE/y2zr2LyIi0yGopeOzR/8Md49lABdwADdirQfgMck1AH0z+/eQskREpPtCu9dewi1p8DRp4IAfcbblt1PGABGRGVc56JB8CmArK+mnX2p6g+Q3JL+uOEFUpBF3Hv9t0lUQmXpBaXCKZpn2Sxt8FlCWiIhMgZCg877k9qcBZYmIyBQICTqJC7ZlmM/fJJlfpvrYrxQaZSjIyjwtIiItFLSIm5+XU9RZ2QJIzpF8AeC5mS2ZWZRmBwCOSe6W3aeIiExOSJbpfQCFlon2gw6qBIh9uMEKFxmkzWzoc63tAVjXUtUiIt2ROnqN5J8A5HVh3SL5C9yQ6CRzcBmm18zsH2Uq5rvPXvsRcEm24ILeLskffKocERFpsawh0/cBbBbcz1bO68ckD83sPxfcHwA8QEbryMyGJE/gAmMf6YFPRERaIqt77RRuYbSbZvZR4OOfSgYcwAWTU79mTpqoFVR5kIKIiDQnK+icwa1VU8sSBRVE3WVZmQyiJQ7SuuBERKRFUrvX/MTPQpM/iyB5p2Tm6YcAdnNWBI2Wqj6qXDEREWlMLVmmC3pRZmM/Si014PiVQ+fgWmMaRCAi0gG1rKdD8g4uu7rS9HJeLysavJC6uJsfTr0OAJ988knNxYuISFmhSxs8RfERbrXxgwv6ADYyhlTDzPbg5vNgeXm5bAYFERGpWUiW6W/gWhnfIT+v2k0AT6uWNVLuHC4nje7VsU8REWlGSEtnBcCnRUe31Zg54BWAPTPbqWl/IiLSkJCBBCclh1On3nspiuQB3MCBvMmoIiLSQiFB57cyG5vZq4Cy4JN7nijgiIh0V2NZpn0ut0pIbgIYJgUcn4layxyIiHRAaJbp/+6HSxexXaUcP1JtIaOFswylwRER6YSgIdNm9pjk9yR7cHna0kaxzSE/Y/U1vgVzz8yy7getoNqyCSIi0rCQIdM3ABziclG1lZy3lJon4wPZKwBHfgBBkh6Anu7ziIh0Q0hLZxvAOYA15CfcvAXgp5L738XlejxZlOxTRKQjQoJOz8zuF92Y5JsyOzezvJaTiIh0TMjotbQurzRrAWWJiMgUCAk6pe7RmNnbgLJERGQKhASdtyTvFt2Y5NcBZYmIyBQImafzI4AVkl8UfEtwGhwREem2kCHTUdboPsl9uNU700aSzaP+9XRERKRjQkavbQC4AYD+/7XO0xERkekTck/nDG4RtY+KPABoSWkRkRkXEnQGcBkJigoaveYTex6H7ENERCYrZCDBfTP7tcT2y1XK8cFmEy5oKZu0iEiHhbR0SiH5Vcnt+75lsw/XqlK6GxGRjmss6KBkJmgzOzSzJTNbM7OXcPeQRESkw0KGTJdpuSzlbyIiItMuZMj0HtwwaOZsF22jIdMiIjMuJOgMATxF+gi2eQALAFYBPEPg6DUREem+kKBzZmbf5WzzCsAeyW+ggQAiIjMvZCBB4fVufHBaDShLRESmQMg8nbLdZXn3fmpHcp3kEcmjd+/eNV28iIiMaHLI9HyDZQEAzGzPzJbNbPn27dtNFy8iIiMaCTokP4aGTYuIzLyQeTrfF9x0Hu5+jparFhGZcaFLG+QZwmUSeGRmfw0oS0REpkDoPJ2vFExERKSooKUNFHBERKSMkCHTlZYqEBGR2dXkkGkREZlxmUHHD3Vuix7gFnWbdEVERKSavJbOxJaHJtkjeUDymOQ5fNABcO6fOyCp1DoiIh2SN3ptgeR/MrP/3UhtYsxsgBL53UREpP2KDJk+IVlm5c4ox5rBtU6i9DcDM3tQsn4iIjJFigSd+2b2quyOST7DZUvlEMDjsvsQEZHpkhd0BmUDDskvAOwDmIObQPrQzH6sWD8REZkieUGn8D0VP9LtBYA+XBfbHoAtM/tQvXoiIjJNMoNO0TVzSH4NYBsu2AwArJnZm/DqiYjINAnJvQaSd+FaNz24gLNVYAlrERGZUZWCju9K2wawDhdsDuFaN+pKExGRVKWDTmygwE00NFDATwJdAXAK4BbcIIVtP5dHREQ6onDQIXkHwC5iAwXM7NF4qnWl3F0A82a2FntuDsAxyQ0zOxx3HUREpB6Fgs7IQIETuNbN2AcK+BbOl2Z2M/68mQ1JbgB4QfJTMxuOuy4iIhIuL+HnXZK/4DLgbJnZcpWAQ/KnCvXbhht6fY1v4ZwBeFJhvyIiMgF5LZ0TuHQ2L+FaN/+oUgjJG3DdcmXeswg3Ku51Tv3WAWxVqZfIrLvz+G8TK/vXZ/88sbJlcop0r731jyck87ZNQgBVskFHQSprsMAAwBzJngYViIi0X17QGQJYxGUSz6rOADwt+Z57/mdWMHnvfy7mbCciIi1QJPdapS61ETsk10u+p8hibdEAgvnMrUREpBXygk7Z1kmWsvddygQSrSY6RSZ5n0FExisv91ptkz4r7GvOv6/IcOhbSU/61lXUwvp/JP9Pyvt/D+C3kvWbhC7Uswt1BFTPOlWqI7fHUJNsXTiWQPvq+R/r3FlQ7rW2M7M9pAy5jiN5ZGbLDVQpSBfq2YU6AqpnnbpQR0D1bIvMeToTNgQusg/keZ+/iYiITFqbg04ZykggItIBbQ460RDorAEFC/7nWWBZuV1wLdGFenahjoDqWacu1BFQPVuBZjbpOiQiuQ1gE8BKWlJPnwx0HcCCJoeKiLRfm1s6z/3PrHs6PQBDBRwRkW5obdAxsxO4ezUrGZv1MeVNURGRadLa7jXgYmmDfQDXli/Iek2kq7qyYKEfVfrKzJYmXZc0vou+j8vekhMAT/0FbauQ3MRl6q+ovtvTuF5Yq4MOkL2IGwAt4iZTowvnuq/POtySInNmFpqXsXa+jvuIBZjYc6twC1BuTLCKF5Lq6p/vAzgA8DJ+PkyD1gcd4MrVX9SiaeXVn0hVUct9dMFC/1ofwAtMsFXv67ANN6r0OVzQWWxp0HkBt/bXte+H2OCjDT95fKJ8XR8mfa6xwVStqGtdOhF02oTkcdu6FLrSNO9Yd0ej3UckT+GuahNzFOa93jSSBwD6bQs6fh2uvpntpLw+B+Dc//fmJLvmY3U5NLNr965jrZ3E17uqtQMJ2sh/uS9Ouh4RknP+SunQzNb8YwXuivTAvzZxsXo+N7MlM1sAEH2ZH/urz1bwdd2EW0Oqkc+65IKFku0B3KKTiXyQiS5ySi0sOQbRHMS0lDfR/MNeA3VpjIJOQSR7aN/S2PtwTfMrLQXfwtkBsFphSYlx2Ifr7riop5kNfV/1HoD1SdeTZJ/kMVxdB2h2faZSCxY2UJ8uWwRw6rsr0xSZeD52vvtvAcCnKZtEn3XregJCKOgUtwGgNd1Vvmm+CtfXn+TA/5zoTUh/Ff864/5b1F20WzDP3liY2aFvha2Z2UuEZ7koo+yChZIu6i67l7FNdJ5N/J6wmQ0yuviiLrXW9ATUQUGnAH/V9BzNfhHl6UrTvEvdHZOiBQvr8xAui0nWva/ob+aogfpUEhsl+LJt92ZDKejk8B/+vbbd7O5Q07wz3R0TpAULa+K7bVO/pP3N+Tm4+6CtnN8XDWIBsDNtw6WBKV9PpyZP2jJiaFTOkPG2NM3j3R1pLZ7WdHdMSPCChVJY9Lfcink6Ed8N/QDuIm0Z7l5tag9Bl6mlk8FfFWWNKGqlljXNp6K7Q7rPt7b7cPNeWnWBY2YnZrZlZit+rtY8yVP/HTRVFHSyrXTtaqNtTfNp6O5ogBYsHLPYzP+tLky09HWMpj60YQRqbRR0Uvi5GpPumiqE5CLJbT9h7y3chMtWdgkmaGV3R4vNamAO9Qou/U3ipNE28oFnCDeyc9IDgmqjoJPAf8CdWTKhq03zNnd3NKzJBQtnjr8YO+zQhVhc1OWcNRCnUzo/kCBkbkdGd85W3QkBx1TPtO33SAKuaV4qb1NT9Qzp7mjyWDYkCjo9pA+miKcOkoJ8touTNgYcn6WjD+DzjNGx0fmQNe+oW8yskw+4yG81PHZH9rsOl8gwqcxdd8gmX8+CZZ/79/baVk+4zMnbbT2WcJNrS33WAefyoq/Xak59zpv+O2vD8Qmo42baOQYXxBP/zhusX3Q+rucdZwCbkz6edT0629Ixs5ckbyJwbofFunX8FfSC1XijcRz1LOEI7kpqFS4tTtb+G6tn1e6OCR/LsTGzE5LRgoVpA1f6yPkM5ZLvul2w9B6LtEnVTToE8CLn+yaq59S0cDsbdICLrpI6u0v6AKIcXEl6gMs0HatDbhbiMdRzLE3zcdRzVGh3RxN1nJCHAPZJblnygoVDAE8nUrOO8XNe7mUEHMAF+EkPFNpGxsrI/nOfQzumPtRGSxuU0KZ07iSjDy71nk1UX7j7JhO/SvYjAm8lBRzfyuxlBNBGTeKz7sIibhG/1EIPE14eYJQfBHSM7DlfPbhzrQ1/x+twgefK+j+xNZSObIqWNQA63tKZcZ1qmneku2OizGyD5KoPPvEFC1cm3SXov8x34bo2e7gc2HBO8gRuVN2uTX5e2y5c3fJGbraii9XcoJ8fAGyTjLqNo2M7lVkJFHS6qzNN8w51d0yc/5Jp3ReND3qtv+LuYqvAtxRnZp6autdK8PdyFuGu2Cd+pdSFpnnXujsibe0+Euk6BZ0cvsXwBNfXMRmgBcsH+z7/bVyO6Iqa5m3o6ojfV8ozMLei6ERkdB8BrnuyLd1HIp2moCMiIo1RGhwREWmMgo6IiDRGQUdERBqjoCMiIo1R0BERkcYo6IiISGOUkUDGwk9Qjea7DDS/JZyfMzYPd0xPJp1pIoufP9aHq+8C3ByniU+olsnrZNDxM/G34LM+e0Mk51Manej3Mp5QUcZmBW5CbR/AHlqY2qVOfnLpFmI503LS/lRxD+54LsItc9DaoAMXbO7BLavRg1sXRkFHutm9ZmZ7fvb6Uuy5m2a2lPC46dOrrMF9IdS+1njbl4WeBPPLZ0+6Hk2IZYJ+EctQsV72vMjb3h/T3KU02sDMBv5YjC1jR1v+7kkYri4AAAiMSURBVNpSj67oZNCJia6ccnNj+e6dNdQcdPwV7kx8uUqqfbguxKjlsYuSSxFM8Xk0lrx1bTlebalHl3Q96JTivwTqbuLX3nKSzllELKGpv8ovu/qszqNy2nK82lKPzpipoOMd+quTuuj+kPQQfkWv86icthyvttSjM6Y66JCc8/3tca9xdWBByP57ANbr2JfMLp1H5bTleLWlHl0z1UEHbjXK0RUpDxF4VeqD2SrczWORSnQeldOW49WWenRVJ4dMl9DDyD0cvyDXRdDx67304IZ47pnZlj+hoi64Bf+++PDXJ3BDVwdw/fmjI5UO4+vsjJSx5Zeo3fQvr8Cv0+KXKf4Sly2xYbR9/HeILTAWbbOUNwciZd/XlsNN27dvMT4BcOqPSQ/A8yrzb/xieNFQ9r34sfVXjwe4nI+ykXZ/hOS2/+d7ALf89lshi6754fhLuPw9Aff5nIxs14dbxygSPwfOCo7cK3UepdQ371xNes84jtsq3Ll8CnfezKHAMPmS51Wl41XnuRtSD1+X3GNf8TspXka0bXzJ8+gzORw9l/17Cp33CfVL/T5Lqht8xTv78AfRAJwnvLboD2A/Zx9RE9ngvkQ2ASyObHMA4Dij/O0CZWz6bdfh13H3rxmA05Htj/3zvYx9ngNYr3DMDvy+54ru2/+euwnbRcOEs8qzlPeupr3mX784XinH89pn6/d5Pvr5lTiXDlLKO8j6jLN+j5LnceZ5NFJeqXN1XMfNv/9Fyme86etpSX+HVc+rMscr5Nyt83Mrc+wR9p20C2A14fm+39/ofkqf9yj5fXbt/SEHfdKP2Idu/iSKHuex5zODTmxf5g9y0gcWfUEmfWBlTn7zf6D92PPrCSdidIJcq0tsm0p/MLF9pwas0X3HT6qEky9vX1mBJeu1xbR9+z/exGOe9sdY4LgcZ32OvszNsr9HyfO4TNA5LnOujvG47Wa9z7+eFnQqnVcl/+4qn7t1fm5Vjn3Z7yRfn2sX4LHXXyR8h1U678t8n117b9UD3oYHamjpjHzAiR9Y7AtwNaX8Ml8WmVcBIx922hVrv+jvVXHfo7/necaxOUf2VWmtQQeXV1iJLbWsL92MOkZXlVkty/W0crN+j4LlVzmPyp6r4zhuRS6Oov0mBZ1K51WZ4xVy7tb1uVU99mU/59hzaeWs42qQquO8L/R9Fn9M7UACc32RuyXfdpTz+nzF6sRd6yNNsQtgMWV495qF5d3K2/dof+wRgLOEkYAAcIaaRgMWtAFgaOn3H6J7W6MDSPL2ObDs+2LR8W7LaKWy5+o4jlt036LqudjEedWGczfk2Bf+nO3y/ssxydWE3/kHXL3HXcd5X/T77MK0DyRIPCAkFy3hBhmayQ1VtIw9uP7cLbiTA8DFTdHTwDqU2rfFboqTjPKpRebRbE6tHoBB7MZlkkO4L5SiFpH/xxPt716J/Y5T2WM+juO2DFwMzimtifOqJeduyLEvW78lAK/gur1A8sTv+3nCd14d533p4zftQSftKmEDsS/bhr0vspGZDUm+hLvCiNd1HS5oVFZl3/4PZgNuRNJFxmCSG6inBZgr1jIbmtlOxqZZr6Xts+iX7WLRfbfFOI6bF42KqqyJ82qS5+4Yj30iH1hu+tFoa3AXBpsANknumB9VV+N5X+j7LG6qg07VK7BQ/opqUEP5uwBWSa7b5bDhWzX9Xkn7Xkjatx8iuYwCQ7PHrMxVeCHmhoMDxb98Gvv9azyPaj9u3gABaWDqPq+Sjtckzt2Reozr2CeVOwe47z3/N73nn48yoG+SfG9mO5M876f2nk6OcV+ZL9dRhl3mitsALuaGHITuN2PfL0a381dMfbjx+LknXo0Zd68dP7ucY1V3vqsiX57R66X7sAPUdR6N67idAJdfdmWM6by6crwmeO5e1GOMxz6t3Gv3XszlAtyAC0Lxno2JnPczF3T8iVhL1M65Aq3rCie66b8IYCVwAEHavnsZ+45yS6WVO3rSlsm4m3X80v4YtgDM+eORiGTPT5ArajdvnwAexLatVUPn0biOG3D1PsmotKBZ+bwqcbzGee6Wqcc4jn2aBxmv7eLq5zGR837mgg7cCRB6Iz4u6WohsZuqoqjraxsV+k9r2Hf0e1z78vNXhtFM5yrOkP+ldIXvNjjB1WwAo7YSRuCl8n3tJ8i+z7cOYGeMXTRjPY/GdNwO4b7Qs45bWkLM0POqyPEa57lbuB7jOPYZFjOCyDxiAXhS533Xg070YeeeOCT7sRQvaQcw7Qswej6pnG24eyNzvpzEusSeX0h6PY0/eV/icgXO2sT2vZqx7+gK58ofjG8dLQJ4Cv85VBhZd+XYxfYdzWQHYgv1xXwOYN6n9rnCv7fKVdkagOVYmpL4Pg8A/GDZKWlCu8EKnUcFyss6V8d13JBy3LZx2S0z+kUYel4VOV7jPHfL1AOofuyrfM4bCX9Tc3AX3KPncKXzvur3GQDQT/LpFJZbrnoO169EbkZXIryaBwx+H7tmtuNPzhe4+gdzAndVcnHF4E+aB/BXEaMfUkoZAyu4sqa/cnliY1hmu8i+o238f1/7n9HNSpCMjtGh7zvOPK4j+16Fu9I6wWVrK7rii38JXMvB5o/7Pbgr2ejKcjfkqszv81bsqTm4CYTXumhSfsehVVzds+J5VOpcjZVT93GL5+4awn0Z7fr6xO8VXnTjVjmvEn6P1ONVRxlFFKnHyLaZx77q5+z3uWyXudBu4fJv6tZoOQn1Cjnvi3+fdTHozBr/xTys+X6OiEjjut69NiseKOCIyDRQ0Gk53z3wOndDEZEOUNBpEZKbJI1Xc6JF4+tFRDpPQad9BvBDPP2wzuNJZVYQEambBhK0jB+6GK0q+LqmsfsiIq2goCMiIo1R95qIiDRGQUdERBqjoCMiIo1R0BERkcYo6IiISGP+P81LN7T7MOYCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fontsize = 30\n",
    "# ytick_pos = np.array([0, 2000, 4000, 6000, 8000])\n",
    "\n",
    "plt.hist(X.flatten())\n",
    "plt.ylabel(\"Number of entries\", fontsize=fontsize)\n",
    "plt.xlabel(\"Entry value of the data tensor\", fontsize=fontsize)\n",
    "plt.xticks(fontsize=fontsize)\n",
    "plt.yticks(fontsize=fontsize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([       0., 10000000., 20000000., 30000000., 40000000., 50000000.,\n",
       "        60000000.]), <a list of 7 Text yticklabel objects>)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAAFMCAYAAACArhjnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dT3cb15nn8d+T09u2ITraJjbUM9uZQNS8gAh0FllNQkovIGPSmW0S0j6zbxny7BNS/QYsyt2rLGJSPfuYZDLLSQ8p91ptiu4z25xnFvcWUYQAVAG49YfA93NOHVJEoeriCqgH99a9zzV3FwAAbfG9pgsAAEAegQkA0CoEJgBAqxCYAACt8jdNFwCBmW1KeizpUtK5uz9tuEgA0AhjVF49zKzn7me5f29KupLUc/enZtZ19wsz60s6cferxgoLAA2iK68GMdgc5v7dkyR3P5Z0FYPWRfY4QQnAKiMw1SAGoIvcnx4rtJYU/96XrltRl/WWDgDaZWUCk5l1Kjx2z8y6uX/3C87X0c0A9F782dXNAAYAK2clAlPsSntjZkdmtm1m3SxwmFknBpZNM9s3s6NZjx/vHfXjcfvxbzN3x7n7U7rxAKy6VRmVl7Ve+nGTJJnZ6H7HkrbmOYG7H5jZvqRTdz8o2P1K0lqubN/Oc04AWEYr0WIq4ULSjrtvzNtiyQ9wyHfrTfCFQred4s/jec4JAMtoVVpMkvRU0hNJ6wrBoKMQkM7yI+LmMTLKTrG78HkW5OKghnUz23T3F+5+ZmbrMZhd5YeRA8CqW4l5TDEwPHD3vabLAgCYjq48AECrEJgAAK2ySveYKvP973/f33///aaLAQC3yunp6b+5+93RvxOYEnj//fd1cnLSdDEA4FYxs38d93e68gAArbJyLaaY8WFb0oPcn7+WdEDWBQBo3koFJjPblnRf0n5+vaOYseGVmX3k7i8aKyAAYKW68rYlyd13Rie0uvuOpANJh3HOEwCgIasSmC4k7U3LYRcn314pBKfCTOQxu8OJmZ28fv06YVEBYLWtRGBy97MSiVUl6Xn8OShxzAN3X3f39bt33xrtCACY00oEphlkS15sN1oKAFhhBKabrpO5ZusqAQDqtVKj8krIDxcvWroCaL33P/l9Y+f+5rOfNnZu3G60mG7KL3d+r7FSAMAKW/rAZGbnZuYlu+bWcr+zqiwANGCpA1MMRlmXXJkl0/PDxBdaPBAAMJ+lDkwads1dKS57XiB/X4nlzgGgAUs9+CEuYX7m7vdLPmUj/jwmbx4ANGPZW0ySdGJmvaKdcsldJYkl2AGgIasQmPYkPTOzouHfz+LPp6O59AAA9Vn6wBS75PYkHU1K0Bqzi28qLH1BawkAGrTU95gy7n5sZlsKLaeBpDOFUXddSX2FQRJbLHkBAM1rPDCZ2TuSHilMaP1W0pm7/3Pq88TuufuxS6+vMDT8XNITuu4AoD2SBiYz+0DSjyY8fOHufx7Z/zeSPsv+GX+6mV1J2nT3/5WyfJLk7hcKay8BAFqoihbTP0h6N/5uCl1mp5K+kHQdmMzstwqj4ExhntEThS62rkIL6tjMNt39nyooIwCgpZIGJnd/ZWbrkv6vQqtk4O6vRvczs48k7UhyhaDVd/fv4sMvFe4FbUr6BzN76e7/nrKcAID2qmJU3kcKq8V+PC4oRfmF+B7mgtK1OBDhpULrCQCwIpIGJjN7V9KGu38+ZZ+HCgMPXGHO0LTW0J7K5bgDACyJ1C2mvqT9gn02cr9/MW3H2OJam7YPAGC5pA5MawpDsKe5Xn5idJQeAACpA1NHxctF9BS68Zg7BAB4S+rAdKWbaxrdYGb5OU5ll5Ww4l0AAMsidWC6kLQ+5fHHud+/LjqYmf1Y0smihQIA3B6pA9OJpo+iy5KoXrn7P5Y43kDDzBAAgBWQNDDF+UhmZv919DEz+7VCVgdXQUogM3vHzP4g6bm7f5OyjACAdqtigu0nkr40s1/FAPNODEoDhaCUpR8ay8yeSHojaW3afCgAwHJKnisvLmf+SNJzSU9zD2U58fr5SbVm9p8VhpBvaDiU3CT1zOyPkh7RagKA1VHJshfu/sLM7ijkw8sGQ3ytsBDfaPqhpxpOov3TyGPfU8j+8MsqygkAaJ/K1mOKAehpif0+rKoMszKzI4XEs2WHsgMAEqt1afW4KGArxWzm2QKCAICGVBqYsoEPZva1mf1V0hsz+8uEfX9jZv9iZr+oskwTzt2R9Kzu8wIA3lZZYDKznyuMrhsopCGy3PaWOALvJ5J+YmZ/MLMfVlW2MQaipQQArVBJYIrDw59rGIi+U8iN90JhZN5Y7n7h7o8kfa6wgu3fVlG+kbL2YpmKcvwBAGqQPDDFfHhPFQLSgaR77r7m7usx6LwpOkYcfPBLhUBWtYG779VwHgBACVW0mJ4ptD7uTVjF1sscJAan72K+vEqY2a5urqYLAGhYFSvY9hQm0U5aVn0WTxTmQiUXBzw8YGg4ALRLFSvYvkiVqcHd/6QQ6KrwTGHyLgCgRapYwfarCo6ZlJn1JX3t7gx4AICWSZ35oaPipdVnVcVCgXvuvlHBcQEAC6piocAHqQ4W50IlDXRmNhBdeADQWqkD05mk7YTH+0zll2AvZGZdSR13P0t1TABAWqkXCnwl6ZWZ/f2ixzKz3yksLLi/cMGG9pWotWRm22Z2YmYnr1+/TnFIAICqWyjwEzP71YTHC+8ZmdlzSR8pLJPxTYpCmdm2pEN3n5h5YhbufhAnDa/fvXs3xSEBAKogMMV5QV9KempmfxyzzPrYCbYx4esTM/tW0s8V7lelat10JG25+9Ql3QEAzatqocCtuLbRQ0kvzEwK958uJK3H5dM7CkPBOwqLCWZJVLOVbjfyK90uaKCKJuoCANKqcqHAjTgC7jfxTz0NJ8vujuye7947k/RwzEq3c8mStDJnCQBuh0rXY4rJUe8rdO3ZlE0KramteN8mSVCKSNIKALdIZS2mTEwrtBXz6PUV5jll3XgXCvOUjhPl1rshtpb6ZlYqcWx0GLseM1vuXkeWcwCAaghMmdgK+jJudZ3zzMzuldw9m8i7o9zcKboAAaBetQWmppQNLLlW0iXBCACaU+k9pkWZ2Qdm9kXT5QAA1KfVgUkh88Ommb3TdEEAAPW4DYFJCeczAQBarvAeU2ytZMuPD6alCDKz3yYqlxRG7W0mPN5YMStEPvHsjpmdcZ8JAJpRZvDDSw0nxj6U9B+n7PtY0ruLFirHNCGF0cIHNtvX+EzofUnn2WAId69iPSgAwARluvLux58mqWjo9YWmT6SddauMu++4uxVtVZYBAPC2Mi2mzzVMK1Q00fRSoYWz5+7/c5GCSZKZHUr62aLHAZr2/ie/b7oIwK1R2GKK6Xw2FJKqPi7Y/SI+Z+GgFKVciwkAcAuUmmDr7i9LHu9CUso8d99KMjN7h5F5ALAaUq9g+1TSBwkPmSV2JSgBwIpInpIoZWbwXH49AMCKaPsEWwDAiqkkMJnZfzOzX1dxbADAckvelWdmXytOyDWzrrv/99TnAAAsr6SBycweajghV5LeS3l8AMDyS91i6kna03CV2ieJjw8AWHKpA9OVJHf3R4mPCwBYEakHPxxL2kp1MBYKBIDVk3qC7StJr8zsx4kO2RELBQLASqligu3HZvaVmXXc/R8XPNxaPGaSzA+5tZc2FBYhXJN0IulM0j5rMAFA8yqZx+TuH0r6L2b2tZn9YoEWT694l3LMbFdhbakrSTvufs/d7ygsgpitwUTSWABoWOrh4u9rGEz+qJCE9VNJB2Z2ppD77iL+vcg9STtKsFCgmQ0kddz9/uhj7n4s6b6ZnUvaNjO5+86i5wQAzCd1V96GpN+N+bspBKxkLaCyzKwrqT8uKI3YkXSkEJyO3L1o7SkAQAVSd+VdavJKtE2tYLsjqWdmp/Ee01ix5ZR/DgCgAakDUzZ44I1C6+nOAts9hXtCi8qCUU9S2flVawnOCwCYQxUTbCXpYIbFBSf5Lg5GWHTo+ZHCSLwrhXlWY420pi4XPCcAYE5JA5O7vzIzSTpPdMhzLbiCrbu/MLM77n5VsGu+lXQ0z7kAAIurYrj4nqa0TGZ0IWlv0XlMJYKSJG3Gn1eSDhY5HwBgflVMsP084bG+k5TseAWyAQ97JQMZAKACrGCr63lOXUlP3Z3WEgA0KHmL6TaJc5z2FEbrbTF3CQCa13hgGkn4euHu31R8vkOF1lE22fdA0gd03wFAO1SxtPqvJb2YIcB8nPu9E1sxbyQ9SZAE9i3ufmNZDjPblPTSzC4VcuiRyBUAGlTFPabsfk0p7v4ot33o7n+n0LX2YR1rMbn7i5iuaE0hketm0XMkycy2zezEzE5ev35dbSEBYIVUEZgWTiXk7q/c/WNJJ2ZW1/LsWUvqsExwcvcDd1939/W7d+9WXDQAWB1tH5V3oOH8okrFLrxs/tWzaXn1AADVaXVgivOY6sxbdxh/ZgsKAgBq1urAZGY/V7156/IDHzZqPC8AIGptYDKznyl05Z0tcIyemZ3HrcxaUPkgWHoABwAgndLDxWcchLBjZvO0ODoKAWFdw+UqFhn88KmGAWZHxess5bsNmdcEAA2YZR7TpaT3FC70+Qmq45Y+X3TAQjayb9fd/7zAcfItoDIZz/OtJOYzAUADSgemcclZ47DqHUkPRx9aoExXkp5L2nf3Py1wHEk6VegK/Mjdy3QJ3sv9XvkcKgDA2xbK/BBzy72I92+eK7Q4XCFYzbNQ4GUciZfKc4UJv2W75bKReBfkzQOAZiRJSeTuZ2b2oaQTSe8qXNhfpTj2Itz9ysz2FCbNPpyWD8/MdjW8r7U1aT8AQLWSjcqLE1Q/S3W8VOIyFvuSXk3K6BCD0iD+c6Nktx8AoAKph4sfKkFKotRicPpA0mMze2Nmh2Y2MLN9MztXCErHku65e6rVdwEAc0iaXdzdX8Wus9aNaIvdeFsx1VBf4X7Yt5KOJB2z7AUAtEOrl1avQgxADGwAgJZqbeYHSTKzH5nZX83snabLAgCoR6sDk+pN4AoAaIG2B6aeJLn7vzddEABAPZLfY8ozsx8rBJd7mr3101HI8D0u5REAYElVEphiZvBnGk5YBQCglOSBycweqqXzmQAA7VfFPaYsKB0odMXdcffvzbpJuiPpywrKBwBosaQtJjP7kUL3XW/B5Srk7t+Z2b6knyUpHADgVkjdYlqXdLBoUMr5VpIxjwkAVkfqwNSR9FXC411I2mG4OACsjtSDH5Lmm4trMz1LeUwAQLulbjFdSHqQ+JgAgBWSNDC5+0tJ91MeEwCwWqoYLj4ws9+mOBBJXAFg9SQPTLHV9CpRcEqaxNXMunGBwNO4YOC5mR2Z2XZcpwkA0LBKkri6+1NJx2b2L2b2qwUOlSyJa1w+/UjSuaSPFLoc9xQGbOxLemNm24ueBwCwmCom2D7K/elLSf/DzJ4qDIy4UPmRe8mSuMag9MDd7408dCHphZltKmSs2Deze+6+t+g5AQDzST1cfF2hFTIaTExhKfNu4vMVMrOupMfuPnFQhru/iEvCDyTtmtnX7s4qtwDQgNRdeZfxp41s4/5WZkthoNB1N1Xsfsxac4NE5wYAzKiKeUxSuI9zf57krSOJXFMkce25+1nJfQ/iz66Z9RKcGwAwo9SBKWtx7Lv7nxIcb3+RJ8fg0o0j8DZLPOU89/v6IucGAMwn9QTbV/HXi6k7lrdoEtfsnlZH0qcl9j/J/T46UAIAUIMqhotvSSrbdVbkQtLWAsPF8+U4mbjXUH5wxvnEvQAAlUm+gq27J1vcLyZxnft47n5hZnckdUveZ8oHpjKBDACQWCUTbCdpIrWQu1/NMPjhOgHtDM8BACRUaWAys3fM7Ndm9rWZ/VUhu8JfJuz7m5gp4hdVlmmSmJIoGyDxtIkyAAAqDExm9nNJbxTmBPVUMD/J3T+X9BNJPzGzP5jZD6sq2wT5dERPaj43ACCqJDCZ2a8lPdcwEH2nMBDhhaakJHL3C3d/JOlzhVx7f1tF+UbF1lI2qXbL3QvTJsXErydmdvL69etqCwgAKyR5YIr58p4qBKQDSffcfc3d12PQeVN0DHc/lvRLhUBWh2yV3KdlUxG5+0F8Tet3796tsGgAsFqqaDE9Uxjmfc/dP87NbcqUSsoag9N3Zvbj1AXMixnFNyUdkLwVAJqXNDCZ2bsK95P6YwLSPJ5I2klwnLFiZoh9haBU2XkAAOWlbjH1Jb1w929SHCymNaokZ128r/RSobwEJQBoidSBaU3SVxUcswqnko7dfaui4wMA5pA6MHU0XPoilVTLXwwPaHYk6YKgBADtU8WyFw8K9yopzoVKmrPOzPYlyd03puzTiQsMAgBqljownenmRNVFfSbpONXB4hLr3WlBKerHDQBQsyqWvXhlZn+/6LHM7HcKSVUXWpMpd7xNhSXWi4KSFFp9qZbuAADMoIp5TJ9I+sTMfjXh8cJ7Rmb2XGE59IMUI/zisPAdd79f8ik9kV0cABpRxbIXx2b2paSnZvZY0hN3/6f8LuOeFzOPf6rQFdhRaLEsPOE13is6lDQoWMV2LZ73nsI8rMK0RACA9JIHJkly96048u2hpBdmJoX7TxeS1s3siUIQyILBevwphRbVlaSNBRYIDAcKc5WONHuXIEteAEBDKglMUhj1ZmYDSb+Jf+ppOFl2d2T3fPfemaSHcZHART3SzcX/yuL+EgA0pNL1mGLuufsKq9DalE0aLqO+nigoZYlWbY6N+U0A0JDKWkyZmFZoK+bR6yuMeMu68S4U5ikdJ8qtBwC45SoPTJnYCvoybgAAjFVpVx4AALMiMAEAWoXABABoFQITAKBVCEwAgFYhMAEAWoXABABoFQITAKBVCEwAgFYhMAEAWoXABABoFQITAKBVVjIwmdmumY2uCQUAaIGVC0wxIA0kvdd0WQAAb6tt2YsmmVlXYS2oPc23oi0AoCZL3WKKXXYu6VTShqQdSS+aLRUAYJqlDkySsqXV77j7lrsfN10gAMB0lQYmM/vMzP5SsM/vzOyPVZzf3a+qOC4AoDpVt5h6ku4V7NOVdL/icgAAboll78oDANwyBCYAQKsQmAAArTIxMJnZj8zsD2b292b24zoLBQBYXdNaTGsKc38+kXRkZn81sz+a2RMCFQCgKhMzP7j7S0nfM7OHkj5UyJywHrddM5PCxNWXko7c/Z+rLy4AYNkVpiSKAepl9u+CQHUm6UjSF+7+v6socFuY2bakbUn6wQ9+0HBpAGB5zJwrryBQ3Y/bXgxU2T4/dPd/Xbi0LeLuB5IOJGl9fd0bLg4ALI2Fk7gWBKofSXJJFzFQHUk6lnRG1x8AYJzk2cXzgcrMvpL0UNIzhYEUH8bNc11/RZkhAAArpJZlL9z9Y0kys3cVWlIPFALVdSoiM/urhi2qY3f/cx1lAwC0S63rMbn7d5K+jNsnkmRmJwpdfq/0douKQAUAK6YNCwVeSpK7/12uRfWhQhcggQoAVkwbAtO1kRaVCFQAsHomBqY4uu4rNTiJdsZAdebuD+osH4DJ3v/k942c95vPftrIeZHOtBbTpaQ/a3y2hy8kvay7lTIhUG0oBKvzOssCAKjGtJREf1IcNTch20Pj3WkxUL2IGwBgCZS6xzQyN6nMfZ8j5Sbdtkxn5CcAoEXmSUlU+r5P9hwz+4VC1983ixd5PmbWUShnP/7pkZkdSjpx96umygUAuClFSqJpgWpL0ruKOeXM7ErSiWKrqupEr2a2KelwwsOdWA7l8vrdd/ezKssEAJiuipRE14HKzD5QaEV9rBCofqQwWGFDw64/SZKZ/afUgcrdX0iywh0BAK1RV0qiZwr58kZH0q1L6il0+501PZgCANC82ifYjo6kK5HolUAFACukNZkfxiR6JdsDAKyg1gSmzIzZHrLBFF9J+rLJUX8AgDRaF5hGTQhUjxW6/R5qOJjiY0n/oaFiAgASaX1gGhUD1fWy5nHkX08Sc5EAYAncusA0yt1fKazlBEzVVFJRALP5XtMFAAAgj8AEAGgVAhMAoFWqvse0L6ko91yZfQAAK6LSwOTu18O8F9kHALA66MoDALTKrR8uPisz60nakbQW/3Qp6dDdj5srFQAgs1KBycz2FTKaf5Stu2RmXUn7ZrYnaYtFAwGgWSvTlWdmRwo59x7mFwN09wt334j/PG2kcACAaysRmMxsWyEoTWsRbUnqxlYVAKAhKxGYJA0kXUxbNj0GrGNJ27F7DwDQgKUPTGa2KamjuDBhgcP4c6e6EgEApln6wKRhkDkvse9F/LlZUVkAAAVWITCtx58XU/e6uU/XzDoVlQcAMMVSB6Z4rygLMIWByd3z+6xP3BEAUJmlDkyS8oMYLhd4LgCgJssemBbpjqMrDwAasFKZH0q6UghK7zVdEACza2ql4m8++2kj511Gyx6Ysnx4miPVEC2mirDEOYBp6MqbbK14FwBAasveYqpMTHO0Hf/5/8zs/zRZnhb6vqR/a7oQLUJ93LR09WGDuZ+6dHUxgx+O++PKBCYz68zYnTd1FJ+7H0g6WKxUy8vMTtydIfcR9XET9TFEXbxt2bvyFlnCguUvAKAByx6YZp27JA3vS32bsiAAgHKWPTDlWz2zTpilxbQYujlvoj5uoj6GqIsR5u5Nl6EyMd/dm/jPjaLl02MKoyzZ6/1py2QAAKqx1C2mONgha/mUaTHlh5eXSfoKAEhsFUblPVcY1l1mTlMWvM7mmJB765lZT2GZkGwO16Wkw6KW5gLn68bz9RXq/lLhC8GhpOdN/x/UXR8FZTmSNGji3LkyNFYfcV21HYX3yZrC++QLd39a9bmnlKn2+ojTVDYk9eJ5TySdSXrS9OclKXdf6k3houeSjkrsexj33W263A3U076kU0m93N+6ko7i1kl8vl2FbtNthQ9ZV2EdrOz/wCVtr0p9FJRlM9bH5qrVR3xvnMdzbGbniX8fSDpdhfqIr/c0vubeyN8Hy3bdarwAtbzI8B/qRW+WuM+bpsvbQP0cxQ//2PrJHk94vl2Fb5aTHs8uxK7QSljq+igoS3aftLHA1FR95N4HEy+48b1U63ukgc9LL74HulP26cZ99pt4jySv46YLUMuLDP+xUy9y8Zt7o99KG6qb7HX3puzTifss/KaPH6DCb7nxguN1/5/UXR8lyrPfRD00XR8a9nRM+8xmLYXaWk1N1EcMghPPN6bOCvdt+9Z4AWp7ocMLXX/MY926LjRt2+K3rMJvd/FboE/71lbyfIdlPzi5lkItrZMm6qPgHFk3zXmDgamR+ihz3lV4fyi0Gkv34sT3yq2/ji31qLw8DzdJNyTtm9nAzPpm1jOzXYWuvh1332m2lPWKN5Q7kl6U2P0w/ly0jnpefhh+Nr+jG280V6qh+phm4O57FR5/qqbqw8z243mLss/tKdz4r6WOGqqPDYUBDmUdawlW316ZwCRJ7n7s7vcUvs30FJq+F+5+x0Puu1WTfWjOp+4VZMPnN+c9WQwuXTN7Ez/kRfLlquPDVmt9TBO/MM2fFjSN2usjzj3clq7zUU7k7gfuft/dywSKFJp6f8z63p8n402rrFRgysQA9TRudb2p2yh7w5eZs5Xt040Xj3lkw/E7kj4tsX/+m+K9Oc85i7rrY6x4vAfe4NDwqIn6eBR/tnFyexP1cS6pE7+olNFX+OJ9q61kYML1HKLsA1P4QXP3/D7ztl7yF5sy3RP5SdFlvqXOraH6mOSZauqemqTB+shaJU0H5RsarI/sMzOIc5gmMrNsPuCt7/1ZhQm2GC9/0Z+16T9r3kFJ4cNqZncUbgiX+UacP88s/ezzqL0+xokXl69HLmxNaKo+snuJ119EzGygEBTW4s8LhftvddZRI/Xh7sdmdqZQL/tmtiHpIx+ZTBtbZYcK98pv/URbWkyra5Huhbmf6+5XMwx+eJB7XtVdO43Uxxh73mA2g5za62NkgMulmXXM7FBhlNmOu2+5+4bCYKXzohZEYk2+Px5qmFptU9Kr+AVG0nVQeqnw3rn1rSWJFhNmc6XwIXuv6hPFD1t247gNF+pxktZHbBk02oW3oEXrI9+yuFJoFW2N7uTuB/H9sW9mhYMkGpTk/eHuV2b2gUKLqB+PeWRmBwr3k3YUWlFtvC83F1pMqyvL76U5mv5Jb/ZPkP82/KSG8zVaH9k9jBZdXJqoj7Xc7z1NuYkfW5UXCsEpWVfqFI2+P2JPw4ZufknbVghWgxa9b5IgMK2uRT4sa8W7zC9+G86GSm/V1GfedH3sq12tpSbqI3/OnRIjZrPH6xhW3/T7I/OtwoCI/GfiKLa2lwaBCW30LP5cieH88V7J4TLctE6ozKi8rEW1mXrIftuYWdfMTiVdxblbd3Rzou+umZ0uSz0QmKA53syVTeCLF+lNSQdNZT2osz7iubZafJ+kzvrIB+ZZJrFK4d5LLer+vMSBDqcaGdwQ779taVhvPUlLEZwITKtrkW/nlXyzj6Oy9hWCUt3poZqqj4GqTWs0r6bfH2WGgucv+A8m7pVGI/URPxNHCl9e3mpFxh6FDzRsYXY1TId0axGYVtc83+Kyb2LfpiyIdGPI64uGchbWXh/xonPVgjlL4zTx/sifs20Dcpr6vBxKOp6WBWTMwIh+fjj5bURgWl35D/6so5qqaDGdKnwA3xoeXJMm6qPRJK0FmqiP/MiyWQNNpQNy1EB95DI5lHqPxPdS1tXX1OcoCQLT6spnUij8UI8MyU2ahSEuG37RYFCSaq6P2Frqm5kXbRpeCA9HHqskgWxU+/tjpOU4a6CputXZxOdlQ5p5cnkWxOoYQl8ZJtiuqDhpL5sAWOZNnP8Gm+wiEJc4UOyKmLRPR9JalV1eddeHu5+ZWdnEtNlAgB3lRqstU33kZOl37pfYNx8gKs2l2FB9dGd9biznmW55YKLFtNqex59luk2yN/pZqmHNMWNyd1pQivqqZ9RVrfXh7hdlttxTLif8vSpNvD/yN/HLnjP/vCrVXR+Xmi/AXKj6FmSlCEyrLRu9UxQYJOlx/PlFihPHbqjHJYKSFEZc1fFBa6w+WqqJ+tiPP8tk5M7KdVZToK67Pk6lt7oFy5iaNeNWqGOZXLb2bgpvfldIhzNtP1fBEhLJ/moAAAufSURBVM8KH9w3Cjf1p+3Xk3Q0QxmPisp3m+ujRJlczS2t3sT7I1uWvF+ybLXVS531odAyc4XJ12XL14/HrOXzUlk9N10AtobfACFI+LSLhUJOrqkXgNw+2dabsF9X4X5ANpF20rYtaVfhG7Qva32ULFOTgan2+ojvEZd0OmWfzVkv2rexPnL7bZcoWzcGpdrfJ8nruekCsDW/xQAw9htq7iKxX3CMwg9a/AZ4PrJfmW3iBeo210eJcnRyZXGFlkR3Wd8fI/tPDDwKrQLXDK3u21wfufMdakJLKB7vTZkAdhu2xgvA1o4tftjPFTIR9OM3w91Z3uwadsGM/VCO+TCW3Wr9VlxXfUx53n6Zeln2+sid802sk0HuGLt1vyearA8NMzq4Qndivj7exJ9zt8rbtll80YCk60l92YJtF74CSVSnoT5uaqI+sjlfdZ6zrLrrI06d6CoMDslW8z32JUsATGACALQKw8UBAK1CYAIAtAqBCQDQKgQmAECrEJgAAK1CYAIAtAqBCQDQKgQmAECrEJgAAK1CYAIAtAqBCQDQKgQmALeSmQ3MzM3sNCY3xZIgiSuAWydm9c4vH37g7jtNlQdp0WICsFLMbNfMBk2XA5MRmADcOu5+LClb++hCYdG8st5LXyKk9DdNFwAA5uHuW3M+tasQzNBStJgArJpu0wXAdAQmAKumV7wLmkRgArAyzGyz6TKgGIEJwCphSPktQGACsBLMrCep33Q5UIzABGDpmVlX0mHT5UA5DBdHEvGD31EY8bQm6cTdz0b26cfHOwrDdc/cvfSw3THnOM4/Pz7ej/ucxbkus7yGfPmuxr2GOcqYvB5KnL8Xj38l6VKhnq4qKPtcdRSP05G0PnKsG2WN5dFo3ZjZrqR7uXJ2JT0fl/khHmNT0qfxPLgFaDFhIXEWvUs6l3Sq8K10X7kuk5jT7FzSloYXhx1J52Z2XnRDeso5erl99uPfMkdmdlqi/B0z2zezNxref7iI5RyY2Zt4ISw6TuX1UHD+fny9p5IeaDhP57GkNzGf3NjRaGa2OaHsj0eOfyRpQyGAXEi6L+k01lGpLjIz68XjHCoElBOFibIXCkHmZaynbixLb+T5HYUgs61Qtz1NCDjx/+1cYfLt6D67Mc/euG03d4zBlP2ybWLAi/U26Xn7k5638tydjW3uTeHisilpV9IbSR63XYWLxqmk7QnP7eWec1hwjm2FC0z+HJvx8SNJ/dz+h7l9xp477rcd9zmS1JlSxvN43l6T9TDl3NnrHUwp23lWnoKyn+fKPoiP70ran1L2G/8fU8q5WXK/3TLHHNnvrfLF19XPbdu5/Q9HHstv3RL/rx7fj0WvpTNSTo/vhd38edhG6q3pArAtzxY/1NcXiniRm/rhix/87AN/OuM5NsddHGKguQ4ME46TXSwKA0G8uGQX7P6MZaykHsa81rGvc+T4ZS72+UAziHU8NijlnrMf939TUIeF5RxzzKIL/8TAVFAHY4N4wfPzdXM+43N3y5aTjcDElnAb+eC7prQwRp7XL3uByV3gstbIWxfxeAE5nHQsDb+5l764jLy2MkGm0nqI+w9mCWS5i71rQgsx7pcPTIV1NFLusYE7d2GeeN4JddiawBSPcVSmDsc8b3vWYLbKG/eYkNJl7vczL3lT3MMghWygwvakeyFj7OjmfaXseGfuvuXjb4Z3NBydVTrxp4cb8FnS0KLRXZXXQ3wsuxfypMzxdfP1bpfYf1tj6neMk9zvk8r8oMRxrsX6bmM+u1nrMLMhaS9xWZYWgQlVOSne5Yb8BbBswOhKej7jeT7NfnH3gxmfm5WxV/Zmv6qrh/zreDFlv2vxYp+Nzns8bd+o4+5PSxw3P+LvXsHus1zMZx7tV7X45SF7vaUm68YvQ72y/08gMKElRj60/ZIrkp55iWHQI7JWxswXPb85/LySb79l6iH+LRvBN+vryFohZVql87RY1ib8/ev4c1BmlGP0xZxlqFrWQu2W/IJStuWJiMCENslfZMt84Ge6aI1cROa94GWBsMoMAkX1sMjryM/7Kgr+KYPCgYZ1N4jD4wfTuivd/UXZbtCa5VvaZb6g7Iw8BwWYYIs2OdHwm/wDDe/pTDLrhTN/EbycuNd0l4pzYsysV9GFs6ge8vdrsvlFZa0p1luJ1uasrdGJ3P3KzB5KeqnhJN5dhflEUgjGx5KOfMaJ0XWLr+WFQqu1b2adSXUZvwzN07JfaQQmtEn+w1tmzZxvZzx+/v7HvBeK0TJWEZiK6iHf0jlw91txU93dz8zsAw0nyOZfRy9uu2Z2JWlvjnuAdXqiYXfqp5rcctqb8hgmoCsPbTJroJnVpPsfbTlepqge8ue9VWl23P3K3ffc/Y5CBoyneju4dyTttzkzQmwpZy32sQM6YvaKbku7I1uNwIQ2ybdo5u1qmybf9TfvBT3/vCrKKBXXQ/51VBUcKxfvIe25+313N4UUR/lRgNszDJRoQjZqsjMhndSOZpiSgCECE9okf5E9r+D4+ZbIvBf0/POqGjFWVA/5v92KFpOZHRXNT4vzz/Yk3dGwFTUoOUIzGTM7LDPabqSr8dMxu2y2vDuytQhMaJP8xaCKG+D5QQRl7mGNk10kryrsoimqh/zf1isqQ2prKjmSMQ4UeKjhvbZHVRUqgSzw9LJs6NL1SrmtHsTRZgQmtEJuyQWpoov+SDaBstklro18i67kolOmHkZeR2eeFkVDXWQbZXeMwSmbL1Q0abdJ+a66/CCHT8XcpbkRmNAWexN+T+36QjJD9oZM/sJaVRnL1sPoRbC0GPyaWGJ81vrOgnKqQTH5+3XTgnm2Tlah+CUhK+e2dGMdKQY9zInAhKqU7n6JH+RsZNNZlf3y8djZBWPW4JKV8cDLL+xXST3EDBE3LogzGKihm/JmNst5s66xJBf42ArLugendeV2Z/j/lXK5Cs1sW+F9VTZ/IcYgMKEqnRkuQlmXR3ZvoWrZ/Yt+2cX54tDlbGXcWVobVdZD9jo6ZlZq2fDYSuw0eFN+d4auxy2F7syU3aZZbsWx9+bil4OZ5rjFLwnZc/YUMqyTF28BBCZU5UDhhvDUb/PxgtpXXBG1xAz5/Ii1ue49xHPcV7iYHJZYQXeg0Co51uyBs6p6yL+OC0mbRfN+YlDaV/luvFIBZCTQFD3nTMXZ2bPBA32F4FRG2VGWWR1NGuI9b2snC/RdFWcsQZGm191gW55NN9dK2o9/21dcRntk377CSp7ZaqLT1gfqKAxW2NTN9XCyNZmyJbZnWhE0Hnc/V4beyOObuTKWWuCuynoo+TrO9fbCidn6VGNX6s3Vb183V/+9Ub+LPCe+xtP4+27891trNylc2A8VFk18ay2mgvNmy60XrZeVrQ11fQ4NV5qdeQXhMf/nrEy74GaxUoGFxW/Ob+I/Dzx2ecX5K58qfHjX4s8rhRbIvhf058cRZGW6w164e9lv2PnjZ4MB+rp57+FE4eL33GfIdVZVPcz5Oi4VWin7PqFLLLYcClsxHibBZs8plTE7e05szQ2y1xjLuqe3u9QuNaXOS5a18H2Q+7/oKfxfXMbnzT2oJctZ6O6lRx9iPAITkpl0QU517GnBoejxOlVZD1UpSETakd5O+jrPc+ooa1Pvg9gdO/ELAMojiStuhaKLTVuC0m01rf4mPTbPc1Jo6rzT5BYDJCglwOAHAFgciwEmRIsJABa3ozBCEgnQYgKABcQBGSwGmBCBCQCmiMvA+5T5bgOR6SEpAhMATBCDUTb0/q3RlfHxCycvXlLcY0JKaxN+XzXUw/LIB5wb86fiSLyBuLeUHPOYsJD44ezGLZvcmdlT+GBfLvs3SuphecXJwV1JW9l9pDhB95mkj/g/TY/AhIWUzcqQzxqwjKiH5Ra77LLlRS4V8hPuMeChGgQmLKypDABtQz0AaRCYAACtwqg8AECrEJgAAK1CYAIAtAqBCQDQKgQmAECrEJgAAK3y/wEV6BvLPaRqYwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fontsize = 50\n",
    "# ytick_pos = np.array([0, 2000, 4000, 6000, 8000])\n",
    "\n",
    "plt.hist(P.flatten())\n",
    "plt.ylabel(\"\\# entries\", fontsize=fontsize)\n",
    "plt.xlabel(\"propensity\", fontsize=fontsize)\n",
    "\n",
    "xtick_pos = [0.2, 0.4, 0.6, 0.8]\n",
    "plt.xticks(xtick_pos, fontsize=fontsize-10)\n",
    "plt.yticks(fontsize=fontsize-10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "prop_B = std_logistic_function(A_true)\n",
    "mask_B = (np.random.uniform(0, 1, prop_B.shape) < prop_B).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_sq = square_unfolding(A)\n",
    "_, A_Sigma, _ = randomized_svd(A_sq, n_components=10, n_iter=20, random_state=None)\n",
    "alpha_of_A_sq = np.max(np.abs(A_sq))\n",
    "theta_of_A_sq = np.sum(A_Sigma) / (np.sqrt(np.prod(A_sq.shape)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.1740064600780293"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha_of_A_sq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7183197378579977"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta_of_A_sq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate_sq(max_iter, fixed_step_size=False, tau_ratio=1, gamma_ratio=1):\n",
    "    \"\"\"\n",
    "    A helper function to run Algorithm 1 with a specified configuration.\n",
    "    \"\"\"\n",
    "    print(\"maximum iteration: {}\".format(max_iter))\n",
    "    \n",
    "    prop_B_sq = square_unfolding(prop_B)\n",
    "    mask_B_sq = square_unfolding(mask_B)\n",
    "\n",
    "    start = time.time()\n",
    "    A_sq_pred = one_bit_MC_fully_observed(mask_B_sq, std_logistic_function, \n",
    "                                          apg_max_iter=max_iter, tau=tau_ratio * theta_of_A_sq, \n",
    "                                          gamma=gamma_ratio * alpha_of_A_sq, max_rank=25, init='uniform', \n",
    "                                          fixed_step_size=fixed_step_size)\n",
    "    elapsed = time.time() - start\n",
    "    print(\"maximum iteration {} finished\".format(max_iter))\n",
    "    return deepcopy(max_iter), deepcopy(elapsed), deepcopy(A_sq_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "change tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau_ratio_all = [0.25, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = mp.Pool(30)\n",
    "result = [p1.apply_async(evaluate_sq, args=[10, True, tau_ratio, 1]) for tau_ratio in tau_ratio_all]\n",
    "p1.close()\n",
    "p1.join()\n",
    "\n",
    "result = [item.get() for item in result]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_true_sq = square_unfolding(A_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_optimization_details(item):\n",
    "    max_iter = item[0]\n",
    "    elapsed = item[1]\n",
    "    A_sq_pred = item[2]    \n",
    "\n",
    "    loss = tensor_log_loss(square_unfolding(mask_X), A_sq_pred)\n",
    "    _, A_pred_Sigma, _ = randomized_svd(A_sq_pred, n_components=5, n_iter=20, random_state=None)\n",
    "    alpha_of_A_pred_sq = np.max(np.abs(A_sq_pred))\n",
    "    theta_of_A_pred_sq = np.sum(A_pred_Sigma) / (np.sqrt(np.prod(A_sq_pred.shape)))     \n",
    "    \n",
    "    return max_iter, elapsed, normalized_error(A_true_sq, A_sq_pred), \\\n",
    "        normalized_error(std_logistic_function(A_true_sq), std_logistic_function(A_sq_pred)), \\\n",
    "        alpha_of_A_pred_sq, theta_of_A_pred_sq, loss\n",
    "#     print(\"iter {}, elapsed time {}, A normalized error {:.4f}, P normalized error {:.4f}, alpha: {:.4f}, theta: {:.4f}\".format(\n",
    "#         ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p2 = mp.Pool(25)\n",
    "result_extracted = [p2.apply_async(extract_optimization_details, args=[item]) for item in result]\n",
    "p2.close()\n",
    "p2.join()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_normalized_error_all = []\n",
    "for idx, item in enumerate(result_extracted):\n",
    "#     print(tau_ratio_all[idx])\n",
    "#     print(\"iter {}, elapsed time {}, A normalized error {:.4f}, \\\n",
    "#     P normalized error {:.4f}, alpha: {:.4f}, theta: {:.4f}, loss: {:.4f}\".format(*item.get()))\n",
    "    P_normalized_error_all.append(item.get()[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbkAAAE8CAYAAABdKc9DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dT2wc150n8O9XluMdKGNTVJxDAjhS00AuAcam6Nx0sEVqc1msE5ASAuQ2Jqn4lIOHtAcDzGEOMuUBdk5JmlIuixwikZr4soAzpJyDbzHZ8mKwh0xG1BhBAkwyEukgQuBY0m8P75VYanX9fdVd3cXvByDaZr/qV02R9ev36vd+j2YGERGRJjpU9wmIiIj0i4KciIg0loKciIg0loKciIg0loKciIg0loKciIg0loKciIg0loKciIg0VnCQI/kGyVeqOBkREZEqHQ45mOQPASwA2AVwrJIzOoC+8IUv2PHjx+s+DRGRkbG9vf1fZvZsVrugIOetA/iwgtc5sI4fP46tra26T0NEZGSQ/DhPu9Ag9+9m9o9FDyL5tJn9IbBvERGRVKH35EjyhRLHbQf2KyIikikoyJnZOwDOl0g8mQjpV0REJI/QxJOnASwDmCe5DKADd39uD8CdhMNeAqD9fUREpO9C78ntYT9gEcB04OuJiIhUJjTI7QBowWVYJo3cuk0A0Lo6ERHpuypGcgtmdrnIQSTvB/YrIiKSKTS7cgtuNFfUjcB+RUREMgWN5MzsfMnjpkL6FWmaj2/fxaUPdvDujd/i7qf3cOSpw3j1xS9h/lQLXzl2RH2r78b0P+i+aaZEx7pNTU2ZKp4cXD//5e/w+o87+Oz+A9x7sP/3ePgQ8eQTh/D970zi5a9+UX2r75Hvv8q+SW7nGTBVugsByVdIXiH5K5K3/eMVkt+ssh+Rpvj49l28/uMO/vTZ/Uf+6AHg3gPDnz67j9d/3MHHt++qb/U90v3X1XclQY7kcZIfAtgAMAeXQXnUP84BWCf5byT/qor+RJri0gc7+Oz+g9Q2n91/gMsf3FLf6nuk+6+r7yq22jkBtwh8AsAlAItwgW3GPy4DuAbgCwA6CnQi+9698dvHPtV2u/fA8NMbv1Hf6nuk+6+r7yp2IdgAcMGX+EpFcgXA+9C2PCIAgLuf3svX7s/52qlv9T2s/dfVd9BIjuRrANbzBDgAMLNlAJdJvhHSr0hTHHkq3+fMI5+r4vOo+j7Ifdfdf119h05XzpjZm0UO8IFuJrBfkUZ49cUv4fAhprY5fIj45otfVt/qe6T7r6vv0CC3W/K4TwL7FWmE+VMtPPlE+p/hk08cwmunTqhv9T3S/dfVd2iQK7vITovzRAB85dgRfP87k/iLJ5947FPu4UPEXzz5BL7/ncm+LJJV3wer77r7r6vvoMXgJK+a2dkSx/3MzP576Y4bRovB5ePbd3H5g1v46Y3f4O6f7+HI5w7jmy9+Ga+dOjGQChTq++D0XXf/VfWddzF4aJD7GwA3zeyfCxzzGoDni97LazIFORGRYvIGudA0llUA2yR3zeznOU5qHsDbAPoz4SwSoO56giJSvdACzZ+QPA/gOsktAFfhFobfgduGZwzAOIBJuEXiLQBnzOwPQWctUrFeNfX++Ok9/OQXv8a17d/0vZ6giPRH8IIEM9skOQXgOoCVhGaEC3pnzOx6aJ8iVYrX1Ot274Hh3gNXU++9753SiE5kxFRSu9LMOmZ2FMB3AXwEF9Sir1twpb1OVB3gSE6SbJNc819tktNV9tGjzwXf102SuyQ3SK6QHOtnv9I/ddcTFJH+qXQXAjNbNbOTZnYo9vW8mb1jZpWujSPZhquV2TazOTObgxtJLvvAU2nQ8QF1G65G5wUzm/CBfdk32SW5VGWfMhh11xMUkf4Jnq70Jbo6ZvZ+BeeTt88NuPt7J81sL/q+me0AmPHPRwGpiv4m4aZjT/o+HjKzDlzh6TZcEs6EmS1W0a8MRt31BEWkf0JrV/4QwEUAa9WcTq4+FwBMA5iLB7gucwBaPvBUYQ3A6e4AF+efmwOw4IOijIi66wmKSP9UMV25DrcsYFBWAOz4EVRPPvhtwgWcVkhnJGcBjKf1F+t3E8AOXCapjIi66wmKSP+EBrl/N7OzeXchiJB8ukxnPuCMwQXWLNHoMjTgzAAoslJ7E0DmAkUZHnXXExSR/gkNciT5Qonjtkv2FwWsmznaRlOLsyX7iisatO5U0KcMSN31BEWkf4KCnB/BnSf5SsFDyyaERMEm8d5YTNSmFZhpeRPAWIHMyWm4jWRlhLz81S/ive+dwre//hw+/9RhkMDnnzqMb3/9Obz3vVNaCC4yokJrVz4NtxZuHm5arwPgQ7iF30mjmZcA/NDMnijYVwv7I7iJtCSQ2DHRm5vx98sK8+vuoqC1aGarOdoeTUmKeYxqV4qIFDOo2pV72N82h3CjmH6JJ5AUnQ4snXziK7p04EqTtUnOAJjvDmJ+tLgGFwhzBzgREemf0CC3AxdA1pE/8EwAKDq9CbiEk7JCF4afhqvcMgZ3j2+a5Fw0OvQB7jqA5bSRnoiIDFYVI7kFM7tc5CCSjxcJ7I+oSPSxkBcxsz2SJ+BGatP+NTdIrsJNTy7Cje4ylxmIiMjghGZXbiFfEki3GyWOGY/+o8R0YHCJLzPbM7MZuMXvkQW4wLeiACciMnxCsyvPlynnledmYQ8hgWo8u0lut+ESbOKBdoNk0g4MIiJSk+CKJyTfKLGEYOSQbPkCzXu+CPVRPLoofYnkdt7lCn43gy2SW7///e/7cs4iIgfdyNWu9P0WHdUFLc72SwO20ZVY4nc+mMP+qG4Srkhz5vn5HRumzGzq2WefDTk9ERFJMEq1K0PS8ksf64stb8AVhH5srZ2ZrQM4AVfOC3DZpgMN+iIi0tso1a4sMxqLRlS3SxwbWQOwmbaYvEdSynS/N28VEZFso1S7Mj4aK7q4u9RIzgeqFvY3Rk1lZssAounMuTJ9iohIdUapdmW87lVmtmTXFjtla2bNAA83Rs0rCohBW/yIiEi4oMXgftpxGcA8yWXkr11ZuGCmX5AdLe7OE0DiyR9l1vLB91PoWH+eHSjIiYjUbpRqVwLAVbgF2HmyK6Mg0wmoJXkH5YLVDipYgC4iImFGqXYl4JJAFuCmES9mtD3nH6+U7Avw9w5JtvLsehAzCaAd0K+IiFRgpGpXxnYEmCY5ljFCm4VbuJ0YDElGtShXfdJIt6twwWoFORNJfLLKOPYTUEREpCajVLsyMu8f30pqQHKhq21Sm1m4acUlvx7uET6ILgKYjb1mIp/ssoYeW/GIiMjgjVLtyujYDlyyy1KvtWg+0LThRmfr3c+X6G/V99cmuZZUzcQHwagqSnC/IiISLnS68hE+27IFAGb2Uez7xwHcMbM/VNGPmV3005ZtkutwFUnuwE09voWMHbz9a6ySnMP+dGXiMgHf3zrctOWu73sLbrp2EsCU///T2o1ARGR4VBLk/Dq5FbgLPuCWErwUbwLgMslfmNk/VtGnr0Ay4UdzUb87vnBy3teYKdB2B8CcH8m14ALbGNyocU7TkyIiwyc4yJG8gP0F0B24e3SPrIMzs1sAzpI8TfKCmSXeTyvKB7vEkltV88Gs479ERGSIhe5CcBouwC2b2SFfVf8sgFu92pvZdQBbB2FrHhERqV9oduUCgJkeBZoTK5qY2TWorqOIiAxAaJA76kdnRTGwXxERkUyhQa5sTcjcySEiIiJlBY/kSh5XZhcCERGRQkKD3C7Jl4sc4LMxB5YNKSIiB1foEoKLcNmSf21mP81qTPJtuFJbJwL7FRERyRQU5Mxsh+SbAK6R3IIboW0BGIvtGN6CWxgebZFzpqrKJyIiImmCF4P78lh34KruT2F/+UC8oDHhSmCdKZmNKSIiUljoPTkAgJmtm9k4gDfhdhhg7OsG3GLxcQU4EREZpEoLNPu927I2MxURERmISkZyIiIiw0hBTkREGktBTkREGktBTkREGktBTkREGqvS7MpBIzkJYBHAuP/WHQBrfiPVfvc96/tu+f53AFzxGaYiIjIERjbIkWzDLT6fN7OO/14LQJvkMoA5v4t31f1OAliDC2ptAJtmtue/f47ktpmdrLpfEREpbiSnK0luAJgGcDoKcIArM2ZmM/5/t/vQ76x/3baZzfhF8Hu+746ZLQO4QnKl6r5FRKS4kQtyJBfgAlzaSG0OQMuP9qrqdxpuBHcxaUrSB7cVf34iIlKzkQtycEFkJz6C6+aD3yaABT+FWYU13+9ySpuoXudYRX2KiEiAkQpyfrpwDMB6juZr/nGxgn7bvt+sachlAB3/KCIiNQsKciRfq+pEcooC1s0cbXf842xIhyTH4EdoZraa1tbMVs3spJnlCcIiItJnoSO5Nsm/rORM8pnyjzuprR5t0/KBqqyz/jFxelRERIZTaJAjgEtVnEhmR+7eWhSsMoOcmcXbTCU2zBaNHvu+9k5ERKpVxTq5syRPwu0Z988VvF6SeALJnYBji5r0jw+nSH0W5RjcIvAxuKC70hVYRUSkZlUknkzCjZRI8irJCySfruB1u4VMOZY61i/wjtwhOUZyDW6d3KKZzfl1edsAbvrlDSIiMiRCR3IzZvaR/+9rAK6RPAHgb/1j28zeD+wjxB5cgDtW8vj4CHAPbrQ2193IzFb9fb82ycwEFRERGYygkZyZXe/xvVtm9qaZnQNw1I/u3qhgdBfVp0TKIvAkZUeB47H/ngSwkdTQLxDfgQt0Va3NExGRAH1dJ2dm18zsLNwo7yLJKyRfLvlyIdOV49lNMvtczLE0IHo+s6wXyQWSWyS3fv/735c8PRERSTOoxeDmH+cAbJL8VUWju0HKk10ZjfRms5Yt+DV1U2Y29eyzz4afnYiIPKavQY7kt0h+CJeZuAB3X+sdAGcAXIcb3f2A5CsFX7foqK5oNmYkPi1aZAE6oPqVIiK1C6148q0e3zvuMyxvw5XWOgngBlxB5XF/v+6Wmd0ws/Nm9l0AJ/1U5vGU7kK2zaliy508ywPiwfSlCvoUEZEAodmVawCeAB4GvEW4EQz986twGYm30l7EzN7xr3GV5JKZ/UePZmVGY9GI73aJY7v7HFSyi4iIVCQ0yJHkr7C/KJpw5a/aZla4EoqZnSX5AwDf7fF0PMi0UKzMVtmRXLyPokGrbLKLiIhUpIqKJxP+cRUuuN0IfL1PEr6/FfvvzADSlca/ldgwhZntkNGgtHDQUvUTEZGaVZF4smJmh/z9tVIBjuQzPgHlBQDP9Grj18ZFI7I869DiI6+QgBON5k7maBsPhHkSVUREpI+Cg5yZvVXBeSzA3c+7DlciK8lV/5hn6jAKhJ0Si8fjoqUDeQJrvI0KOouI1Cw0yP3vtCdJnvbJJD/wGZdJSwVW4TYaPWtml1NeMtoIdSbHuZ3zj1dytE3T9o95djKIzqujYs0iIvULDXL/I+1JM7tuZmf9MoG34cp8XeleBG5mn5jZO73KhHW124SbPpzOsVZuFsCeL7fVE8k1krt+V4GkPnfgRmVjJLPWvkXPX8hoJyIiA1DFfnK5+EB2zde0zCx7lWLePyZOk8Z2A5jPaDMLN/W51LXjQLdoT7nE8yY5C1ffcl07g4uIDIfQIGfZTaplZh24qc2lXiMrn1XZBrBaVbDxo7k5AJN+q53uPqfhplI3e+1SICIi9ahiCUEhJJ9B2E7dMLOLJDtwFf/X4WpG3oGbLnwLrphy6nY3fnucOX/Mqg+eae3XSc74PnfhkmD24EZv03CbxiZOjYqIyODlCnI+tT/pub9CvmnLFlypqwXsZ0mW5u/PTfhRVDTVuGNmRwu8Rp4Ell59RoENcGsDC72OiIgMRt6R3Azcou8puIASn6YsUnkkqoiyXOCYVD7wDDRd34/6irxvERGpQa4gF9WWBB6553UaLtjlLd+1B2AjK4NSRESkKoXvyfkkjBmSGwAmzex89aclIiISLiS7srIpRxERkX4oHeT8fanCOw2IiIgMStA6OTN7s6oTERERqVoVuxAU1mtHcRERkarVEuSwX2hZRESkb/IuBn/FzN7v8f3jyLcFTVzR9iIiIqVkBjmSWwBeJNk2s9e7nj4P4G8K9knUUPNSREQOnjwjuaNwgalX6a7b/vuryL8T9vMAXsvZVkREpLQ8QW4SwLSZXevx3A5c5f1CC8JJJm6BIyIiUpXMIGdmnwDoFeAAF+TKJJEMtNakiIgcTEFb7ZjZDQA3Shx3JqRfERGRPOpaQiAiItJ3fd801S/8bmH//t0f+t2niIgIEBjkSP4QwInYt/bM7Jx/7gSAf4ELcJ/ABblLJJfM7Ech/cb6nwSwCGDcf+sOgDW/x9xA+V0ZVuroW0REegudrlyB2+17E8D5KMB5a3AB7ryZjZvZFNzGq+f8IvIgJNtwBaLbZjZnZnP+fJZJbpAcC+2jwLnMwu0UPrA+RUQkW2iQOw3gtJm9Y2a3om/6JQKTANbN7OFOBWa2B2AOgdv0+FHTtO/74Q7dZrZjZjP+f7dD+ihwLmPQbgwiIkMpNMg97zMsu83BVTW50P2EX5LwSdkOSS7ABbg5HzR7mQPQ8qO9fluBRnAiIkMpNMjdTvj+FACY2UcJz4eU9VoBsBMfwT324i74bQJYINm3Wpn+nuAe3P1GEREZMqFB7rFgRfJFuJFNYhBCyZGPv/c1BmA9R/Nokfpimb5yWjEz7ZAuIjKk+rFO7hxc8EvLMuxVBzOPKGDlqZMZja5mS/aViuQS3KhSRESGVGiQ+wLJp7u+FwWVnvfDSP4A5YPDlH/MMz0YtWlVnWnpX+8lLRcQERluoUHuAoA1kn9F8gWSP4NbNrBuZv8Rb+if/xBuLd2tHq+Vyt9bi4JVZpAzs3ibqcSG5VxCYIaoiIj0X2jtyk9IvgVXwDlK8FiPLQh/Bi4gnIbbsgcAJknCzN4q2F08geROwLFBSE4D+LAriIqIyBAKLuvlsxyfJ/mMXx4Qf+4Tkr1GPEWDFBCWpl/ldOVybC2eiIgMscpqV3YHuNj3C09NVmgPLsAdq+LFSK5A05QiIiOjll0ISF4pcVhUnxIpi8CTBI/konuCaevzRERkuPR9F4IEZdL6QwLVeHaTTG24SioiIjIiKgtyJF+BC0RZAWWiqj4HxZcSWysxgsx6zQUAeO6556p6WRERiQkOciQvAFgqcgjCynqB5FjBgFMm0eVhX3B1MitNNjGzVQCrADA1NRX08xARkd5C95Obh0vE2IQr45VUyzLueQCvleguZBQVcuwK+lsaTERE+iR0JLcIYMbMrhc5yAfHokKWHeQJvo+JCjBrTZyIyGgKza68UzTAeWXKYcVHY0UXd5cdyakAs4jICAsNcqXS6c3sTInDtmL/nZkt2bXFzlZiw+TjJwFMk7SsL+wH3bWu5/pSHFpERPIJna4c2GahZrZHMlrcnWckFz+3wtONZtYhmTcTNNoVYRGxUaqmOUVE6hUa5NZIvmZml4scRPIHZvbdEv1dhUu7zxNco0DYKZv6nzdIkQ93DrqjwCYiMjyCpiv9/TiSLJotuVCyy2gj1Dzp/Of8Y5nqKiIi0gChSwheAPAhgBm/jc4WgG2kZ0KW3hHAzDZJduDulWWtlZuFy4y8mNSA5BqAaQCrSjAREWme0OnK9wE84/+bAE4ie6F36GLwebhA+hYSiiX7aiJR294n4dpEiSFLJK+oLqWISLOEBrloxLZa4JgJAN8q26FPCFkGsEJyo3t3bp9V2YYbna2X7ScPXw0lPvW6SLKj+3IiIsMhNMjtAXi7ROLJg5BOzeyin7Zsk1wHsAEXcKfhRniLvmxW2muskpzD/nRl7lEcyTZ631ecBnAzSkQxM/ZoIyIiAxIa5LZQIj0f5RaDP8KP4Cb8Tt2T/ts7ZnY05bDu1yhVj9LMFqFSXyIiQy8oyJnZ+ZLHlVkMnvRam6ggaIqISPNUumkqyadJvuCzLuPfP07y6Sr7EhERyVJJkCP5il9CsAuX+XipuwmAyyTfqKI/ERGRPIKDnN9PbhNu+cANANfQdZ/OzG6Z2VkAN3x7ERGRvgsKciRPw61VWzazQ2Y25YPZrV7tfYWULb+LuIiISF+FjuQW4PaTe6fr+4mLvc3sGoC5wH5FREQyhQa5oyX3k9P6MRER6bvQIFe2skfutWwiIiJlBY/kSh6Xd582ERGR0kKD3C7Jl4scEMvGFBER6avQsl4X4bIl/9rMfprVmOTbcDsDnAjsV0REJFNoWa8dkm8CuEZyC26EtgVgLFb1pAXgJezv6H3GzP4Q0q+IiEgeoSO5qJr/Hbjtdqawv3wgXqWfcDsWnCmZjSkiIlJYJWW9zGzdzMYBvAlX9YSxrxtwi8XHFeBERGSQgkdycWZ2Ee4+nYiISO0q3YVARERkmCjIiYhIY1U6XQm4bXfgduo+BuAm3G7d71fdj+9rEm6H7nH/rTsA1vxGqv3or+X7m4bLGr0DV/VlDcBVM9vrR78iIlJOZUHOL/JeSnhuFy755EcV9teGy+acN7OO/14LQJvkMoC5KoMOySW4ALcC4ApctugkgHMA2r7fRTNbrapPEREJExzkSD4Dtz5uEvuFl3f8V8t/jQNYJTlnZt+ooM8N/7on44HMzHYAzPjnt1FR+TAf4F4ys+7X2wGwTnIWbjTXJjlhZstV9CsiImGquCcXbZh6DS7oHDKz583sjH88BDfieh/AGZI/CemM5ALcdGHaSG0OQMuP9oL40eE5M0vcHsjM1uH21QOAJR/0RESkZqGbps7DBbhJMztrZjd6tTOzjpnNADgD4GzRepddVuDu83WSGvjgtwlgwQepECtwpchS+eUTUdBdCexTREQqUMWmqSfN7KM8jX1CyFkA58t05kdIYwDWczRf84+LZfqKmUwLqF2i+3EtnxQjIiI1Ct6FIGn0lsRP7ZXdNDUKWDdztI32uis9degDVYvkbs4pyPh5TZXtV0REqlHXpqllj4sCR57jozYtkmMl+4umOscAvJWj/Vbsv7VnnohIzepaDG7ZTR7l761FwSozyPlMy0jZUVV8mnIrsdW++P2/PKNNERHpo9Agtx3bUicXkk/DLaJOev5CwlPxAJJ4fI5jc/OB8ijcfcc89/bi/eQJiiIi0kdBQc7MLgH4W5J/WeCwFbjF00mS7n2VnXIMOtbM9goknrwUOy7vMSIi0idBi8H9KK4NoON3/b6V0jy6r9UGMEX2zD2JFo9XZc/3e6zC1+zJ3/eLArR2YhARGQKhFU86cPfXiP30+TRE+iguTVSfEimLwJOEjALzim8SmzTlKiIiAxQa5KLSXesofp+slwkAryQ8FxKoxrOblOdHcdEC8EprZoqISHmhQW4PwIKZXa7iZACA5P2qXmuALvnHi34dYCZfnmwBAJ577rl+nZeIyIEWml25hfJr3pJkLi4vse6tilFmTz5YzQJYLVKY2cxWzWzKzKaeffbZfp2eiMiBFjSSM7NS5bkyXjNpTVvIFGBfpg99RZQ2XIALLR8mIiIVG6WdwcuMxqIR3+0qTwR4OJq8DmBdAU5EZDhVGuRIvkLyCslfkbztH6+Q/GYFLx8fjRVdZtCPkdw2gM20LXhERKRelewMTvI4XNX/+MapgKsWMgFgluRNuMzD/1uym3gFkcxsya4tdiqtPuI3Zd1RgBMRGW5V7Ax+Am5UA7gsw224qcVoIXYLrhLINNyi8ckygc7M9kjGXzNLPDmlsuSYaCNWvz9eUpsxAONd9TNFRGTAqhjJbQC4YGbvZDUkuQK3Q3jZCiRX4dLu82RXRoGwU9W6NZJLAFppAc6bhhtt5lkgLyIifRK6M/hrcIkXmQEOAHyK/WWSb5TsMtoINSvIAMA5/3ilZF+P8PvJncsR4AA3ctUoTkSkZqGJJzNm9maRA3ygyxMoeh27CVdKbDrHWrlZAHtmllhHkuSa3xB1JamNbzcJYNHMTuY81UloFwIRkdqFTlfuljzuk4A+5+Hu+70FoOfia79AO2rbU2wRNwAskbzSa+cAn8CyBmAlY3fwcbhp1AkA0yrtJSJSv9AgV3jz08DjYGYdkstwQWfDj+4e8kEpWqCdq8RWEj9a3IC7v1eksLS22RERGQKh05VlE0iCdgXwU5AzANokV0hOk5z0iSHbcFOLqQu0zWwVQBQgVxP2fzuLclv/6H6ciMgQCB3JfUjyW2b2z3kP8MkqmfUps/gR3ATJabh7YIBbu3a0wGuk3hv0gVAZkiIiIyo0yK0C2Ca5a2Y/z2pMch7A2wBOBPb7kA92m5kNRUTkwAkt0PwJyfMArpPcglvH1sGji8HH4UZai3BTf2fM7A9BZ33AfXz7Li59sIN3b/wWdz+9hyNPHcarL34J86da+MqxI3WfnojI0KBZ6RyQ/RdxKfbXATyd1AQu6M2Z2fXgDhtmamrKtrbyrTj4+S9/h9d/3MFn9x/g3oP9f7vDh4gnnziE739nEi9/9Yv9OlURkaFAcjtl15qHKinQbGYdfy/suwA+ggtq0dctuFT/EwpwYT6+fRev/7iDP312/5EABwD3Hhj+9Nl9vP7jDj6+fbemMxQRGS6VFGiOKFGjvy59sIPP7j9IbfPZ/Qf4Xxv/hs//t8OazhSRAy94JEfyDZKvVHEyku7dG799bATX7d4Dw7sf/RY/+cWv8cdP78EA/PHTe/jJL36Nb/zTB/j5L383mJMVERkCQSM5kj+EK5i8i/Jr5iSnu5/ey92213TmvQduOvO9750a2hGdkmpEpEpV3JNbh1sWIH125Knw2eXP7j/A5Q9uVXA21fv5L3+Hb/zTBxqFikhlQoPcv5vZ2by7EERIJmVhSopXX/wSDh9idsMU9x4YfnrjN6ltPr59F3/37r/ia3//M5x48//ga3//M/zdu//a14QWJdWISD+EBjmSfKHEcdvZTaTb/KkWnnwifPB998/J0551jabyJtUM6yhURIZT0BXTj+DOl0g8mQjp96D6yrEj+P53JvEXTz7x2IiuyAjvyOd6T3vWOZrKm1STNQoVEYkL3TT1abg1cJMkf0byAslvkXyF5AsJX/MI2IXgoHv5q1/Ee987hW9//Tl8/qnDIIHPP3UY3/76c/ifL2RPZx4+RHzzxS/3fK7O0bdsvQAAABCiSURBVFTepJq0UaiISLfQTIY97AcsApgOfD3J4SvHjuAfXv0a/uHVrz3y/Y9v38W//L//xL0H9xOPffKJQ3jtVO/SoUVGU919hzry1GH8MUegSxqFioj0EnrF2IGrR7kOV68yjwkAWlfXB9F0ZlbZr6RU/DpHU6+++CX85Be/Tg2yaaNQEZFeqhjJLZjZ5SIHkUweakiQaDrz8ge38NMbv8HdP9/Dkc8dxjdf/DJeO3Uida1ZnaOp+VMtXNv+TelRqIhIL6FXqy2U2yA0eD85SZY0nZmlztFU6ChURKSX0OzK82b2fonjMitHy+DlWaLQz9FUWlLNe987pd0VRKSwSrbaqYvf4mcRbs86wN0XXPMbqY5Mn0W22uk3beUjIqMg71Y7lQc5v2ZuEq6W5U0AO2VGezn6aQOYAjBvZh3/vRaAtm8yZ2Z7o9DnMAU5wGVplrmnJyIyKAMPciQvAFhKeHoXwLKZ/aiivjbgsjpP9goq0fNmVtmi8372OWxBTkRk2A1s01SSz5D8EC7AxTdKve4fCTe1t0ryvQr6W4Bbj5c2apoD0PIjr2B19CkiIuGq2IVgE8BJANfgRjmHzOx5MzvjHw/BTfG9D+AMyZ8E9rcCNwXaSWrgA9EmgAU/nRiqjj5FRCRQaFmvebgAN+l3I+i5NMDMOmY2A+AMgLMkXy7Z3yyAMbjF51nW/ONimb7q7FNERKoROpJbgBu9fZSnsc9APAvgfMn+ouBxM0fbaP3ebMm+6uxTREQqEBrkdpNGb0nMbB3uPl0Z0U3GPAvQozYtkmMl+6urTxERqUBokCtT7aTUcf4+VxQ4Mo83s3ibUovP6+hTRESqU0XiSRll1i3EkznyFoPudeyw9ykiIhUJDXLbRXcG93vQJQYMv96ul5Dpv7LH1tGniIhUJLR25SUAf0vyLwsctoL9CiG9VJm0Ea1pO1bhaw5jnyIi0kPQLgR+FNcG0CH5Ntzi7yRjAN7y7afInrknLSRP80W1IlGidFbZUVUdfYqISEVCt9rpwN1fI4DVHO2J9FFcmpCgMZ7dZGj6FBGRitSxM3iaA7NruC8VtuD/948kf5nQ9AsA/mswZzV0Dup71/s+WA7q+wbC3vtX8jSqZWfwNHl2DSc5VnD6MDgAV92nma0ix+iX5NZB3X/voL53ve+D5aC+b2Aw7z00u7LszuBpkhaXh2ybU/bYOvoUEZGKBI3kzKxsea6010yK6mVGY9E9tdslT6eOPkVEpCJ1LQYvIz4yKrrQuoqR3KD67JYnoaepDup71/s+WA7q+wYG8N4r3xm8X3wtyF3/vzO+2HNa+xb2iyqfTNsmZ5j6FBGR6ozMSM4nfUSjozyjqnj6f6n7hnX0KSIi1QnNrhy0q3Bp93nWr0VBqVNiIXdtfZKchNveJ1pndwfAWtYoUmRU+T0bF+H+fsbhPiBeMbOLtZ5Yn/jlQzMAJuHe7xbcmuMLgdeqWpFcAoAy/279vO6NzEjOizYlncnR9px/vDIqfZJsA7gEoG1mc2Y2B1cGbZnkRpO37yHZIrlCcpvkLsmb/j0vNPl9Z/E/g+m6z6MfSE6SvAl3cWvDTfEfBTAP4BjJ7VpPsGL+/W7DrQe+YGYT/v0u+ya7UaAYNf68V1CinGHfr3tmNlJfALbhqqyMZbQzuP3u0tqswd1zWxlUnynHbsDdz+vZR/R83T//Pv2bLvn3vgD36bYFV8N0zf9MDW49Zu3nOuCfy6x/77N1n0sf39tSxu9F6t/mqHz53+tdAK2UNi3fpl33+eZ8Ty3/N3sz9nda6N9rENe92n9QJX9ZUn+Y/gefenGItYm+JvvdZ45j085hzLcZiT+AAu99CW5aIun5WZT8AxrlL//vvdvEIAdgOsff04pvs133+Vb0nm+m/X33+Nlktq3xvSxFH+j9B9Hp2AfS3H+jg7ru1f4DC/whT/d4rpXnh1IkyFXVZ8pr7+b5tOI/1Vjap8FR+vI/t8yLWOxn37gLfsp7bjf1Pef5fY8F+JGfvfAf1HLP8PiAOLQfZtFj1FUyyA3kujdq9+QAPLyxOQOg7e/jTPv57iW4qcVFM1vMeI1VANFNzVXLSPevos9e/E33Mbj6n1mi+4OF+xlSK3D3X1L5n310Q36lr2c0BPxN+D00MEPX338ZQ/a/4zJcMsZyRrtRMAOXXJLXJoChLfNlFSTHDPK6N2rZlQ+Zy7qZ8DflJ/23d8zdyM37GnmSSSrts4foH+5maivfl3+cRTP++CezPlzErMKN6Fokixw3ilbMbMZfCBrDJxAsAA8/ZCaynLVdR0jRoFVFwfthNrDr3sgGuYgPPANNr6+4z+iXP8+n9qhNq0TB6KHiRystkrsA5s0s6xNd/I9hCu5TfuPEstSa6Kx/bOS/XYqbAMZILlm+9PpplN+SbFQM7Lo3ktOVTeErpETpsZn/2GYWbzO00xk5RWsKo810s8SneyaqP536+ZHOS9bcNZHRp/emvr8kUVBf8WvkEvlZohaaNYp9xKCveyM/khtx8SoqRacnitbSHDbxT/N57lfE32+eKY5RdAnNmIZOEk3xP/z3I7kCd8Eb9487cNO1jbkfaWabJDtw779NcgZu9uKREYn/kLMGd39/ZGdpchjodU8juXqFLHIc6QXS/iJ2FG4BcJ4byvFf7iI38UeC/wT/YZMu7nF+ejpyh+QYyTW4LMJFc4uAZ+CSuG5mjXhG0GnsJ0/NArgVX+TvA9x1AMtZ9ysbYKDXPY3kRs8e3D904coCw8Z/Ws17f+al2HFNvKezXDQRasTEP6TswY3W5robmdmqv+C3SWYmqIwKM9sjeQL768rGAGyQXIVLkV+EG9018Xe7CqWvexrJ1Suq01YmLXekR3JF+ItelGnYuHqGfsquydOUQOx3HW7abiOpoU/O2IELdKM+Lf+Qme35DzLx3+EFuMC3coAC3ECvewpy9QoJVOPZTRojPnV1obaz6IPoJvwBuMDFf9cXc2TTRs83MdP0NtwMRvwCv+E/7BwEA73uKcjJUPOjuOiPf66BN+TbaP4orlue7MpopDfblALdvgj5NoA9M4uKUceD/ZIvUN6I9zssFOSGRIlf7KYvFo1c8o8Xc3z6Hyk+uWKtgYG7l/h7LLIAGHD3sEaaTzLZRldiib8vOYf9n88kgAMT6AZx3VOQq1fIxa3xF0YfBGbhyq41arTj/7jnmpJYUVCeDNL4xeylxFYjwGeWbsD9ez82ivUf3k5gf4Tbwn4pqyYa6HVPQa5eZUZj0Sef21WeyLDxF4Y2XIBrSq3OuBU0pwZpHvHf9YOWZLUGYDNtkX+PpJTppu4jiAFf9xTk6hX/Yy+aRdbYkVxszdB6EwNcVIC5qWviEsQTa4oGrZFNsopVMMk1E+FnLKLR/WNLLBpioNc9Bbl6xRc1Z/4hd6VTN25BdMw23Cffpv6RrzRt+jVLV0AvGrRG+cPADFB4bWf0u9GY5RNdBnrdU5CrkU84iD6Z5PmFjn8CHuU//EQkN+B2dmhkgPOjuGmSlvWF/d+Jta7nRnV3guhCfzJH2/jFb5TLuLVQ8G81ViShkUFu0Nc9VTyp31W4dWB5pnCiX4hOEzPy/F5jqVsg+anM8VGd6jOzDsm8Baaji/siYmn3o/re4d7DJPJd2OJtRrmg8x2UC1Y7GP17kWkGdt3TSK5+URZVnpJO5/zjlT6dS238FjOtHKWtpjHiKeVmtpPnK3bInYTvj5po+5g8leSj34POiL/nbeCxKbc8UqvCNMDArnsKcjXzGVcduCmsrE81s3AJC40qbeWn387lrN34Eho6Vdt0Plhtwu2tlvVBJXp+1CvcXPWPuauZ+J/NOBq83c4gr3sKcsNh3j8m7qsWq8o+n9RmFPl7VItmluc+DeA+4TY56abpomzZxIu+/9AzCZddO9IFAPz02iJc5ZbMnRX8iG8NPbbiGXJjXY95DOS6RzMre6xUKLYj9Ez3ehr/i38TDVsz5t/XBtz7Tls7E+01NgFgwcw4gNOrnU8+Adwi4pG+2Mf5ILYGF8Tmup6bhvud2GzSrgyxv+91JAQwf0FfwQhtt+NHYdPYn37cg1v6sJUnSA/iuqcgN0T8H3gb7g9hA+7CPw33SWdkfvHz8H8c2yh+U75TYNQ3kvzPJrrgAW6Kb3HE7009Iva7Pg43pbcHN3qbhvtdb9SUPPDwor0CN/3WgZuRiN73lP//5WEv1h37kJLXybT31O/rnoLcEPL/6NEmkztN+hQf8Z9a25kNH/fYp/+m8NmlmVNaTRrJRksq/P828ne9m/8Q04ILbNFu6JsjNj1ZuX5d9xTkRESksZR4IiIijaUgJyIijaUgJyIijaUgJyIijaUgJyIijaUgJyIijaUgJyIijaUgJyIijaUgJyKl5Ck4LFI3VTwRkcJihZSPHvRyVDLcFOREpDCSawDGmrRTgDSTpitFpIyilehBcoHkhv/a9sWZRfrqcN0nICKjxW+1Auzvep3VfgwuIO5FIz8/3blNcqJJWwjJ8NF0pYgUUnSqkuQ2AHTvA0hyF8DVJm0ELMNH05UiUlTuqUqSG3B7p53u8fQWgLMVnpfIYxTkRCS3IlOVfolBtNN3rwzMPQBjujcn/aQgJyJFnEOOXaz9fbgVuPtwqwnNxvzjeIXnJ/IIJZ6ISBGzAOZytFuBC2LLKW2i4NYKPSmRJEo8EZFc/FTlmpkxo90YgF0ASGtLMrr4LKaM9kSCaLpSRPJaBLCeo91b/jGxbdd9uDshJyWSRkFOpAFILpG8SdIKfE0XeP0xuCSSKzmaRzUt09pOxf5b6+Skb3RPTmTE+TT93AHL20Gx4HIWAMwsdSTnA2eUULKZ0jS+xk5BTvpGQU5khJFswyVuzJjZZtdzN81soqKu5pBvqjKelLJLpt6+A1z2pQo8S98oyImMKD9qmu4VyHySSCUjpNhUZZ6symhEuQqXYdlLC24HAyB9tCcSTEFOZHSt4NFpv7gZAJ2K+sk7VTmG/eUAG0k1KbvuBW70aiNSFSWeiIyu0ynFjc8C+LCifvJOVcaTSfLej8tV5FmkLI3kREZU0r0sn54/hgpGcrGpyjzFmKNRXNZ9tmgk19H9OOk3jeREmmcaACrawuYsXNDKc+8sujeY2G8sAANAO/DcRDIpyIk0zzlUdz9uEfmnFKPglRZco2110mpailRGQU6kQfz04iTcNjZVvVahHcCRXsEkWiieVtNSpDIKciLNEt3v2q7gtRaQf6oScFvnxB8f4bfeAYAdjeJkUBTkRJrlnH8MHsn51yqS/ZiVzRmN3rQTuAyMgpxIs8wCgJkF3ZMj2ULxqcpoxDfW/QTJJbjsy9UCI0ORYApyIg0Rq+xfRVblLIpNVUZLGjp4dL1cFDBX4DZb1ShOBkpBTqQ5oqnKKkZKRacqI8sAJqOA6wPcBoB1M8uz1k6kUloMLtIct+GSPoLWn8WmKueLHmtmmyRnAFwiGWVZLmeVBBPpF+0MLiKP8PfP3jKzo3Wfi0goTVeKSLdFuF0EREaegpyIPOSnKlvItwO4yNBTkBORuEW4rMqqyoKJ1EpBTkTiZqGpSmkQBTkRAaCpSmkmBTkRiSzC1ZXUVKU0hoKciERmkW8HcJGRoXVyIgIAILkL4GRFm62KDAUFORERaSxNV4qISGMpyImISGMpyImISGMpyImISGMpyImISGMpyImISGP9f2rxUB0Q0xfpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fontsize = 40\n",
    "markersize = 100\n",
    "plt.scatter(tau_ratio_all[:12], P_normalized_error_all[:12], s=markersize)\n",
    "plt.ylabel(\"propensity error\", fontsize=fontsize-5)\n",
    "plt.xlabel(r'$\\tau/\\theta$', fontsize=fontsize-5)\n",
    "ytick_pos = [0, 0.2, 0.4, 0.6, 0.8]\n",
    "plt.xticks(fontsize=fontsize)\n",
    "plt.yticks(ytick_pos, fontsize=fontsize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "change gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma_ratio_all = [0.25, 0.5, 1, 2, 3, 4, 5, 8, 10]\n",
    "p2 = mp.Pool(30)\n",
    "result_2 = [p2.apply_async(evaluate_sq, args=[10, True, 1, gamma_ratio]) for gamma_ratio in gamma_ratio_all]\n",
    "p2.close()\n",
    "p2.join()\n",
    "\n",
    "result_2 = [item.get() for item in result_2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p3 = mp.Pool(25)\n",
    "result_2_extracted = [p3.apply_async(extract_optimization_details, args=[item]) for item in result_2]\n",
    "p3.close()\n",
    "p3.join()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_normalized_error_all_2 = []\n",
    "for idx, item in enumerate(result_2_extracted):\n",
    "    print(gamma_ratio_all[idx])\n",
    "#     print(\"iter {}, elapsed time {}, A normalized error {:.4f}, \\\n",
    "#     P normalized error {:.4f}, alpha: {:.4f}, theta: {:.4f}, loss: {:.4f}\".format(*item.get()))\n",
    "    P_normalized_error_all_2.append(item.get()[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcsAAAE9CAYAAABkyUlDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dz3Mc533n8c+XkpxyqUyBkLWHpEomQVXtMRIIHn2wBLD2slWxCiQvvkUE6Jz2oAVE/wMUqBz2ZBugfNnaQyhA6xy9ISgdfCQAKseUTTA5JIetkIRccW0ltvTdw/M00QJnpn88PT3TM+9XFWpATHc/Tw9n+jvP08/zfczdBQAA+js16goAADDuCJYAABQgWAIAUIBgCQBAAYIlAAAFCJYAABR4cdQVSGFm85JWJc3GPz2RtO3uuy2VvyZJ7n6rxr4jrTsAoDzr6jxLM9uUtCDpmrsfxL/NSdqMm1x296Mhlr8maUPSLXdfr7jvSOsOAKimky1LM7sraU7ShXxQcfdDSUvx+X1J5xsud07SoqT1WH6dY4yk7gCA+pLvWZrZ+2b2dhOVKVneikLAGtT6uixpLrbgmihzzcxcIYgtKXSf7tQ4Tut1BwCkS+qGNbOfS1qR9NTdX22sVoPLfCrpibsPbHnFFtqipPOx1ZZS5szJ4GZm25KWVaEbdth1/+53v+tnz54tuzkAQNL+/v6/uvtrg7Zpoht2R9L9Bo5TyMyWJc1I2iqx+bZCwFlV6DatrYn7h23U/ezZs9rb26tXQQCYUmb2T0XbpHbD/tbdr7j7R1V2MrPTNctbjY8PS2ybtciWa5bVtC7XHQCmWmqwNDN7s8Z++zXLW4iPZboms23mzGymZnlN6nLdAWCqJQXL2KK8XmOAT+WRnnEkahY4CgPOiXt9C303bEGX6w4ASLxnGbtT1yVdM7N1SQcK9y+PFCbZ93JRUp1RRfmpGv2OXWbfUehy3QFg6qUO8DnSceAzhUEpw5LSHTnqrswu1x0Apl5qsDxUaPnsqHyL6byktuZlHikEm1amtTSsy3UHgInSRMtyxd0/rrKTmX1Vo6wsh2qdqRyjbp11ue4AMPVSR8PuqdzozpMe1NgnJWjMFm8yVF2uOwBMvaSWpbtfr7kfIzwbElPorUjS66+/PuLaAMBk6uR6ljXmHlYdgTo0Tdfd3bfcfcHdF157bWC2JgBATY0GSzN728zumNlvzOxxfLxjZj9s4PApKedGvdxVl+sOAFOvkSW6zOysQj7TeYUpJJkzCqNfl83socJqG39fs5g6rcOsFfe4ZplN6XLdAWDqJQdLMzun4/R1t+PvT3Q89WFOIRHBoqQDM5uvGTDzLaw5hQQIdfYdhS7XHQCmXhMty7uSbpZJpm5mG5I+U725g/nlNApHiMYUc732HYUu1x0Apl7SPUsze0/STtlVR+K6jx+b2ftVy4rzE7NWVpkUcPmBNEnrWabqct0BAOkDfJbc/YMqO8SAuVSzvE/iY5kRpVlQOmhiPcoGdLnuADDVUoPl05r7fVlzv+34WCbYXo2Pd2qW1bQu1x0AplpqsKyzekjt/dx9V2FwzGKJ+YrLko7c/Va/Dcxs28yexnupQ9V03QEA7UkNlnWTfKekf7sWH2/02yBmtclv22+b5ViXNTObr1iPmROPZTRSdwBAu1KD5X0ze7fKDnFQUJ3csJIkdz9QWENzzcyeWxIsjiTdlLTl7jt1y+nHzGbMbFnHy5FdMbMyrcWR1x0AUE/q1JEtSftm9tTdPy/a2MyuSfpQ0rmUQt39lpkdSNo0sx2F6StPFALYDUmr7r5VcIwtM7sc99mKgaxfvZd1fM/xpJlYvsye5WO40O94TdQdANCu1ETqX5rZdUn3zGxPYcTngb6ZlGBWIbPPqsIoz0vu/rukWuvZPcDzsYWWdaEeuvuZCscoNSo3tvKscMPy5SbXHQDQHnOvO0Ynd5Bwv++epNP9NlEInpfd/V5ygehpYWHB9/bIYQAAVZjZftFqWI0kUnf3g9gq+rGkLxSCY/bzSOE+3TkCJQCgixpJpJ6J99q43wYAmCjJLUsze9/M3m6iMgAAjKOklqWZ/VzSikImn7pzLgEAGGtNdMPuSLrfwHEAABhLqcHyt+7+11V3MrPTTUwfAQCgDan3LM3M3qyx337xJgAAjIekYBnXsbxeY4DP+ZRyAQBoU+oAn9MKcyivmdm6Qvae+woJCJ702e2i6q9WAgBA61LvWR7pOPCZjpOLAwAwMVKD5aFCvtcd9W9JnnReEvMyAQCd0UTLcsXdP66yk5l9lVguAACtSR0Nu6fQuqyq9nqWAAC0LXWJrus19xuY3R0AgHFCblgAAAqQGxYAgALkhgUAoAC5YQEAKEBuWAAACpAbFgCAAuSGBQCgALlhAQAoQG5YAAAKkBt2yv3T49/r9q8P9bcP/kW///c/6uU/eVF/8daf6tr35ySp73Pfe/XlpGOX2X9Y50XZk1f2NJvW173t8zb3+rcPY1KCT9z9s4r77ZHyrnkLCwu+t7dXevvP/+H/6q/+14H+8NXX+uPXx++DF0+ZTplJJn39tT/33EsvnNJPfzSvH/zn/1Tr2GX2T0HZ01X2NJvW173p8zaz/aKYlBQsMV6qBMt/evx7/Zf/8Wv9vz/Ua+R/+6UX9Kv/9v2e3+DKHHvQ/ikoe7rKnmbT+roP47zLBMvk3LAnCjxtZm+enHtpZmfjyFmMidu/PtQfvvq69v5/+OprffzrR7WPPWj/FJQ9XWVPs2l93Ud13o0ESzN728zuK+SI3Zd0++Qmkj42s/ebKA/p/vbBv3yj+6KqP37t+uWDf6597EH7p6Ds6Sp7mk3r6z6q825i1ZGbknYlXVBYp/JTnVjj0t0fufsVSQ/i9hix3//7H9OP8R+9j1H22P32T0HZ01X2NJvW131U550ULM3sHYWkBOvufsrdF2JQ7Nn+dfd7kvZY0mv0Xv6T9Bz6L3+r9zHKHrvf/ikoe7rKnmbT+rqP6rxTW5YrkpZi2ru8vm1kd/9U0uXEcpHoL976U714ymrv/+Ip0w/f+rPaxx60fwrKnq6yp9m0vu6jOu/UYHkmtharqn+VRiOufX9OL71Q/7//pRdO6b3vn6t97EH7p6Ds6Sp7mk3r6z6q804NlofFm/R0JrFcJPreqy/rpz+a17dfeuG5b2kvnjJ964VT+taLp3o+9+2XXtBPfzTfd1h20bGL9k9B2dNV9jSb1td9VOedmpTgjrtf7fH3m+5+Y8B+JCUYgqpJCaQwZ+njXz/SLx/8s37/H3/Uy996UT9868+efSvr91zZDD4p+6eg7Okqe5pN6+ve5HkPPSlBzOBzx90/P/H3vsEyjoY1d/+gdsHoqU6wBIBpVyZYpg4XuqUwuvUv3f2XJSr0oaRrkiarEx0AMNGSgqW7H5rZB5I+NbM9hfmWe5Jmcll85hTWsFyRNCPpkrv/LqVcAADalDwRxd23zOyJpC1JCzqeNrKS28wUVii5VHP0LAAAI9NIujt333H3WUkfKGTxsdzPA4WkBbMESgBAFzWa4sDdbyncxwQAYGI0uuoIAACTiGAJAEABgiUAAAUIlgAAFCBYAgBQgGAJAEABgiUAAAUIlgAAFCBYAgBQgGAJAEABgiUAAAWSgqWZvddURQAAGFepLctNM/tOIzUBAGBMpQZLk3S7iYoAADCumrhnecXMfmNm7zZwLAAAxk4TwXJe0oIkM7NPzOymmZ1u4LgAAIyF1MWfl9z9i/j7p5I+NbNzkn4SHzfd/bPEMgAAGKmklqW73+vxt0fu/oG7X5V0JrY236e1CQDoqtSW5UDunm9t3jKzM5J+7u6fD7NcAACa1FZSAo+PlyXtxgFBtDYBAJ0w1JZlHCF7Q2EQkEk6krQlaVPSjEJr0yVtc28TADCuUjP4PDddxMzOxhGxjyVtS7og6YGky+4+G+9nPnL3B+5+3d1/LOmCmd0xs7Mp9QEAYBhSW5bbkl6QngXOVUmLCq1IKbQiN9z90aCDuPtH8RifmNmau/9jYr0AAGhMarA0M/uNpFmFblWTdKAwZaRyZh93v2JmP5P048R6AQDQmCbuWZ6Pj1sKQfJB4vG+TNwfAIBGNTEadsPdT8X7j7UCpZm9YmY/M7M3Jb3SQJ0AAGhMcrB09xsN1GNF4X7nPUn7DRwPAIDGpHbD/s9BT5rZOwpB8LHCtJG7faaIbMXHg15ZgQAAGKXUYPlfBz0ZA989KXS1Slo0szuSrrn773LbfSnpo8S6AAAwFE2sZ1mKu3/p7p/GnLEbieUCANCa1GDpxZsAANBtbeWGfSZ2xy60XS4AAHWVumcZp3T0e+7PVa47dk7SRYWRr5+Uqh0AAGOg7ACfJYXkAwsKSdHz3a8HFcrLMvysV9gHAICRKhUss9ytkmRmcwqrhryjEDTLprXLpo5M9dQQM5tx96NR1wMAUF7lqSPufihpyczuSpp39+vNV6scM5tXmMc5G//0RGG5r91xLNPMFiXdNbNdhST0u5KeuPuRmc0odFXPKbTk59x9aRjnAACoJmWe5bqku01VpCoz21ToFr7m7gfxb3OSNs1sXWFJsEZbcA2UORMfF+NPdtyT2+0qLJQNABgDtUfDxmBReWWRJsRW7aKkd7KgFet0mGuNNZo2r6UyDyWtuvsSXbUAMD6SMvi4+wdNVaQsM1tRCFoXBgSUy5Kemtmmu6+OWZm3JN1UaKHOKbQ2DxVS/R2m1hUA0LwmluiqzMzedff/XXP3DUmH+dbdSfEe4K6kFTPbaCAINVpmDLhDu68KAGhW60kJou06O5nZskJLbKdCGUkty1GUCQAYL2WTErzda7UQMzur0JVYRdXt87Ig9LDEtlnLbllp8zpHUSYAYIwUBksz25P0VrwX91cnnr4u6b9XLNNUP6dsliavTLdqts1c4tzGUZQJABgjZVqWZxQCXK+Udo/j37dUruUlSW9Ieq/kts/EKRrZ1IvCwOXuh7kpGQuqcY9wFGUCAMZPmWA5L2nR3T/t8dyhpN2qiQnM7FqV7aN89+2ThH3HvUwAwJgpDJZxYeZegVIKwbLOYJ06La6Z4k0a33doZcaMPSsKyeUz9yVt0X0LAOMldZ7lA0kPaux3KaXcCo4UgtarLZVXqsw4b/OCpE13v5X7+6akR2Z2zd3LjL4FALRgVFNH6shysapGy6tuC3EYZa7E462enLcZkxlsSdqOU1YAAGNg6MHSzN41s/fj4+mEQ6V0ic4Wb9JKmYeS1t19q99O7r6u0Drdjl21AIARS+qGNbOfSzqX+9ORu1+Nz52T9HcKA12+VAgUt81szd1/kVJuV8WWZJn1Pz9RaIFuqCDBQezSXZGk119/PbWKAIAeUluWGwoDVHYlXc8CZbStECivu/usuy8oLCB9NSYzqK1Gi6vqSNZRl5mt5rJStKG7b7n7grsvvPbaawlFAgD6SQ2W7yiswvGRuz/K/hinhsxL2nH3ZyuTxPt+l1Uvu03KCNG6+46iTCk3pzOugQkAGKHUYPlGHBF70mWFLD03Tz4Rp6J8WaOsOi21rDX4uMa+oypT+magZb4mAIxYarDsFxAWJMndv+jzfJ10dykBpImWZVtlSt8M0ucTjgMAaEBqsHwu6JnZWwqtq0EDWeqM8tzL/V44ujWmquu170jKNLOHZuYlu1XzZaW0UAEADRjG1JGrCkF0UJaeXnlmB4r3O7PWWplWXj4g11rPsqkyY4DM9r9c9zgAgNFIDZbf7TF3MptMv9lrBzP7mcIo2jo+iY9lWqZZcDpITB/XRJlZt+qRyqUHzAdmkrEDwIilBsubCpPn/9zM3jSz/6Nwod9x93/Mbxifv68wF/NRj2OVkQWapRLbZtNY7tQsq7Eys/mV7n7G3csEv6ysXfLEAsDomXvdpSXjAczmFVpfWWtox92vxOdekXRbYYrJmfi8S7rl7jdqlrevMC3lzKBAYmauEJjPDNhmW9KiQvLyvtNZmigz5n3dPJnirsd2M5Kexn9eKNo+b2Fhwff26t6eBYDpZGb7MRdAX8n3LN39wN3fUAgkp7JAGZ/7UmFOZZaQ4LzCepYfJhSZLe/VN9jGrDb5bftts6zQvboWg/4wy1xXyGBUdO8zm5d6q0qgBAAMT2MDfGJg7PX3Rz1+6syzzI53oBB41nqNLI3BaFOhtdjIyh1NlBlbpOuS7vZLkh5bn8sqaOkCANqV3A1bq1CzOydS49U5xqJCgNpRSA/3RKFL9YYKkpXnjnFXx92wA3OwNljmvELrMZtec6jQhb0Yj7deN8jTDQsA1ZXphk1KpJ4gefmpOFDmfAxgWRfq4aB7lD2OUWbQTtNlHki6EFujiwpB86Gkm3S7AsB4aixYmtnbChf+osn7jWakiQGs1ekVTZTp7ocKa1cCAMZccrA0s5uS1qrsonrp7gAAGInU9SyvKQxa2VW4/1YmNdsbkt5LKRcAgDaltixXJS25+70qO8UgCwBAJ6ROHXlSNVBGpHADAHRGarCsNXrT3S8llgsAQGtSg2WdpbYAAOiU1GC5bWaVB+vElUcAAOiEpGAZ71dajYC5UrwJAADjIXXqyJuS7ktaistv7Una1/H6jb2UWUQZAICxkTp15DNJr8TfTdIFFSccICkBAKBTUoNl1oKskrbtvKR3E8sFAKA1qcHySNKH7v5xlZ3M7OvEcgEAaE3qaNg9hSWmqiIpAQCgM5Jalu5+veZ+JCUAAHRGasvyG8zstJm9GUfJ5v9+1sxON1kWAABtaSRYmtnbcerIU4WpI7dPbiLpYzN7v4nyAABoU3KwjOtZ7ipMG3kg6VOduI/p7o/c/YqkB3F7AAA6IylYmtk7CutZrrv7KXdfiEHxUa/tY8afPTN7O6VcAADalNqyXFFYz/KjE3/vm3TA3T+VdDmxXAAAWpMaLM/UXM/SEssFAKA1qcGyzhxLSTqTWC4AAK1JblnW3O98YrkAALQmNVg+NbMfVNkhN3oWAIBOSM0Ne0thdOtfuvsvizY2sw8lXZN0LrFcAABak5ru7tDMPpD0qZntKbQY9yTN5LL4zEm6qDBydkbSJXf/XUq5AAC0KbVlKXffMrMnCst0Leh42shKbjNTWKHkUs3RswAAjEwj6e7cfcfdZyV9oJDFx3I/DxSSFswSKAEAXZTcssxz91sK9zEBAJgYja46AgDAJCJYAgBQoNFuWCks1yVpXtKrkh5KOnT3z5ouBwCAtjQWLGOygbU+zz1VGOTzi6bKAwCgLU2sZ/lKXPh5TccjYB9JuhcfTdKspC0z+1VqeQAAtK2Je5bZws+fSroQ17V8w90vxcdTCvMvP5N0ycz+poEyAQBoTVI3rJldUwiU8+7+Rb/t3P1A0pKZLUr6OzPbdPfPU8oGAKAtTSz+fGFQoMxz911JVyRdTywXAIDWJK864u4Pquzg7jti8WcAQIeMavHnuvsBANC6USUl8OJNAAAYD6nBcj+3FFcpZnZa0pMBz99MrBMAAI1KCpbuflvST8zsOxV225C0OeD55ZQ6AQDQtNSpI28qBL4DM/tQIQlBPzOSbsTtF8x6jvGZiz8AAIyN1HR3Bwr3H01h8ecipsGtSgAAxk5qsDxUaAnuaMB9yArOS3q7geMAANCY1GB5JGnF3T9uojKSZGZfNXUsAACakDoadk/Nz5mslOQAAIBhS2pZunvjaevcfaHpYwIAkGJUSQkAAOiMRoOlmb1tZnfM7Ddm9jg+3jGzHzZZDgAAbUod4CNJMrOzkrYlzeubSdLPKIxwXTazh5Iuu/vfN1EmAABtSQ6WZnZO0n785+34+xOFkbIzClNLLkpaVEheME/ABAB0SRMty7uSbrr7R0UbmtmGpM8kvdpAuQAAtCLpnqWZvSdpp0yglCR3X5f0sZm9n1IuAABtSh3gs+TuH1TZIQbMpcRyAQBoTWqwfFpzvy8TywUAoDWpwbLuIs4s/gwA6IzUYFl3oM5MYrkAALQmNVjeN7N3q+wQBwWR/xUA0BmpU0e2JO2b2VN3/7xoYzO7JulDSecSywUAoDWpidS/NLPrku6Z2Z6kTxQWhM4nJZhVyOyzqpCg4JK7/y6p1gAAtCg5KYG775rZgqR7kjb6bGYKwfOSu99LLRMAgDY1kkjd3Q/c/YykH0v6QiE4Zj+PJK1LOkegBAB0USOJ1DPuvqVwHxMAgImR3LI0s/fN7O0mKgMAwDhKalma2c8lrShk8iE5OgBgIjXRDbsj6X4DxwEAYCylBsvfuvtfV93JzE4zfQQA0BWp9yzNzN6ssd9+8SYAAIyHpGAZ17G8XmOAz/mUcgEAaFPqAJ/TCnMor5nZukL2nvsKCQie9Nntolh1BADQIan3LI90HPhM0mLi8QAAGDupwfJQId/rjvq3JE86L4l5mQCAzmiiZbni7h9X2cnMvkosFwCA1qSOht1TaF1WxXqWAIDOSF2i63rN/RZSygUAoE2NJlJvm5ll62TOxj89kbTt7rvjXuYo6g4AqKfxYBnnXM4r5Ip9KOnQ3T8bQjmbkhYkXXP3g/i3OUmbcRrLZXc/GscyR1F3AEB9jQVLM7spaa3Pc08lrbv7Lxoq667CKNwL+aDi7oeSluLz+2ow+UFTZY6i7gCANE0s0fWKmd1XCJT5BZ/vxUdT6GrcMrNfNVDeisJ8zkGtr8uS5mILLllTZY6i7gCAdMnBUtKupAuSPlVoLZ1y9zfc/VJ8PKXQ5fiZpEtm9jeJ5W0odO0e9NsgBqJdSSuxezNVU2WOou4AgERJwdLMrikEynl3v+LuPaeEuPuBuy9JuiTpipn9oGZ5y5JmFJIgFNmOj6t1ymq6zFHUHQDQjNSW5YpCa/KLMhvHkZ5XJNWacqLj4PGwxLbZ/M/lmmU1XeYo6g4AaEBqsHzarzXZj7vvKNzHrCObn1kmEUK2zZyZzdQsr8kyR1F3AEADUoNlnew9tfaL9++ywFG4fxxdmqmVBKGpMkdRdwBAc5oY4FNHnSW68oNdyiZt77XvKMocRd0BAA1JDZb7ZvZmlR3iGph9A0acr9lLSndk3X2bKnMUdQcANCQpWLr7bUk/MbPvVNhtQ9KgOYRNDmrJ5jK+2uAx2ypzFHUHAPSQlMEntio3JR2Y2YcKSQj6mZF0I26/YNZzjM+c+nc7ZjlUNWBC/6Cy62iqzFHUHQDQkNR0dwcK9x9N0laJ7U2DW5WDpASN2eJNhlrmKOoOAGhIarA8VGgJ7qj6wJVezkt6u4HjTI2YQm8l/vPfzOwfBmz+XUn/OvxajR3Oe7pM63lL03vuqef9vaINUoPlkaQVd/848TjPmNlXJbaZqdidmRzImyqz6bq7+5bKteplZnvTuJYo5z1dpvW8pek99zbOO3U07J7qz7Xsp1+Sg5Qlq+ru21SZo6g7AKAhSS1Ld6+btm7QMft9O6jTOszuFT6uWZ2myhxF3QEADRlVUoI68i2sqhP1m2hZppQ5irr3Uqq7dgJx3tNlWs9bmt5zH/p5m3udZDrtizlSn8Z/LsWk7IO2n9Nx0vILg5bFGnaZo6g7AKA5nWlZxkExWSurTOssP12j1n3VpsocRd0BAM3pTLCMPomPZeYtZkHpoEYigGGU2XrdzWzezDbNbDv+bJrZYt3jAePMzJbN7K6ZPTSzp2a2b2Zro67XMJnZSvxsZ+d818w2ur5akZmt1f2/G9Z1r2vBMlsUeanEtlfj450xKbPVupvZpqTbkjbd/bK7X1ZINbgeP1Cd/jD1Y2Zz8WKxHy8eD+P5rkzqOZcRX4OJ/KIUL44PFdaM3VS4dXFG0jVJr5rZ/kgrOATxnPcV5qbfdPfz8ZzX4yZPu/pFIdZ7QzVSfQ71uufunfqRtK+QNWimYDtXWG9z0DbbCvcSN9oos8m6F+x/V+GeZ89ysudH/X85hPfGWjzvFUnzCi305fj/7PFnZdT1HMHrshzPfXnUdRniua0VvC8Gfsa79BPf208lzQ3YZi5usznq+pY8p7n4uX2Y+6xW+j8b9nVv5C9SzTfKwBcyvugDLw65bbKf+RbKbOQ4Ba/PSonzmYnbdOKDVPK81yRtD3h+WTU/hF3+if/XTycxWEpaLPF52ojb7I+6vg2e98NBn+8er0/htiM8l7WscRC/1C7mvtyW/py2cd0b+YuV+AIv9nhurswLUiVYNlVmk8cZcPynZb49xW9ZrgHfTrvyE1+3woth7rWfuMAx4Jw3J/Wcy7zXc18UJqInJX7pK93rFAPr2H4pVo9WYM1gOfTrXtfuWUqS3P2Wwr2/zXh/ajH24a8pdHWuuvtqwTG2JGVTOLa8YHpGE2U2eZxezGxZ4dvTTonNs3uotcoaMxsK96cGiq99NmBqY6g1GgNmNq9wvhM3ojrem5pR8f/jusKCD+sF23XFkkLmtLJ2JY1t+jtPG3wpqb3rXmpu2JHxMFfxfBy0MB//fOjhJnfZY5QZbNNomU0ep4fsDfBw4FaxvPi4rO5fSOaLvuzkbCm0MOfMrMp+XbTh7kvxYjIx4iCNFenZl96+vELu5A6pGvyaWORinLVy3etssMzEwDNwkv+4ljmEumcfojItiWybuRrJ3cdGbD3NmdlTSdfcvejbZf4DtaDQ6pg4uRGFk+hKfJzI/7sCDyXNmNla7Ckpsqj6yyJ2RSvXvU52w+J5MetPNiy68E3j7vltxrabpoRsTmq2uHiRfBfW+earM3qx5XXRCzJFdVjWkpjU8xsk+4KwEZfn6yv2XM1p8lrWz7R53et8yxLP5DMDVe12qZqvdpzkWxdl7uXkz7VMt00X3Vb3u9YHyW5dPPv/M7MNhYvmbHw8VOiGnqj7te6+a2YHCq/BppktKfSofKOFFL8wbSuMgehkr1FJrV33aFlOjpQJ952drB8vhmcUJqKXuWmf/4BUGSjRCbE1cX/SgkQmdrtnnpjZjJltK4z4XPUwEX1JYbDcw6LWV0e9o+OBasuSHuUTTsRAeU/SetE93QnQ2nWPluV0O1J4w1TOlDFO4jfnsvevLub2m8R7XutVB651TP7LzpFC6/HyyY3cfSsGjU0zKxwI1CXufmRm53Q8L3FG0l0z21KYGrGq0NqcxPd3E2pd92hZTo7Z7Jca3S6dbVlWES+e2cjQMoMjOiV2RU5y96uUe58rdEXe7VZmlbsAAAiASURBVLdhHABzqBAwu3yr4TnufhS/FOXfxysKAXRjigJla9c9guXkSAl4s8WbTIR8l9zNkdViCLKBDlNwkcy/z1dLjH7Onp/UkcGPFXpV8oHibvziNA1au+4RLDEVYqsyu4BcnsBBD5ua/FblSWVGw2Ytz+VJSqQfFwzYl3Tk7lni+PwXh7W4mMDEnPOoESwnUI0PyKRPWpbCCFFJulWiNdIpcRDL9gR+Aeglf45VJqFL4f5e58XBPPs6MYAn3ru9rOPXaF7S1ATMYV/3CJaTI+VCOdEX2RhMlhXSGk5U6yteIC5P0gCWCsqM+M1fEC/23aoj4mjguwr/58+1rOMXwXM6bnXP6TjF2yRq7bpHsJwcdVqH2Texx01WZJzEi8umQqCchDy4J21oMvL7lpV/n0/jQLZtSbuDEk70GPyzOKlrmarF6x7BcnLkLxxVR/5NZMsyN99sZxIDZZYofVLnVPaRH8BUNfh1eiBbLiNPqd6R2IuS9Tg8N71mQrR23SNYTo78BPvCi8KJofQTNzk/2lf4Fj6pF4qNSetWLnLii0HV4Nf1LxVLUuX5wdn7Y6KmzuS0dt0jWE6IOLgj+6ZU5oOR/1be9YvIc8zsrsJKLhMZKGOrctHMvOhHx++H7RPPdXU1kixYXCixbf4C2vX0hnOq+FnNJeyYyGDZ5nWPDD6T5ROFuYRluqeyN9bBpI2ijGsdDlyCLXbRzna1C9PdD8ysbCL4LEisKjfdoqvnrnAO8yp3ccxv0/XE609UL+gdajLu1/bTynWPluVkyUa9lUl3djU+3hlSXUYiLk01VyLl26I6PpXA3Q/L/OR2edLn712TLTlVZtWI7H1w0PFzlsJthZNdiWUMzHQ0AVq57hEsJ0gcIXeg0D1X9C1rWWFwyMSkfYvdildL5ka9qAnsfp4GMejtKqzrWPSFJ3t+EjI2fRIfS2fnia/PrCZ4ma62rnsEy8lzLT72XdsxtxLDtX7bdE28h7fq7mXuY0nh2/akDmyaBtno5r6BI355mlcYDd35RBSx23BVIRtR4WoqsQW6rR5LeI25mROPZQz9umfuXmc/jLHYFbkhaenkfKz4AXqoCZp3GM/prsI5D5p3la11eF7SirtbC9UbuTjIRwoT2TsfNDIxGG4rBMPLJ55bVHhP7E7aKiy5z/eO+gTCGBg21KFlumKrcFHH3apHClNe9soE+2Ff9wiWEypeLDYVPlB3FYLIosI3r858gIrED9i+qg98OKjQCu2k+NpkF00pdF2uTsC9u2dy7/NZhW7KI4XW5KLC+3xibjPkxYv/hkK34oFCL0l27gvx3+vjnlg/94WnrAuDzmmY1z2C5YSLb55swdzDSWpZSM++QW8Wbvi851ojkyKOBi7sppuklnU2lSb+c+Le5/3EL0RzCgFyRuE+/G7Hul0bN4zrHsESAIACDPABAKAAwRIAgAIESwAAChAsAQAoQLAEAKAAwRIAgAIESwAAChAsAQAoQLAEMBJlkoED44IMPgBal0t0fmbaU7OhGwiWAFpnZtuSZiZtRRBMLrphAYxC1dUmgJEiWAJoVVyWSQpLagGdQLAE0LarYhkpdAzBEkDb6IJF5xAsAbSGLlh0FcESQJvogkUnvTjqCgCYKsuSLo+6EkBVtCwBtCLrgnX3nVHXBaiKpAQAWmFmdyUduTstS3QOLUtgypnZnJltmtlDM3Mz24/p6Abts1glt6uZzUhalHSnRv0WzWw71iv72TSzuarHAuoiWAJTLHaNPpS0IikLPvOS7prZ2oBd1919q0JRV6RqXbBmNhNbo9uS7rj7hewnbvLQzOb77VuhbkAhumGBKZVLZi5JO5Lux98vKgzEkaTz7n54Yr95SavuvlqhrEpdsLHVuC/piaSlk3WI2zyUNOvuZ078fTn+vUowBwYiWAJTKLa8HinMd1w/OZUjF6x2Twa4mAR9vVcAG1DWU0mXy7Qsc3WbkXTB3Q/6bLchaS3W5Vbu73djWUxPQWPohgWm04akT9x9tVdQiYHwgo5bmJKeBdGZsoEyqtoFe08hUN7qFyijrCX8rIUbA+0sgRJNo2UJTJkYUPbd/XyJbTcl3c0CXfz3trvvViivdBds7ELNUuENXOsydgfvS5K7W/zbmqRDpqegabQsgekzJ2m95LbbkvJrTi5UDJTZKNjNkrtsxMcyWX6ePZ8b0HOVQIlhIFgCU8bdD8oGlBgYF6Rnrb6yQS9zRaFVWRhgYxdvNiK3TKL1J7nfZ2Orsmr9gFJIdwegSBaUVt19aeCWz1tV+aTp+bmdpVuv0ZzCqNmq9QNKoWUJoFC8PzhosE2vfWYU5myWXY7r2T3UigOIpNB9S2YgDA3BEkCRI0k3VD37zopKdsGeUDZQzuZ+32QELIaJYAmgyBNJ8wXTOHq5qmrrVj6Mj2WDXv6eZdWADFRCsARQZFYhw09pcbBOlS5Y6TjglU1VdyP3O+ntMFQESwBFZlS9C3ZZFbtg433KHUlz/XK+St/IGftQx/dRZ3PPb/TcEUjAaFgAReZa6ILNXFMY2bptZhfy9yFja3U9Pr/q7odm9kSh9TovaTdOH3lco1xgIDL4AOgrtvBu51b6KLPPnEKrr29e1xLHWFOYdnKk43uTBwqrjxz02PaGwsCgXXcvm3ABKI1gCaCvGIguVlmwOQteJ1cDAbqMe5YABrmq42W8ylqVxPJYmCgESwCDzEvaK7txLmVd1QFBwFgjWALoKS4OrYr3HVcVRsHWulcJjCuCJYB+5lV9sv+y6ILFBCJYAujnoiokFaALFpOMYAmgnxlVmyu5qrDwMl2wmDgkJQDQU43lrpZVMS0e0BXMswTQCDN7qpCIoOryWsDYI1gCAFCAe5YAABQgWAIAUIBgCQBAAYIlAAAFCJYAABQgWAIAUOD/A3UfkKAN1pqTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fontsize = 40\n",
    "markersize = 100\n",
    "plt.scatter(gamma_ratio_all, P_normalized_error_all_2, s=markersize)\n",
    "plt.ylabel(\"propensity error\", fontsize=fontsize-5)\n",
    "plt.xlabel(r'$\\gamma / \\alpha$', fontsize=fontsize-5)\n",
    "# xtick_pos = [0.2, 0.4, 0.6]\n",
    "ytick_pos = [0, 0.05, 0.10]\n",
    "plt.xticks(fontsize=fontsize)\n",
    "plt.yticks(ytick_pos, fontsize=fontsize)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}