MCM 2021 Scipy and Pytorch demo

mcm2021_scipy_pytorch_demo.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "mcm2021_scipy_demo.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/tupui/4310485c9920868c404728cdf8cfe93f/mcm2021_scipy_demo.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AgylrjFWsC3P"
      },
      "source": [
        "# QMC in SciPy\n",
        "\n",
        "Tutorial from Pamphile T. Roy and Max Balandat"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZFifc6a33gql"
      },
      "source": [
        "Since SciPy 1.7, we introduced a QMC submodule in `scipy.stats.qmc`.\n",
        "Here are some things we can do with it!\n",
        "\n",
        "A tutorial can be found in the documentation at:\n",
        "https://scipy.github.io/devdocs/tutorial/stats.html#quasi-monte-carlo\n",
        "\n",
        "Let's check our version of SciPy and play with the QMC submodule."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "czfyTmTVr0gJ",
        "outputId": "45234c53-bce6-44c8-e88c-e91b8e1ef3f9"
      },
      "source": [
        "# ensure the most recent version of SciPy is installed\n",
        "!pip install scipy==1.7.1"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: scipy==1.7.1 in /usr/local/lib/python3.7/dist-packages (1.7.1)\n",
            "Requirement already satisfied: numpy<1.23.0,>=1.16.5 in /usr/local/lib/python3.7/dist-packages (from scipy==1.7.1) (1.19.5)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ogrNY4q5sPIk"
      },
      "source": [
        "from scipy.stats import qmc"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5bqzmTYkwJzh"
      },
      "source": [
        "First let's define a `rng` object to make the following reproducible. Please don't use this in production code. See the section bellow."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zagi1rpyqf7P"
      },
      "source": [
        "import numpy as np\n",
        "\n",
        "# seed = np.random.SeedSequence().entropy\n",
        "# print(seed)\n",
        "seed = 292114020772849599029278515437886320941\n",
        "rng = np.random.default_rng(seed)"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zzhh2ede1fjL"
      },
      "source": [
        "## Basics"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GwCVF0rOwhk6"
      },
      "source": [
        "There are multiple `QMCEngine` that can be used to sample points in $\\mathcal{U} \\sim [0, 1)^d$."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "d0zhYGo5sk1z"
      },
      "source": [
        "n, d = 64, 2\n",
        "lhs_engine = qmc.LatinHypercube(d=d, centered=True, seed=rng)  # centered to make nice bins\n",
        "sample = lhs_engine.random(n=n)"
      ],
      "execution_count": 14,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "id": "CUU5orYwtea9",
        "outputId": "f805ffbe-64a6-4a4a-dfea-00277e7aa2bb"
      },
      "source": [
        "import seaborn as sns\n",
        "import pandas as pd\n",
        "sns.pairplot(pd.DataFrame(sample), diag_kind=\"hist\", corner=True, diag_kws={\"bins\": 4})"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<seaborn.axisgrid.PairGrid at 0x7fc227a38d10>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 15
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAAFlCAYAAADcR5KFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYk0lEQVR4nO3dfYxcV3nH8d+zm01s8kLAdhLwSxy3TlNDK0hXToCmjQK0TirFf1Cwjao2KGDUNoQqFVL6okBTVYKiUjXFLSyQBvIHTkrVym2N04oWQUts2S0hYAcqYyBeJ+C1gRCwnNjep3/sbDw7nt2583LmPufe70eKtDtzvXu4jH56zrnnOWvuLgBA+UbKHgAAYAaBDABBEMgAEASBDABBEMgAEASBDABBnFf2AHrAPj0MgpU9AKAVFTIABEEgA0AQBDIABEEgA0AQBDIABEEgA0AQOW57m9fylav01OThUscwOnaBzpx6jjEEGIMkvXzFSh05/GTZwwAKsQyP35x3wGamTR/90jDHco6H3vlaxhBkDLPjmOczzj5khMOSBQAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEQSADQBAEMgAEYe5e9hi6Yma7JC1t89ZSSceGPJyouBdztbsfx9x9QxmDAeaTXSDPx8z2uft42eOIgHsxF/cDuWDJAgCCIJABIIgqBfJE2QMIhHsxF/cDWajMGjIA5K5KFTIAZI1ABoAgzit7AN3asGGD79q1q+xhIH/WxbWs62EQOn7msquQjx2j3wFANSULZDO738yOmtnX5nnfzOw+MztoZo+b2bWpxgIAOUhZIT8gaaHW1JslrW38t1XS3yYcCwCElyyQ3f0Lkr6/wCUbJX3KZ+yWdKmZvSzVeHCu6WnXoakf69FvHtOhqR9repqlUqBMZT7UWy7pcNP3k43Xni5nOPUyPe3atf+7uuvhx3Ty1LQWjY3oQ295lTa84gqNjHTzvAvAoGTxUM/MtprZPjPbNzU1VfZwBqbMCvXbx3/yQhhL0slT07rr4cf07eM/GdoYAMxVZoV8RNLKpu9XNF47h7tPqNH+Oj4+Xol5ddkV6vd+dPKFMJ518tS0jj57UmuWXZT891fF8pWr9NTk4c4XJjQ6doHOnHqOMQQYgyS9fMVKHTn8ZE//tsxA3iHpDjPbLuk6Sc+4e22WK+arUK+584ahBOLllyzSorGROaG8aGxEl128KPnv7tX0tOvbx3+i7/3opC6/ZJFWL7mw9OWVpyYPa9NHv1TqGB5652sZQ5AxzI6jV8kC2cw+LelGSUvNbFLSeyWNSZK7f0TSTkm3SDoo6YSkt6UaS0RlV6irl1yoD73lVedU6KuXXJj8d/ei7BkFMAzJAtndt3R43yX9bqrfH13ZFerIiGnDK67QNXfeoKPPntRlF8eoOOdT9owCGIYsHuoNw7AfsM1WqIvGZv4vKKNCHRkxrVl2ka5fs1Rrll0UNoylhWcUQFVkd5ZFCmVMh3OrUMtW9owCGAYqZJW3BSynCrVVHWcUQGpUyCr/AVtumFEAaVAh6+x0uNmVSxZr8dgobcVtMKMA0iCQde50+Moli/Wum9Zq08RubfnYHt1y3xe1a/93CeUGHrABabBkoXOnw4vHRrVpYveCW6wiNikMCw/YgDSokBuap8Mnnj+zYAU4u4Z6y31frGUFzQM2IA0q5DY6VYA5NSmkqOR5wAakQSC30amtOJddGSl3Q8zOKCL97wVyRyC30akCzGUNNadKHkBN15CLNDUstMUqlzVUdkMAealdhTyIaXwua6i5VPIAZtSuQh5UU0MvTQq0GwNYSO0q5LIeyNFuDKCT2lXI7dqkBz2Nb1cJ024MoJPaVcip/1LGfJXwS140lsVWOQDlqV0g9zuN79RoMV8l/NDW1/CADcCCahfIUu9NDUXWgedboz515kxWf8MOwPAlDWQz2yDprySNSvq4u7+/5f1Vkj4p6dLGNXe7+85+f2+qg3+KNFrMt9XspRdeoGtXvZQHbADmleyhnpmNStom6WZJ6yRtMbN1LZf9saSH3f3VkjZL+pt+f2/Kg3+KNFostNWMB2wAFpKyQl4v6aC7H5IkM9suaaOkA03XuKRLGl+/WNJT/f7SlO3CRRot2GrWvTofZQo0SxnIyyUdbvp+UtJ1Lde8T9K/mdm7JF0o6Q39/tKU+4yL7tDg4J3iytifDURV9kO9LZIecPe/MLPXSHrQzF7p7nMS1cy2StoqSatWrVrwB6ZsF6b6HTwOQALOStkYckTSyqbvVzRea3a7pIclyd0flbRI0tLWH+TuE+4+7u7jy5YtW/CXpm4XZh24N/O1jXMAEnBWygp5r6S1ZnaVZoJ4s6S3tlzzpKTXS3rAzH5WM4E81c8vpYqNZ6FlCQ5AAs5KViG7+2lJd0h6RNITmtlNsd/M7jWzWxuX/b6kd5jZVyR9WtJt7t73dgiq2OEfZLSQhdrGOQAJOCvpGnJjT/HOltfuafr6gKTXpRxDHUV7UNbpQSszGmBG7Q4XqoNUBxn1WnV3OtCJGQ0wg0CuoBQPyvppuGFZAiim7G1vA0FjwVwpHpT1sz2NB61AMdkHcrT10ghSHDHab8MNzTJAZ9kHcl0bCxaaFaSoSNmeBqSXfSCX9SeZylRkVjDoijT1wf4AKhDIdazcypgVsA4MpJf9LosqPMHvdjtZWe3GbE8D0sq+Qs69cuvloWQdZwVAHWRfIUvFK7dI7cSzemniqMKsAMC5sq+Qi4q6Pa6Xh5K5zwoAtFeJCrmIVO3E/erUVjyfQa7nRpw5AHVUm0COeu5u2csPKf8GIYDu1GbJot8HYanas8tefqhrYw0QUW0CuZ/GhtTrz2W2FdexsQaIqjaB3E8lWuUqsqwtdBwIBZyrNoEs9V6JVrmKLKMlOuqOF6BstQrkXlW5EaOMNewqzziAftRml0U/yt4JkdrszGH96iWSpD3fOp50+1vUHS9A2ZJWyGa2QdJfSRqV9HF3f3+ba94i6X2SXNJX3L31L1OXruydEMMwzGWEKs84gH4kq5DNbFTSNkk3S1onaYuZrWu5Zq2kP5D0Ond/haTfSzGWQTQ+VP1gnWE2zlR9xgH0KmWFvF7SQXc/JElmtl3SRkkHmq55h6Rt7v4DSXL3o4MeBA+Qihnmg8s6zDiAXqRcQ14u6XDT95ON15pdLelqM/tvM9vdWOIYqDJbpnNqSe61hbtXVZ9xAL0o+6HeeZLWSrpR0hZJHzOzS1svMrOtZrbPzPZNTU119QvKeoCUW0syywhA+VIuWRyRtLLp+xWN15pNStrj7qckfcvM/k8zAb23+SJ3n5A0IUnj4+NdJVpZD5By29qVahmBBhCguJQV8l5Ja83sKjM7X9JmSTtarvknzVTHMrOlmlnCODTIQZRV+eW4tWvQywi5zRKAsiWrkN39tJndIekRzWx7u9/d95vZvZL2ufuOxnu/YmYHJJ2R9B53Pz7IcZT1AKnuW7ump11fPfLDrGYJQNmS7kN2952Sdra8dk/T1y7prsZ/yZRxeE+d/0rzbGX89e/+qLIt50AKtE4nUuetXbPr52+/YU2tZwlAt8reZVFpw25JjmJ2/fwf/mdSd960lp0bQEFUyInVsTFldv386WdO6sHd39Htv7hGoyPS66+5TD+3/NLK/u8G+kWFnFjUv+WXUvPOlqefOalP/NchXXPFJUMJ45yacYBWVMiJVfks5fmUtX5ex9kIqoUKuU+dKrJhtyRHUUZrdB1nI6gWArkPRRofaEkenhybcYBmLFn0oUh7dBW2v+XS/lz3Zhzkr+cK2czeNsiB5KhoRZbzyWY5tT8zG0Hu+qmQ/0TS3w1qIDmqQ0UW+ZCkdpV77rMR1NuCgWxmj8/3lqTLBz+cvNShPTrqLpGFdlQMu00eGJROFfLlkn5V0g9aXjdJX0oyooxUYX24kzJmAUXWrCNX7kCvOgXyv0i6yN0fa33DzD6fZESZKePgomEa9iyg6F7iqJU70I8FA9ndb1/gvXB/HRqDN+xZQNHKtw7r96gf9iGXILf23mHuEim6c4UdFagi9iEPGe29Cyta+dZh/R71Q4U8ZKnae3OruufTTeWb8/5uoB0q5CFL8TCqSlU3lS/qjAp5yFIcNlSVQ3Vmq/w935r5s4rrVy+h8kWtJA1kM9tgZt8ws4NmdvcC173JzNzMxlOOJ4IUD6OqcKhOTi3aQCrJlizMbFTSNklvlDQpaa+Z7XD3Ay3XXSzp3ZL2pBpLJCmm5EUfhEU+JIhGDyBthbxe0kF3P+Tuz0vaLmljm+v+VNIHJOVTzvVp0A+jilTd0SvQKlT5QL9SPtRbLulw0/eTkq5rvsDMrpW00t3/1czek3AslVak6o5egfbT6BG58ge6UdouCzMbkfQhSbcVuHarpK2StGrVqrQDy1SnFu7orca9tmhXaYcJkDKQj0ha2fT9isZrsy6W9EpJnzczSbpC0g4zu9Xd9zX/IHefkDQhSePj4zHm2JmJ3mrc69p69Mof6EbKNeS9ktaa2VVmdr6kzZJ2zL7p7s+4+1J3X+3uqyXtlnROGNfdoBo+cmg17mVtnbVnVEmyCtndT5vZHZIekTQq6X53329m90ra5+47Fv4JGOR0vKoNF9Erf6AbSfchu/tOd7/a3X/K3f+s8do97cLY3W+kOp5r0A0fKVuNy2rdzqHyB4qidTqw6A/iZpX5YK2qlT/qqdat09EP5EnRZp1C2a3bHDKEqqhtIEdvlJDymY7zYA0YjNouWeSwXSqX6XiqB2s0fKBuahvIuazP5vA3+1L83T0aPlBHtQ3k+aq6Ky5ZpENTP6Yq60KKSj6HGQwwaLUN5HZV3Yff+modePpZqrIeDLqSz2UGAwxSbQO5XVXnLv3aX3+RqiwAGj5QR7XdZSGdu13q6LP57xaIvpWvqFx2mACDVNsKuZ3cq7IqPQjLZYcJMEi1rpBblV2V9Vvdlt2gMWg0fKBuqJCblFmVDaK65UEYkDcq5BZlVWWDqG5zabUehqqspaNeCOQgBtF+XPaSSxQ5tMUD7bBkEcR8DxQXj41qetoLVeo8CJtBUwlyRYUcRLvq9s6b1urO7V/uqrrjQRiHHSFfVMhBzFa3y7der899/ajOTEsP7v6Onn7mZOjqLuIBQLlvX0R9EciBjIyYTjx/Rvd97uCc16PulIi67znFYUfAMBDIweRU3ZW1VtupKmctHblKuoZsZhvM7BtmdtDM7m7z/l1mdsDMHjezz5nZlSnHk4OcdkqUsVZbdAcFa+nIUbIK2cxGJW2T9EZJk5L2mtkOdz/QdNmXJY27+wkz+21Jfy5pU6ox5SCn6q6Map4dFKiylBXyekkH3f2Quz8vabukjc0XuPt/uvuJxre7Ja1IOJ6elNFgkEt1V0Y1zw4KVFnKNeTlkg43fT8p6boFrr9d0mcTjqdrUR9aRVFGNZ/TGjvQrRD7kM3sNySNS/rgPO9vNbN9ZrZvampqaOOKdFhP1FbgYVfzOa2xA91KWSEfkbSy6fsVjdfmMLM3SPojSb/s7s+1+0HuPiFpQpLGx8eHlkRRDuuhUj8rpzV2oFspK+S9ktaa2VVmdr6kzZJ2NF9gZq+W9FFJt7r70YRj6UmUw3oiVer9GkSln8saO9CtZIHs7qcl3SHpEUlPSHrY3feb2b1mdmvjsg9KukjS35vZY2a2Y54fV4oo0+OqPMji0B9gYUkbQ9x9p6SdLa/d0/T1G1L+/n7beqNMj6vyIIsta8DCKtupN6h110H/NeVeVKUVOMqaPBBVZQO5StVYlEq9X71W+hEPMAJSqGwgV60ai1Cp96uXSp8dJqiTygZyVdZdq6SXSr9KMx2gkxCNISlE2SExaFEbRIrqdstaVXaYAEVUtkKuyrprszpO35npoE4qVSG3Vo+SKtVAUKUGkaKqOtMB2qlMhVyH6rFqDyqLqOJMB5hPZSrkOlSP7Vq5r1yyWIvHRrNdUy6CVmnURWUCuQ4Pf1qn71cuWax33bRWmyZ204oMVEBllizq8PCndfq+eGxUmyZ2D3xLGI0YQDkqUyHX5eFP8/T9xPNnBj4r4AAgoDyVqZDr+PAnxayg30YMqmugd5UJZKka7cXdSHHoUD87Oeqw0wVIqVKBXDcpZgX9VN20OQP9IZAzN+hZQXPV/ZIXna83j6/Q1ZddLPeZCnihsK/jPmlgkAhkzDFbda979w363yd/qD/8x68WXn6ow04XIKXK7LIoW+6H/jQbGTFNu14IY6lYo01ddroAqVAhD0AVH2b1svxQx50uwCAlrZDNbIOZfcPMDprZ3W3ev8DMHmq8v8fMVqccTyq5tW0XqeZ7/YvbtDkDvUsWyGY2KmmbpJslrZO0xczWtVx2u6QfuPtPS/pLSR9INZ6UcmrbLtr4wfIDMHwplyzWSzro7ockycy2S9oo6UDTNRslva/x9WckfdjMzN2zWoDN6WFW0a1pLD8Aw5dyyWK5pMNN3082Xmt7jbuflvSMpCUJx5RETtVkN9U8yw/AcGXxUM/MtkraKkmrVq0qeTTnyqmazKmaB+omZYV8RNLKpu9XNF5re42ZnSfpxZKOt/4gd59w93F3H1+2bFmi4fYnl2oyp2oeqJuUFfJeSWvN7CrNBO9mSW9tuWaHpN+S9KikX5f0H7mtH+cmp2oeqJtkgezup83sDkmPSBqVdL+77zezeyXtc/cdkj4h6UEzOyjp+5oJbSRWt0OYgFwkXUN2952Sdra8dk/T1yclvTnlGAAgF7ROA0AQltuSrZlNSfpOm7eWSjo25OFExb2Yq939OObuG4r8YzPb1fgZRX5uXXEv5urpM5ddIM/HzPa5+3jZ44iAezFXqvvBfT6LezFXr/eDJQsACIJABoAgqhTIE2UPIBDuxVyp7gf3+SzuxVw93Y/KrCEDQO6qVCEDQNayCuS6HHhfVIH7cZuZTZnZY43/3l7GOIfBzO43s6Nm9rV53jczu69xrx43s2sL/lw+c034zJ2V5DPn7ln8p5n2629KWiPpfElfkbSu5ZrfkfSRxtebJT1U9rhLvh+3Sfpw2WMd0v34JUnXSvraPO/fIumzkkzS9ZL2DOge85mbew2fubPvd/2Zy6lCfuHAe3d/XtLsgffNNkr6ZOPrz0h6vZlV9dScIvejNtz9C5o5D2U+GyV9ymfslnSpmb2sw4/lMzcXn7kmKT5zOQVybQ68L6jI/ZCkNzWmS58xs5Vt3q+Lover23/DZ+5cfOZmdP2ZyymQ0b1/lrTa3X9e0r/rbCUHpMJnrg85BfLADryviI73w92Pu/tzjW8/LukXhjS2iIp8fnr5N3zmmvCZm6Prz1xOgfzCgfdmdr5mHqDsaLlm9sB7qfoH3ne8Hy3rVbdKemKI44tmh6TfbDz5vl7SM+7+dId/w2duLj5z3en6M5fF39STOPC+VcH7caeZ3SrptGbux22lDTgxM/u0pBslLTWzSUnvlTQmSe7+Ec2cy32LpIOSTkh6W6efyWduLj5zc6X4zNGpBwBB5LRkAQCVRiADQBAEMgAEQSADQBAEMgAEQSBXTKfTuIBB6nTiGbpDIFeImY1K2ibpZknrJG0xs3XljgoV94CkQn+9G50RyNXCaVwYqgInnqELBHK19HKiGYAgCGQACIJArpZeTjQDEASBXC1FTicDEBSBXCGNv1gxexrXE5Iedvf95Y4KVdY48exRST9jZpNmdnvZY8oZp70BQBBUyAAQBIEMAEEQyAAQBIEMAEEQyAAQBIEMAEEQyAAQBIEMAEH8P5imN0gx71ETAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 360x360 with 5 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_jVyie7LxL8r"
      },
      "source": [
        "With the samples, you can further compute quality metrics such as the discrepancy and also scale to bounds."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "VmpnCrwstedZ",
        "outputId": "a19e2287-ed48-442c-85b8-9a6b9eb1a239"
      },
      "source": [
        "qmc.discrepancy(sample)"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.0003283487484500025"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 8
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "id": "WmYpgzuExg4o",
        "outputId": "8ffbc5d2-05ca-473a-c3bf-a86b09e1dcf4"
      },
      "source": [
        "scaled_sample = qmc.scale(sample, l_bounds=[-1, 10], u_bounds=[4, 22])\n",
        "sns.pairplot(pd.DataFrame(scaled_sample), diag_kind=\"hist\", corner=True, diag_kws={\"bins\": 4})"
      ],
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<seaborn.axisgrid.PairGrid at 0x7fc227837850>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 16
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWUAAAFlCAYAAAAzhfm7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYs0lEQVR4nO3de5BedX3H8c9342oil6oJBM2FNQ5KI7UwXRkGZcQwY1Okxc600Dg6qLShVgWUDlWcKfWPzjjV4qV2qqukoZZB6IiXcRyUSbXgSGIjoohQLxggEMgmXogwEZL99o/nWbLZ7HM/v3O+v3PerxmG5NnL+T3u8vH7u3zPMXcXACCGsaoHAAA4hFAGgEAIZQAIhFAGgEAIZQAIhFAGgECeVfUA+sS5PaRiVQ8AmItKGQACIZQBIBBCGQACIZQBIBBCGQACIZQBIJBcjsR1tGLVaj2y86Gqh6FF48/Rwad/W/UwJMUZS5RxvGjlKj380INVDwPoi2Vy686OgzQzXfipb5c5lgXdeMmZIcYhxRlLpHF0+T3nnDJCYfkCAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIhlAEgEEIZAAIxd696DD2Z2S2SllU9DrXGsKfqQSTUxPe3x93XVzEYYCFZhHIUZrbd3SerHkcqvD+geixfAEAghDIABEIoD2aq6gEkxvsDKsaaMgAEQqUMAIEQygAQyLOqHkA/1q9f77fcckvVw0A9WZePsbaHVDr+3mVRKe/ZU+d+BgA4JFkom9kqM/uGmf3IzO4xs8var3/IzO4zsx+Y2RfM7HmpxgAAuUlZKR+QdIW7r5V0hqR3mNlaSbdKOsXdXyHpx5Lel3AMAJCVZGvK7r5L0q72n/eZ2b2SVrj71+d82lZJf5ZqDOhuZsa1Y+8Teuzx/Vp+7GJNLD1KY2PdllgBpFbKRp+ZTUg6TdK2eR96m6QbyxgDDjcz47rlnkf1npvu0v6nZ7R4fEzXXHCq1r/8BIIZqFDyjT4zO1rS5yVd7u6Pz3n9/WotcVzf4es2mtl2M9s+PT2depilm5lx3T/9G93xsz26f/o3mpkpd6N/x94nnglkSdr/9Izec9Nd2rH3iVLHAeBwSStlMxtXK5Cvd/eb57z+FknnSTrHO7QUuvuU2m2xk5OTtTqaFKFKfezx/c8E8qz9T89o9779WnPc0aWMIVcrVq3WIzsfqnoYkqRF48/Rwad/W/UwwoxDijOWF61cpYcfenDgr0sWymZmkq6VdK+7XzPn9fWSrpT0Gnd/MtX1I+tUpZ586VmlBeLyYxdr8fjYYcG8eHxMxx+zuJTrj6LqtfBHdj6kCz/17dKu182Nl5wZYixRxiHFGcuNl5w51NelXL54laQ3S1pnZne1/zlX0ickHSPp1vZrn0w4hpC6VallmVh6lK654FQtHm/9CsxW6xNLjyptDMOYnWWc+/HbteHT23Tux2/XLfc8WvryD5BKytMX39LCXStfTXXNXESoUsfGTOtffoJOvvQs7d63X8cfk8fpiwizDCClLDr6Uqlqsy1KlTo2Zlpz3NE6Y80yrTnu6PCBLMWYZQApZXHvixSq3GzLtUqNIMIsA0ipsZVy1UfCcqxS52r6LANIpbGVMkfChscsA0insZXy7DR4rsXjYzrh2MWVNnXkgFkGkE5jQ3mhafAn3niafrRrH8etemCzDUinscsXC02D3aXX/8vtXY9bVd24EAGbbUA6ja2UpSOnwbv3da8AaVxoYbMNSKexlfJCelWAOTYupKjs2WwD0iGU55itAOefKpitAHM7sZHylMTsLCPi+wZyRijP0asCzG0tNcfKHmi6Rq0p99Pw0O24VW5rqZySAPLTmEq5iKl8bmupuVX2ABpUKRfV8DBs40IVbcm5VfYA0t7kfpWk/5C0XJJLmnL3j5nZC9R6Lt+EpB2SLnD3X6Yax6wqN+mqakvOrbIHkLZSPiDpCndfK+kMSe8ws7WS3itpi7ufJGlL++/JdWqrLnIq36karrItmZZkIC/JQtndd7n7ne0/75N0r6QVks6XdF37066T9IZUY5gr9VS+W2MJG24A+lXKRp+ZTUg6TdI2ScvdfVf7Q4+qtbyR3KhT+V5NGN2On7HhBqBfyUPZzI5W64nWl7v7463nqba4u5vZgjteZrZR0kZJWr16dSFjGbbhoZ814W7V8OkTS7s2pQDArKShbGbjagXy9e5+c/vlx8zshe6+y8xeKGn3Ql/r7lOSpiRpcnKy51GFlDcK6qcJo1s1zIYbgH4lW1O2Vkl8raR73f2aOR/6sqSL2n++SNKXRr1W6hsF9bMm3GvNmg03AP1IWSm/StKbJd1tZne1X7tK0gcl3WRmF0t6QNIFo14odTtxP2vCVMOj4ZaoQEuyUHb3b0nq9F/VOUVeK/UZ5F43KprFTXqGU+XjpYBoatFmnfp0A1VwWtw4CTikFm3WZbQTsyY8uk7NNZzjBg6pRaVMJRtftyUKznEDh9SiUpaoZOer4gZI3XRrNefGScAhtaiUcbiIG2e9NmOZ6QAttamUcUjKGyANW4H3uiEUMx2ghVCuoVQbZ6M06bBEAfQn2+ULmg06S7VxNsrRNTZjgf5kGcoR10wj6bfZZVCjNunQXAP0lmUoN73ZoNcsIVVVytE1IL0sQ7nKRztVrd9ZQoqqNFUFDuCQLEO5yRVblbME1oWB9LI8fVGnnfxBj5hV3ZLM0TUgrSwr5bpUbMNsWDZ5lgA0QZaVsjRYxRat5XjWME0edZolADhSskrZzDZJOk/Sbnc/pf3aqZI+KWmxpAOS/sbdv5NqDFLs43PDbFjWZZYAYGEpK+XNktbPe+2fJH3A3U+V9PftvyeVsuV4VL1ajzspel036kwCaKJkoezut0n6xfyXJR3b/vPvSHok1fVnVb0x1k2EpYjUzzcEMJiyN/oul/Q1M/uwWv+HcGbqC466MZaynTvCUkTTG3GAaMre6Hu7pHe7+ypJ71bradcLMrONZrbdzLZPT08PfcFRqtEyqsiqj5hFnkkATVR2pXyRpMvaf/4vSZ/p9InuPiVpSpImJyeHTsFRqtEmVJFVHrHjplLAkcqulB+R9Jr2n9dJ+kkZFx22Gm1CFVnVujZr2cDCUh6Ju0HS2ZKWmdlOSVdL+itJHzOzZ0naL2ljqusXoQmNGlWtazdhFgIMI1kou/uGDh/6g1TXLFpTbsAzO5OYWHqUdux9Qtt+vjf5ckKTbyoFdJNlm3VZIpyOKEvZTTZNmIUAw8i2zbofRTRFVH06oixlN9lEOKMNRFTbSjlye3VEZS8nNGkWAgyitpVyle3VObYtD9vyPYqmzEKAQdQ2lKs6zpbrUS+WE4AYart8UdVGUq5HvVIuJ9AkAvSvtqFc1XG2nI96pXiuH2v7wGBqG8pVbSRx1OuQmRnX3Q//KsuZA1CV2q4pS9VsJLE22zJbIW+5b3ftW9WBItW2Uq4KR71aZtfW//KsNcwcgAHUulKuymyFfvrEUknStp/vzeZoXFFm19Y//92dunTdSY2fOQD9olJOpOkbXLNr67t+vV+f3fqALn71Gi0ak845+Xj93ornNeJ/A2AYVMqJRH42YBnmrq3v+vV+Xfut+3XyCceWEsg5Nu8As6iUE8n5aFwRqlpbb/oMBflLVimb2SYz221mP5z3+rvM7D4zu8fMkj/NOqVuFVkVbcvRVHH6pekzFOQv5fLFZknr575gZq+VdL6k33f3l0v6cMLrJ9WrnZqjcdVowtNiUG8pb3J/m5lNzHv57ZI+6O6/bX/O7lTXT61XO3Xdjsbl0ipN8w5yN3SlbGZvHeLLXirpLDPbZmb/Y2avHPb6VeunIqvLXdByuskSMxTkbpRK+QOS/n2I671A0hmSXinpJjNb4+5H/NdtZhvVfobf6tWrRxhmGk2qyKLfZGl+Ff+6312ur9ZkhoLm6RrKZvaDTh+StHyI6+2UdHM7hL9jZjOSlkmanv+J7j4laUqSJicnw5VkTXl+nxT7JEm30xZVjw0YRq9KebmkP5T0y3mvm6RvD3G9L0p6raRvmNlLJT1b0p4hvk/l6rZm3E0Vs4J+17CjV/HAoHqF8lckHe3ud83/gJl9s9sXmtkNks6WtMzMdkq6WtImSZvax+SeknTRQksXuUhxq8uIyp4VDHLWOHIVDwyjayi7+8VdPvbGHl+7ocOH3tTHuBBI2bOCQarfJq3toxlosy5Zri3AZZ4kGeSsMactUDe0WZeIFuD+DFL9NmltH81ApVyilC3AuVbgCxm0+q3LeXBAolIuVapNqbpV4FS/aDIq5RKluklRnW7CM1vxb/v5XknS6RNLqX7RKIRyiVJtStXlJjw5tXMDqbB8UaJU0/J+NsZyuKEQjSAAoVy6FA0nvZo7cllzphEEIJRroVcFnksFOmojSA6zAaAXQrkmulXguVSgo7Rz5zIbAHohlBsgl1bkUdbcc5kNAL1w+iKwohpCcmpFHrYRpC4nUAAq5aCKnI43oRkjl9kA0AuVclBFN4SkbkWuus07p9kA0A2VclC5bM5JMTbZmjAbQDMkq5TNbJOZ7W7f0H7+x64wMzezZamu303VVV0/UrVkpxClzZsbE6EOUi5fbJa0fv6LZrZK0uskPZjw2h3l0sqb03ScTTagOMmWL9z9NjObWOBDH5F0paQvpbp2N7kcncppOp5yk42GEDRNqRt9Zna+pIfd/ft9fO5GM9tuZtunp4942PXQcqrqcpmOp6rqc5nVAEUqbaPPzJ4r6Sq1li56cvcpSVOSNDk5Wdh/hQtVdScuXaIl44t0x8/2UI0NIVVVn8usBihSmZXySyS9WNL3zWyHpJWS7jSzE0ocwxFV3YlLl+hd607ShVNbqcZGkKKqz2lWAxSltErZ3e+WdPzs39vBPOnue8oag3RkVbdkfJEunNpKNRYQDSFoopRH4m6QdIekl5nZTjO7ONW1BjW3qnvyqYO1qsZyOO7Xr5xOoABFSXn6YkOPj0+kuvYg6lSNRWjiKFJOJ1CAojS+zTpKNVZEhRuliaNIuZxAAYrS+DbrCNVYURVuTq3ZABbW+EpZqr4aK6rCzak1uwx1Wl9HcxDKARR19CvKUkwENJ4gV41fvoig02bjkvFFmpnxviv3CEsxUdB4glxRKQewUIV76bqTdOnnvjdwdVf1UkwUNJ4gV1TKAcxWuCs2nqEt9+3WwRnps1sf0K5f7w9f3UW9YVCdjjqiWQjlIMbGTE8+dVAf3/LTw16PfHoi8rnoUZ6MDVSJUA4kt+qu6nXbblU66+vIFaEcSG7VXZXnovup0mfX1yPOMoBOCOVAcqvuqqzsq67SgVQ4fdFBVY0HOZ2eqPJcNKcrUFdUyguIvIEVSZWVfW7r70C/qJQXEO3GPpHbhauq7OleRF0lq5TNbJOk8yTtdvdT2q99SNIfS3pK0s8kvdXdf5VqDMOKdGMfqvaF5bb+DvQrZaW8WdL6ea/dKukUd3+FpB9Lel/C6w8t0o19olXtRSiq8s9p/R3oV7JQdvfbJP1i3mtfd/cD7b9uVes5feFEmhrXbUOLGwUB3VW50fc2STcW/U2LaPuNNDWu24YWR9mA7irZ6DOz90s6IOn6Lp+z0cy2m9n26enpvr5vkVVYlKlxpKq9CHWr/IGilV4pm9lb1NoAPMfdO6alu09JmpKkycnJvlK1jlVYpKq9CMNW/lFvfAQUrdRQNrP1kq6U9Bp3f7Lo7x/p1ESR6tQuPEwrOSdQ0CQpj8TdIOlsScvMbKekq9U6bfEcSbeamSRtdfe/LuqadVt/raNhKv86zoCATpKFsrtvWODla1NdT8rvhj6DqssUftDKv64zIGAhtWqzrtv661xNnsIzA0KTZN9mPb8RQVKIUxNFq2MTSb/qdgIF6CbrSrlJ1WOTp/B1ngEB82VdKTepelyo9fvEpUu0ZHxRyBsVFS3KuXEgtaxDuUmNCPOn8CcuXaJ3rTtJF05tpV0ZqJGsly+atAE0fwq/ZHyRLpzamuSYWF1OeQA5yrpSbtoG0Nwp/JNPHUwyS+CGQUC1sq6Um7wBlGqWUESjBpU2MLysQ1mqVwvyIFI1yox6yqNJJ2KAFLIP5aZKNUsYtQKnJRoYDaGcsRSzhLkV+POf+2z9+eRKvfT4Y+TeqoJ7hX6Tz1MDRSCUcZjZCnztZWfpzgd/pau+cPdAyxBNOhEDpJD16YsoIj9tehhjY6YZ1zOBLPXfmNO0EzFA0aiUR1TXja1hlyGafCIGKEKyStnMNpnZbjP74ZzXXmBmt5rZT9r/fn6q65clt1bvfqv6UZ7oTUs0MLyUyxebJa2f99p7JW1x95MkbWn/PWs5tXoP0hjCMgRQjZQ3ub/NzCbmvXy+Wk8jkaTrJH1T0t+lGkMZctrYGuS4GssQQDXK3uhb7u672n9+VNLykq9fuJwqykGrepYhgPJVttHn7m5mHY8pmNlGSRslafXq1aWNa1A5VZQ5VfVAU5VdKT9mZi+UpPa/d3f6RHefcvdJd5887rjjShvgMHKpKHOq6oGmKrtS/rKkiyR9sP3vL5V8/UbLqaoHmipZKJvZDWpt6i0zs52SrlYrjG8ys4slPSDpglTXx8KaegMnIBcpT19s6PChc1JdEwByR5s1AARi7vHv02Bm02otd1RtmaQ9VQ8ioSa+vz3uPr/JSZJkZre0v6bX96irprzXKt5n59+7HEI5CjPb7u6TVY8jFd5fOd8jF015r9HeJ8sXABAIoQwAgRDKg5mqegCJ8f7K+R65aMp7DfU+WVMGgEColAEgEEK5D2a23sz+z8x+ambZ3wN6LjNbZWbfMLMfmdk9ZnZZ1WNKwcwWmdn3zOwrA3xNIx7U0OF9/oOZPWxmd7X/ObfKMRah0+96tJ8podyDmS2S9K+S/kjSWkkbzGxttaMq1AFJV7j7WklnSHpHzd7frMsk3Tvg12xWAx7UoIXfpyR9xN1Pbf/z1ZLHlEKn3/VQP1NCubfTJf3U3e9396ckfU6tm/XXgrvvcvc723/ep1Zwrah2VMUys5WSXi/pM4N8nbvfJukX814+X60HNKj97zeMPMCKdXiftdPldz3Uz5RQ7m2FpIfm/H2nahZas9pPijlN0rZqR1K4j0q6UtJMr0/sQ+0e1NDFO83sB+3ljeyXaeaa97se6mdKKEOSZGZHS/q8pMvd/fGqx1MUMztP0m53/27R39tbR5fqenzp3yS9RNKpknZJ+udqh1Ocbr/rEX6mhHJvD0taNefvK9uv1YaZjav1S3q9u99c9XgK9ipJf2JmO9RaelpnZv85wvfr+0ENOXP3x9z9oLvPSPq0Wst42evwux7qZ0oo9/a/kk4ysxeb2bMl/YVaN+uvBTMzSddKutfdr6l6PEVz9/e5+0p3n1DrZ/ff7v6mEb7l7IMapBo/qGE2pNr+VNIPO31uLrr8rof6mdI80of2caCPSlokaZO7/2PFQyqMmb1a0u2S7tahNderarLbfhgzO1vS37r7eX1+/jMPapD0mFoPaviipJskrVb7QQ3unvUmWYf3ebZaSxcuaYekS+asu2ap0++6WuvKYX6mhDIABMLyBQAEQigDQCCEMgAEQigDQCCEMgAEQijXSJ3vZoeYFrrDHEZDKNdEA+5mh5g2a+E7zGFIhHJ91PpudoipKXeYKxOhXB+NuZsdUGeEMgAEQijXR+3vZgc0AaFcH7W+mx3QFIRyTbj7AUnvlPQ1tR5zc5O731PtqFB37TvM3SHpZWa208wurnpMueMucQAQCJUyAARCKANAIIQyAARCKANAIIQyAARCKANAIIQyAARCKANAIP8P1NoXD9WiJi4AAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 360x360 with 5 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8ClgsrjFyuCt"
      },
      "source": [
        "The start of the show is the Sobol' sampler:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "M8NA1eA51fpk",
        "outputId": "ef25eb46-a83e-41a2-a53c-2c34556fbf0e"
      },
      "source": [
        "sobol_engine = qmc.Sobol(d=d, scramble=False, seed=rng)\n",
        "sample = sobol_engine.random(n=n)\n",
        "qmc.discrepancy(sample)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.00024001962608743987"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 13
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "id": "sxNasrNb134A",
        "outputId": "4c6bf5b7-78ba-467f-a778-8bda2ebad4c1"
      },
      "source": [
        "sns.pairplot(pd.DataFrame(sample), diag_kind=\"hist\", corner=True, diag_kws={\"bins\": 4})"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<seaborn.axisgrid.PairGrid at 0x7ff6b6f25990>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 14
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAAFlCAYAAADcR5KFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAXx0lEQVR4nO3da4xdV3nG8ef1ZIitXECyJwF8yUV1S02FIIySCEQbcVEnrhR/oCUxQm3SUCO1CVXTIqVqFWj6BVqJihQLGGhIQSqJy4fKbY3TSgURNRfZEmmEHaWauAGPA3hsohCInMTx2w/n2HPm5Nxnr73evff/J0WambM9s7R1eHj22mvtY+4uAEB+a3IPAADQQiADQBAEMgAEQSADQBAEMgAEQSADQBDn5R7ABFinhyJY7gEA3WjIABAEgQwAQRDIABAEgQwAQRDIABAEgQwAQVRx2VtfGzdv0bOLR7OOYWr6fL36ykuMIcAY3rxps44d/WHWMQDjsAo+frPvgM1MN37p4TLH8hoPfOxdjCHQGAa8v1mHjHCYsgCAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIAhkAAiCQAaAIMzdc49hLGa2X9KGPi9vkHSixOFExrlY1utcnHD3uRyDAfqpXCAPYmYH3X029zgi4Fws41ygKpiyAIAgCGQACKJugTyfewCBcC6WcS5QCbWaQwaAKqtbQwaAyiKQASCI83IPYFxzc3O+f//+3MNA9dmIxzGnhyKM9H5L1pDN7F4zO25m3+/zupnZPWa2YGZPmNlVo/zeEyfY6wCgnlJOWdwnadBOqOslbW3/t0vSFxKOBQDCSxbI7v5dST8dcMgOSV/zlkclvcHM3pRqPAAQXc455I2SjnZ8v9j+2Y/yDKc8Z864njn5C/3kZ6d06cVrdfn6C7RmzahTmgDqqhI39cxsl1rTGtqyZUvm0azOmTOu/Yd+rDv2PK5Tr5zR2uk1+uyH3q65t76RUAYaLueyt2OSNnd8v6n9s9dw93l3n3X32ZmZmVIGl8ozJ39xLowl6dQrZ3THnsf1zMlfZBvTmTOuI0s/1yNPn9CRpZ/rzBkWFgA55GzIeyXdZmb3S7pG0vPuXvvpip/87NS5MD7r1CtndPyFU7py5sLSx0Njn9zGzVv07OLR4QcmNDV9vl595SXGEGAMkvTmTZt17OgPJ/73yQLZzL4h6TpJG8xsUdInJU1Lkrt/UdI+SdslLUh6UdItqcYSyaUXr9Xa6TUrQnnt9BpdctHaLOPp19jf8vH3ZPs/iKrMrz+7eFQ3funhrGN44GPvYgxBxnB2HKuRLJDdfeeQ113SH6X6+1Fdvv4CffZDb39NI718/QVZxhOpsdPW0XSVuKlXJ2vWmObe+ka95ePv0fEXTumSi/K2wEiNPVpbB8rWiGdZRLtptWaN6cqZC3XtlRt05cyFWdvf2ca+drr1VsjZ2Ae1daAJat+QuQweLFJjj9TWgRxq35AjLjOLpl9jL/vKIlJbB3KofUOOdNOqSnJcWURq60AOtW/IZy+DO1XlMjjn3HeuK4tI8+tA2WofyFW9DD7bULff85B2fvkxbb/nIe0/9OPSQpkbbED5aj9lsdrL4FwbFXIvAeMGG1C+2geytHwZPG6Q5VyhkXvuO9oGFqAJGhHIk8rZUnM31Ag32Kq0jRooAoE8QM6WGqGhTnplUQTWj6OJCOQBcrbUCA01p9xz6EAOtV9lsRq5V2g0eQnYOKs8om2NByZFQx6g6S01p1GvTpjaQJ3QkIeIsq24aUa9OmFrPOqEhjwBWll6o16d5F4eCBSJhjyBiK2sjo19lDn0Km+NB7oRyBOItq049zbrnHLfeAWKxJTFBHJv2uhW9hKxSBs2uPGKOknakM1szsyeMrMFM7uzx+tbzOzbZvY9M3vCzLanHE9RorWyMht7xDbe5OWBqJeUnzo9JWm3pA9IWpR0wMz2uvvhjsP+UtIed/+CmW1T65OoL5/0b5bV3KK1sjIbOxs2gHRSNuSrJS24+xF3f1nS/ZJ2dB3jki5uf/16Sc9O+sfKbm6RWlmZjT3a/DlQJynnkDdKOtrx/aKka7qO+ZSk/zCz2yVdIOn9k/6xJje3Mhv7pG080rwzEFXuVRY7Jd3n7pskbZf0dTN7zZjMbJeZHTSzg0tLSz1/UdObW1mNfZI2HnHeGYgoZUM+Jmlzx/eb2j/rdKukOUly90fMbK2kDZKOdx7k7vOS5iVpdna25/+Ko618qKtJ2niTr16AcaRsyAckbTWzK8zsdZJukrS365gfSnqfJJnZr0paK6l3BR4i2sqHaIrcODJuG2/61QswqmQN2d1Pm9ltkh6UNCXpXnc/ZGZ3Szro7nsl/amkL5vZn6h1g+9md58oKaKtfIgk91Zvrl6A0STdGOLu+9Raytb5s7s6vj4s6d1F/b3cD1SPetMq95RBhIftA1XATr0C5G6gw+R+AA9XL8Bocq+yqIWIDxvq1O8BPG+8eG1pDySKtG4biIpALkD0m1a9bnh+/sPv0OEfvcBSNCAQpiwKEH2zRK8pA3fpt/7+IZaiAYHQkAtQhc0S3VMGx1+I3eqBJmp0Qy6qoVZxswRL0YB4GhvIRa+MGHfJXe6VDyxFA+JpbCA3vaGyFA2Ip7FzyLlXRkTY6p1zKVodPwMQWK3GNmQaaj7RN9IAuTS2IUduqHVvj9E30gC5NLYhR22oTWiPuW9oAlE1tiFLMbfzRmyPRTf2flu5WXKHpmt0IEeU+2ZjtxQbWCJMFwERNXbKIqrcNxu7pVgeuNrposiPOgVWg0AOJtqGjVTzvZM+u7oJc+xoLgK5LUrrinazsQmNHYiCQFa81pXzk0+6NaWxAxEQyKJ1DUJjB8qTdJWFmc2Z2VNmtmBmd/Y55kNmdtjMDpnZP6UcTz/jrGyo+6aNXiItD2SFBuosWUM2sylJuyV9QNKipANmtrf9waZnj9kq6c8lvdvdnzOzS1KNZ5BRW1e0qY0mitbYgSKlbMhXS1pw9yPu/rKk+yXt6DrmDyTtdvfnJMndjyccT1+jtq6ImzaaqKlbzlF/KeeQN0o62vH9oqRruo75ZUkys/+WNCXpU+6+P+GYehq1dXFDKS6uXlAHuXfqnSdpq6TrJO2U9GUze0P3QWa2y8wOmtnBpaWlJAMZZZ6ULb/D5WqpXL2gDlIG8jFJmzu+39T+WadFSXvd/RV3/z9J/6tWQK/g7vPuPuvuszMzM8kGPAw3lAYr+3MCO0Xbcg5MIuWUxQFJW83sCrWC+CZJH+465l/UasZfNbMNak1hHEk4plWJcEMpygaWXnIuH2Q5HOogWUN299OSbpP0oKQnJe1x90NmdreZ3dA+7EFJJ83ssKRvS/qEu59MNaYi5P6UjVwNdBQ5WypXL6iDpBtD3H2fpH1dP7ur42uXdEf7PwwRfQNLr5Z62fp1Wjc9pUeePpG00Ue4egFWi516FRJ9lUf3NuvL1q/T7e/dqhvnHy1l5UOkLefAJAjkCok+T9rdUtdNT50LYyleoweiyb3sLYm6bhCowjxp5xz7iy+/ysoHYAy1a8h13iBQtXnS6I0eiKZ2DTniBoEiG3ukB/0MU4VGD0RSu4Yc7cZXnRv7MFVr9EButWvI0bY3l93Yo82fV6nRA7nVLpCjXSaXuVki+sYRAIPVbsoi2mVymTe2om8c6SXyVnCgbBM3ZDO7pciBFCnSZXKZjb1qD9ih0QMrraYh/5WkrxY1kDKV2crKbOxVW2ZWxUYPpDQwkM3siX4vSbq0+OGkl2PVQ1lbeqN9QvQw0VbEALkNa8iXSvpNSc91/dwkPZxkRIlFbGVFNfZo8+fDVK3RA6kNC+R/k3Shuz/e/YKZfSfJiBKL1sqKbuxVesBO1Ro9kNrAQHb3Wwe81v2w+UqI1soiNvayVK3RA6nVbh3yME1epyyxcQSIrHbrkIeJ1srKbOxN3sYNVEHjGrI0uJWV3SDLbOwRH7wEYFnjGvIguZbEldXYo93QBLBS0oZsZnNm9pSZLZjZnQOO+6CZuZnNphzPMLkaZFnzqNEevDSuaPPfQNGSBbKZTUnaLel6Sdsk7TSzbT2Ou0jSH0t6LNVYRlW1rcfjinZDcxxss0YTpJyyuFrSgrsfkSQzu1/SDkmHu477a0mfkfSJhGMZSbQlcUVb7fRIzgcBNXl5IJoj5ZTFRklHO75fbP/sHDO7StJmd//3hOMYWZUb5KgmnR7J3VDrfvUCSBlv6pnZGkmflXTzCMfukrRLkrZs2ZJsTLmXxEV+FGXuhlr3qxdAStuQj0na3PH9pvbPzrpI0q9J+o6ZPSPpWkl7e93Yc/d5d59199mZmZmEQ863USF3Ax0md0NtwtULkLIhH5C01cyuUCuIb5J0bru1uz8vacPZ79vPxvgzdz+YcExh5W6gw+RuqLmvXoAyJGvI7n5a0m2SHpT0pKQ97n7IzO42sxtS/d2qGreB1nkDSz9ss0bdJZ1Ddvd9kvZ1/eyuPsdel3Is0Y3TQOu+gQVoqkZunY5onAYabQMLGzaAYrB1OohxGmikLdA8sAgoDg05kFHnSCNtgY74wCIaO6qKQK6gCDfYzsq9HK5b9OWDwCBMWYwo0qaNSDfYci+H6xZ9+SAwCIE8gojzpFE+Oy/a5+JFml8HxkUgj6CKrausRh+prUvxGjswDgJ5BFVrXWU3+ihtXYrX2IFxEMgjqFrrqmKjL0q0xg6Mg1UWI5h0VUOu5VfRVj6UjS3WqCoa8ggmaV05bwRWrdEDaKl1Qy6yoY7bunJumIi0ThnA6GrbkHMvVct5I5B5VKCaatuQc2/pzb29Ofc8KtuXgfHVNpBz39hq8rQB25eBydR2yqLXja3L1q/TuukpPfL0ieTbn5s8bdDkZXfAatQ2kLs3CFy2fp1uf+9W3Tj/aKkPdY+yYaJMVdtIA0RR20DubqjrpqfOhbFEa0tpkmV3kR7eBORS2zlkaeWNrRdffrXRmyXKNO78OXPOQEvShmxmc5I+J2lK0lfc/dNdr98h6aOSTktakvT77v6DFGNhs8RwRbXUcefPmXMGWpI1ZDObkrRb0vWStknaaWbbug77nqRZd3+bpG9K+ptU42nyqodRFN1Sx1l2l3tFDBBFyoZ8taQFdz8iSWZ2v6Qdkg6fPcDdv91x/KOSPpJqME1e9TCKnC2VqxegJeUc8kZJRzu+X2z/rJ9bJX0r4XjYLDFAzpbK1QvQEmKVhZl9RNKspN/o8/ouSbskacuWLSWOrDi5t3IPk7OlcvUCtKRsyMckbe74flP7ZyuY2fsl/YWkG9z9pV6/yN3n3X3W3WdnZmaSDDa13Fu5h8ndUgddvUS+sgCKlLIhH5C01cyuUCuIb5L04c4DzOwdkr4kac7djyccS3bRN0tEbanRryyAIiVryO5+WtJtkh6U9KSkPe5+yMzuNrMb2of9raQLJf2zmT1uZntTjWdcRbey3A8bGkXuOfZeol9ZAEVKOofs7vsk7ev62V0dX78/5d+fVIpWxme9TSb6lQVQpBA39aJJsQRstVMCTd1azJI4NAmB3EOqVjbpw4aaPI/KlQWahEDuIVorK3vTRqQ2HvVmI5ACgdxDtFZW5jxqxDbe1MeYonkI5B6itbIyGzsP+gHyqfXjN1cj0hKwMjdtjLOFmg0bQLFoyBVQZmMftY1HnNoAqo6GXBH9GnvRLXXUNs6GDaB4NOQKS9FSR23jbNgAikdDXqWc86ipWuoo8+dV2AoOVA2BvAq5PwuOZxgD9cKUxSrkXiLW9GcYR9rAAhSBQF6F3POouTew5NywwSoP1BGBvAq9Gupl69dp3fSUHnn6RPLWFqGl5pL76gRIgUBehe6Getn6dbr9vVt14/yjpbW2pm4rzn11AqRAIK9Cd0NdNz11LowlWltK0R4ABRSBVRar1LlE7MWXX8226qGfum5vZpUH6oiGXKBora3ON76aPH+O+qIhFyhaa4u4vbnIxh7pAVBAEZI2ZDObk/Q5SVOSvuLun+56/XxJX5P0TkknJd3o7s+kHFNK0VpbtBtfdW7sQBGSNWQzm5K0W9L1krZJ2mlm27oOu1XSc+7+S5L+TtJnJvlbkeZJI7W2aNubIzZ2IJKUUxZXS1pw9yPu/rKk+yXt6Dpmh6R/bH/9TUnvM7OxEiz39uXIok2h5NzqDVRByimLjZKOdny/KOmafse4+2kze17SekknRv0jbBDoL9oUSrSbnkA0lbipZ2a7zOygmR1cWlpa8Rqta7BIUyjRGjsQTcqGfEzS5o7vN7V/1uuYRTM7T9Lr1bq5t4K7z0ual6TZ2dkVcxG0ruqI1tiBaFI25AOStprZFWb2Okk3SdrbdcxeSb/X/vq3Jf2Xu481+UvrqpZIjR2IJllDbs8J3ybpQbWWvd3r7ofM7G5JB919r6R/kPR1M1uQ9FO1QnsstC4AdZF0HbK775O0r+tnd3V8fUrS76z27zT1ATsA6qUSN/UAoAkIZAAIwsa8h5admS1J+kGflzdojDXMNce5WNbrXJxw97lh/9DM9rf//ai/t6k4Fyt1n4/R3m9VC+RBzOygu8/mHkcEnItlqc4F53gZ52KlSc8HUxYAEASBDABB1C2Q53MPIBDOxbJU54JzvIxzsdJE56NWc8gAUGV1a8gAUFmVC2QzmzOzp8xswczu7PH6+Wb2QPv1x8zs8vJHWZ4RzsfNZrZkZo+3//tojnGWwczuNbPjZvb9Pq+bmd3TPldPmNlVI/5e3nNtvN+WJXm/uXtl/lPrmRhPS7pS0usk/Y+kbV3H/KGkL7a/vknSA7nHnfl83Czp87nHWtL5+HVJV0n6fp/Xt0v6liSTdK2kxwo6x414z/F+S/9+q1pDLuVTSCpklPPRGO7+XbUeUtXPDklf85ZHJb3BzN405NfynlvG+61Divdb1QK516eQbOx3jLuflnT2U0jqaJTzIUkfbF8yfdPMNvd4vSlGPV/j/pumvOd4v41n7Pdb1QIZ4/tXSZe7+9sk/aeWmxyQAu+3VahaII/zKSQa9CkkNTH0fLj7SXd/qf3tVyS9s6SxRTTK+2eSf9OU9xzvt/GM/X6rWiCX8ikkFTL0fHTNWd0g6ckSxxfNXkm/2777fa2k5939R0P+De+5ZbzfxjP2+y3pA+qL5iV9CklVjHg+Pm5mN0g6rdb5uDnbgBMzs29Iuk7SBjNblPRJSdOS5O5fVOvDErZLWpD0oqRbhv1O3nPLeL+tlOL9xk49AAiialMWAFBbBDIABEEgA0AQBDIABEEgA0AQBHKNDHsSF1CkYU87w/gI5JowsylJuyVdL2mbpJ1mti3vqFBz90ka+knKGB2BXB88iQulGuFpZxgTgVwfkzzJDEAgBDIABEEg18ckTzIDEAiBXB+jPJUMQGAEck20P6ni7JO4npS0x90P5R0V6qz9tLNHJP2KmS2a2a25x1R1PO0NAIKgIQNAEAQyAARBIANAEAQyAARBIANAEAQyAARBIANAEAQyAATx/+9+WffLCwnsAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 360x360 with 5 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "19rER8C12Ito",
        "outputId": "5841785f-4400-4f03-daa1-050abed83ead"
      },
      "source": [
        "sobol_engine = qmc.Sobol(d=d, scramble=True, seed=rng)\n",
        "sample = sobol_engine.random(n=n)\n",
        "qmc.discrepancy(sample)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.0001551092547196209"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hgTaZwB128hH"
      },
      "source": [
        "We spent quite some time on Sobol'. First to ensure it starts with 0, and then to have some warnings if users do not ask for $2^n$ points. Who would dare doing this!?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ff5Prbx22seV",
        "outputId": "0ece6166-3217-46c7-c59c-9b5325ca95bf"
      },
      "source": [
        "sobol_engine = qmc.Sobol(d=3, scramble=False, seed=rng)\n",
        "sample = sobol_engine.random(10)\n",
        "sample"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/scipy/stats/_qmc.py:1078: UserWarning: The balance properties of Sobol' points require n to be a power of 2.\n",
            "  warnings.warn(\"The balance properties of Sobol' points require\"\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[0.    , 0.    , 0.    ],\n",
              "       [0.5   , 0.5   , 0.5   ],\n",
              "       [0.75  , 0.25  , 0.25  ],\n",
              "       [0.25  , 0.75  , 0.75  ],\n",
              "       [0.375 , 0.375 , 0.625 ],\n",
              "       [0.875 , 0.875 , 0.125 ],\n",
              "       [0.625 , 0.125 , 0.875 ],\n",
              "       [0.125 , 0.625 , 0.375 ],\n",
              "       [0.1875, 0.3125, 0.9375],\n",
              "       [0.6875, 0.8125, 0.4375]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 16
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fGR-VQjO1j1k"
      },
      "source": [
        "## Sample from another distribution\n",
        "\n",
        "The following is also new to SciPy. We can sample from any arbitrary distribution using a MC or QMC engine."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nVRBvo431PAq"
      },
      "source": [
        "Just defining a helper function to do a nice plot of the PDF"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Y5CZ2NgUtefs"
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "\n",
        "def plot_pdf(dist, sample):\n",
        "    fig, ax = plt.subplots()\n",
        "\n",
        "    x = np.linspace(dist.ppf(0.01), dist.ppf(0.99), 100)\n",
        "    pdf = dist.pdf(x)\n",
        "    ax.plot(x, pdf, '-', lw=5, label='fisk pdf')\n",
        "\n",
        "    delta = np.max(pdf) * 5e-2\n",
        "    ax.plot(sample, -delta - delta * np.random.random(99), \"+k\")\n",
        "\n",
        "    ax.hist(sample, density=True, histtype='stepfilled', alpha=0.2)\n",
        "    ax.legend(loc='best', frameon=False)\n",
        "    plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6dXW9GrB14qJ"
      },
      "source": [
        "First we sample from an arbitrary distribution-here Fisk-using MC."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "uKbfC9RHteiK",
        "outputId": "3d2a5d83-cf47-4368-979d-b99656480540"
      },
      "source": [
        "from scipy.stats import fisk\n",
        "\n",
        "c = 3.9\n",
        "dist = fisk(c)\n",
        "\n",
        "sample = dist.rvs(99, random_state=rng)\n",
        "plot_pdf(dist, sample)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXwb5ZkH8N8jyfJ9JXYOX3FiEhIn5DQ5SEkIRwnQJhxJy1WWUqBsN2zvLmxbCJRuge6yWwqUBUopLDQk4TLkKjQJIZDLOUjiXDjOYTtOfMW3ZVnSs39YBlkzumxJI42e7+fjT+R3RppHivXT6J133iFmhhBCiOhn0LoAIYQQwSGBLoQQOiGBLoQQOiGBLoQQOiGBLoQQOmHSasNZWVlcWFio1eaFECIq7d69u4GZs9WWaRbohYWFKCsr02rzQggRlYjolKdl0uUihBA6IYEuhBA6IYEuhBA6IYEuhBA6IYEuhBA6IYGuoXOtFhw524o2S4/WpQghdECzYYuxqvp8J1bvrsb6g2dx5GwbAMBAwPgRabikaCj+6ZJC5A9J0rhKIUQ0kkAPo5VlVfjVuwdhtTn6tTsYOFTbikO1rXh9x2k8dv0k3DQjT6MqhRDRSrpcwoCZ8cd/fIFfrN6vCHN3XT12/HTV5/jJyn3otNrCVKEQsenpp5/GhAkTcNttt6G0tBSPP/64x3VfeeUVLFu2LCjbvfPOO7F69WoAwCeffIKJEydi6tSp6OrqGtTjyh56iDEzHi4tx6vbPJ7cpertPTU409yFV++aBbNJPneFCIXnnnsOH330EfLyer8RL1q0KOw1vP7663jwwQdx++23D/qxJNBD7LXtp7yGeVZKPBrau1WXba9swq/ePYAnbpoMIgpViUJoqvCBNSHfxsnHr1O03XfffaisrMQ111yDu+66C5mZmSgrK8MzzzyDVatW4ZFHHoHRaER6ejq2bNnS775r1qzBY489hvfffx9ZWVlfti9fvhzHjx9HRUUFGhoa8Itf/AL33HMPmBn3338/PvzwQ+Tn58NsNgMAXnrpJaxcuRIbNmzAunXr8Prrrw/qeUqgh9Cpxg78bu0R1WXfnJKDf1t4IfIyk3Cu1YLH1x3BO3trFOutLKvGuOGpuPvSMaEuV4iY8vzzz2P9+vXYtGkTsrKy8Morr3y57NFHH8WGDRuQm5uL5ubmfvd755138NRTT2Ht2rXIzMxUPO7+/fuxfft2dHR0YNq0abjuuuuwfft2HD16FIcOHcK5c+dQXFyMu+66C3fffTe2bt2Kb3zjG1iyZMmgn5PP7/JE9DIR1RHRQQ/LiYieJqIKItpPRNMHXZUOOByMn6/aj64eu2LZPZeOxh++PRV5mb2jWYanJeCpb03B75dMhsmg3BP/7drD+PhYfchrFkL0mjt3Lu688068+OKLsNu/eg9v3LgRTzzxBNasWaMa5gCwePFiJCYmIisrCwsWLMDOnTuxZcsW3HLLLTAajcjJycHll18ekrr96Zx9BcBCL8uvATDW+XMvgD8Nvqzo98pnJ7HzZJOi/dZZBfjldcUwuAU3EWFpST4eWTxRcR9m4IG39qO9Ww6SChEOzz//PB577DFUVVVhxowZaGxsBAAUFRWhra0Nx44d83hf9+7RcHaX+gx0Zt4CQJlMX1kM4FXutR1ABhGNDFaB0aiu1YInNyi7WvIyE/HLayd4ve9ts0bhu3MLFe21LRb819+PBqtEIYQXx48fx6xZs/Doo48iOzsbVVVVAIBRo0bhrbfewh133IHy8nLV+7733nuwWCxobGzE5s2bcfHFF2PevHl48803YbfbUVtbi02bNoWk7mD0oecCqHL5vdrZVuu+IhHdi969eBQUFARh05HpxU8qYelRDk/8/ZIpSI73/ZL/8toJOHq2DZ8db+zX/tfPTuKGabmYnJcRtFqF0JraAUut/fznP8cXX3wBZsYVV1yBKVOmYN++fQCA8ePH4/XXX8fSpUvx/vvvo6ioqN99J0+ejAULFqChoQG//vWvkZOTgxtuuAEbN25EcXExCgoKMGfOnJDUTczseyWiQgAfMPMklWUfAHicmbc6f/8HgH9jZq9XrygpKWE9XuCiqcOKuY9vVPSdf2f2KPzmesXL51FVUyeu+u+PFR8MxSPTULpsLkxGGcooRKRZvnw5UlJS8LOf/Sxk2yCi3cxcorYsGKlQAyDf5fc8Z1tM+sunJxRhnhBnwA+vHBvQ4+QPScKPrhynaD9U24oVu6pU7iGEiHXBCPRSAHc4R7vMBtDCzIrulljQ0tWDVz49qWi/deYoZKXEB/x43/vaaIwfkapo/+PGL2BRGT0jhNDW8uXLQ7p37os/wxb/BmAbgAuJqJqIvkdE9xHRfc5V1gKoBFAB4EUAPwhZtRHutW0n0eY2EsVsNODeeQMbQx5nNOB3N16kaD/X2o3/2x7YmadCCP3zeYSOmW/xsZwB/EvQKopSNrsDf1U5I3RJSR5GpCcM+HGnFWTimkkjsO7g2X7tz20+jptnFiDFj4OsQojYIEfWgmTz0XrUt/U/hd9oIPzz/CIP9/DfT64aB/ehrE0dVry89cSgH1sIoR8S6EGyarfyQOVVE4YHZW7zscNTccPUXEX7i1sq0dIlF8cQQvSSQA+ChvZu/ONwnaL9WxcHb07zH105TjEtQFu3DX/beTpo2xBCRDcJ9CB4d28NbI7+4/mHpcZj3tjsoG2jYGgSlpbkK9r/8ukJdNtkxIsQQgJ90JgZK8uU3S03Ts8L+sk/3583RtGXfq61G+/tOxPU7QghopME+iDtr27BsXPtivalJcG/hFxhVjIWThyhaH9xSyUcDt9n/Aoh9E0CfZDU5jAvGZWJouyUkGxPbUz7F3Xt2HxM2YcvhIgtEuiD4HAw1ruNDwdCs3feZ1pBJmaOHqJof3GLDGEUItZJoA/C/poWnG219GszGQgLJ4V29uDvq+ylb6tsxLFzbSHdrhAiskmgD4La3vklF2QhPTEupNtdcOEwFGUnK9pfC/BC1EIIfZFAHyBmxvqDyjnI1A5aBpvBQPjO7FGK9rf3VKPNIicaCRGrJNAH6Ni5dpxs7OzXRgRcVTw8LNu/cUYekszGfm0dVrvqQVohRGyQQB8gte6Wi0cNQXZq4NPkDkRaQhxunK6cDuDVbafgz0VLhBD6I4E+QOvLlYF+9aTQd7e4umNOoaKtoq4d2yoblSsLIXRPAn0ATjd24nBtq6L96onh6W7pM254KmapDGF8Y4fM7yJELJJAH4BNR5Un8VyUm468zMHPrBgotb30v5efQ1OHNey1CCG0JVdHGICPj9Ur2q6c4LZ3fmZvWGq5KpORlWhAQ9dXF5O22h14e0817r50YFdKEkJEJ9lDD5Clx45tx5V91JddGLyZFQNhNhJuKlZ+M3hj52k5OCpEjJFAD1DZyfPocrtA85BkMy7KTdeoIuDmicpAr6zvwK6T5zWoRgihFQn0AH2sMgnWpWOzYHC7+EQ4jc4wYXauWdG+Qi5+IURMkUAPkFr/+fxx2nS3uLplknIvfc2BWrR0ypmjQsQKCfQAnGnuUp37/NIgXplooK4uSkRGUv85ZLptDpTul4tfCBErJNADsEVl73xiTlrYzg71JsFEuGGa8szR1SpXUxJC6JMEegAitbulz9IZymuOfl7dgiNnlSdBCSH0RwLdT3YHY2tFg6I9kgK9OCdNdbTNqrJqDaoRQoSbBLqfys+0oM1i69eWEm/C9FGZGlWkTu1qSe/urYHV5lBZWwihJ34FOhEtJKKjRFRBRA+oLC8gok1EtJeI9hPRtcEvVVufqZxMNHP0EMQZI+szcdGUHJhN/Wtq7LBi4xG55qgQeufz1H8iMgJ4FsBVAKoB7CKiUmY+5LLarwCsZOY/EVExgLUACkNQ78AE4TT8bYeUgX5JliVsp/j7KyPJjK8XD8cH+/tffGNVWRUWhnk2SCFEePmzezkTQAUzVzKzFcAKAIvd1mEAac7b6QB0NVaux87YdUY52dXsPO1Ht6j5Vony4OjmY/Woa7OorC2E0At/Aj0XgOvYt2pnm6vlAG4nomr07p3fr/ZARHQvEZURUVl9vXLESKTaX9eDzp7+86KkxxOKsyNzbrO5F2QhJz2hX5vdwXhvr64+Z4UQboLVAXwLgFeYOQ/AtQBeIyLFYzPzC8xcwswl2dmRMzrEl23V3Yq22XnxMJB2p/t7YzQQbpyuPDi6ene1TNglhI75E+g1AFy/w+c521x9D8BKAGDmbQASAGQFo8BI8FmVMtAvidDulj43zVAG+tFzbThYI2PShdArfwJ9F4CxRDSaiMwAbgZQ6rbOaQBXAAARTUBvoEdPn4oXFhtjd62y/3xOnnIyrEgyOisZJSpDKt/aI2PShdArn53AzGwjomUANgAwAniZmcuJ6FEAZcxcCuCnAF4koh+j9wDpnayT7/Z7z1rR3X+2XGQlGjB2SAT2n7uNuFlS5EDZqf6rvLvnFB6c2o14UxC7i3KmBe+xhBAD5lcqMfNa9B7sdG17yOX2IQBzg1taZFDvPzeDIrT/3NW1YxOxfEsrLLavPlubLYyNJy245oJEDSsTQoRCZJ0VE4F21Kh1t0R2/3mftHgDFhYlKNrfOtypQTVCiFCTQPei28bYezZ6xp+rWTJBOU/6ppPdqO+0q6wthIhmEuheHKizwqrSfz4mw6hNQQMwJ8+MkSn9/5vtDLx3tEujioQQoSKB7sVOlbNDL86Jjv7zPkYD4cbxyr301Yc6ZUy6EDojge7FTpX+84tVrt0Z6W5S6XY50mhDeb1NZW0hRLSSQPfA7lAffz4zJ3r6z/uMyTRh+og4RftqOTgqhK5IoHtwuMGGNmv/LolUM2FCVgSOP/fDkmLlXnrp0S5Y7dLtIoReSKB7sOuMcvz59JFmGA3R03/u6rqxiYh3O5bbZHFg40mZgVEIvZBA90BtutyZOdHXf94nPd6ArxcpTyZafUi6XYTQCwl0FcysekB0ZhQeEHW1dIIy0Ded7EZdh4xJF0IPJNBVnGi2o6Gr/zU4zUZg8rDoDvS5+fEyJl0IHZNAV6HWfz51uDm4E1ppwNOY9FUyJl0IXZBAV1GmMlzx4ijuP3elNtrlWJMN++t6NKhGCBFMEugqdqscEC3RSaCPzjChZKTyuaySg6NCRD0JdDeNnXZUNisPEk4foY9AB4ClxcqDo6XHuvpNsyuEiD4S6G7Uzg69cKgJ6Qn6eamuG5uIRLfjAa3djA3H5eCoENFMPykVJGqBPkOliyKapZgNuHascp70N8ul20WIaCaB7kbtgKhan3O0+5bKwdHPqq2oapEJu4SIVhLoLiw2xoFzytEeejkg6mpmjhmjVeZ1XyUTdgkRtSTQXRyos8La/3wiZCcZkJ8WPRe08BcRYanKXvqqQ52wO+TgqBDRSALdRZnacMWR0XVBi0AsmZAEo9tTq2134JPTyhOrhBCRTwLdheoBUR12t/QZlmzEgkLl/O5ycFSI6CSB7sTM2F2r0n+uwwOirtQOjn54wiIXkRYiCkmgO1U223He0r8DPcFEmJitvNKPniwoTEB2Uv8/A5tDptUVIhpJoDupne4/ZXgc4tw7mXUmzkiqe+kryjvhkAm7hIgqfgU6ES0koqNEVEFED3hY51tEdIiIyonojeCWGXpq/ed6727pc/NEZaCfarFjW7XyNRFCRC6fgU5ERgDPArgGQDGAW4io2G2dsQAeBDCXmScC+FEIag2pWDhD1JP8dBMuLVAeHH3jYIcG1QghBsqfPfSZACqYuZKZrQBWAFjsts49AJ5l5vMAwMx1wS0ztJotDlScV54hOU1HE3L5cusk5V76349b0CAHR4WIGv4Eei6AKpffq51trsYBGEdEnxLRdiJaGKwCw2GPyt55UaYJmYmxc4jhytEJyHI7ONrjkGl1hYgmwUosE4CxAC4DcAuAF4kow30lIrqXiMqIqKy+vj5Imx68WO4/7xNnJCydoNxLf/2gnDkqRLTwJ9BrAOS7/J7nbHNVDaCUmXuY+QSAY+gN+H6Y+QVmLmHmkuzs7IHWHHRqE3LFSv+5q1snJcF9TE91qx1b5MxRIaKCP4G+C8BYIhpNRGYANwModVvnXfTunYOIstDbBVMZxDpDpsfO+FxlQq7pI/U9/lxNfroJl6mcOfrafjk4KkQ08BnozGwDsAzABgCHAaxk5nIiepSIFjlX2wCgkYgOAdgE4OfM3BiqooPpcEOP4ko9GQmEokyTRhVp6zsXJSvaNp3slml1hYgCfqUWM68FsNat7SGX2wzgJ86fqKI6XHGEfifk8mX+qHjkphpR0/bV6BYG8MbBTvzb3DTtChNC+BQ7wzg8UOs/nx6D/ed9jAbCbRcpD46+Wd4p1xwVIsLFdKAzs+op/7F4QNTVt4qTYHb7y2iyOPD+MbnmqBCRLKYDvabNjrMd/SfkMhl653CJZVlJRlw3NlHR/pd97WCZ30WIiBXTga7Wfz4xOw6JcTH9sgAA7pyqPDh6qMGGXSrfaIQQkSGmk0tOKPJsynAzpo9QflN55XMZwihEpIrpQFe95JyOr1AUqDunpijaNhy3oKZNhjAKEYliNtDbrQ4caVQGU6wfEHV1TVEChif3/xOxM/Cq7KULEZFiNtD3nrXCfYqS/DQjhiUbtSkoAsUZCbernGj0xsFOtFsdKvcQQmgpZgNdtbtF9s4Vbp2UBLPbZ1ybleVC0kJEoJgNdNUzRKX/XGFokhE3qczC+PK+DthkFkYhIkpMBrrNwdh7Vjkhl+yhq7t7mvLgaE2bHWsrLBpUI4TwJCYD/WijDR09/fcuU82EcUNjc0IuX4oyTbhytHIWxhd2y4lGQkSSmAx0tf7z6SPNMMTohFz+uGe6ci/9YH2PXEhaiAgSk4G+64zygg0yXNG7mTlm1SkRnitr06AaIYSamAt0ZlY9ff1iOSDqFRHhXpW99K1VVuyratagIiGEu5gL9OpWO865TcgVZwCmxviEXP5YWJSAMRnKcfrPbqrQoBohhLuYC3S1vfNJw2RCLn8YDYT7SlIV7R8eOoejZ6XrRQitxVyKqQX6TOlu8dv1FyYiJ0W5l/7cZtlLF0JrEugASnKUQ/KEOrOR8P0Zyr709z8/g+P17RpUJIToE1OB3tRlR8V55YRcckJRYL49MQlZif3/dBwMPP2PLzSqSAgBxFigq40/HzvEhMzEmHoZBi3BRKrj0ks/P4OKOulLF0IrMZVkMlwxeL4zWbmXzgz8z0eyly6EViTQJdAHJCnOgPtU+tLXHKiVES9CaCRmAr2zx4GD9coJuSTQB+62i5KQnaTcS3/qw6MaVSREbIuZQN9T2wOb2zUZRqYYkJsqF7QYqMQ4A/65ROUydeXnsOf0eQ0qEiK2xUygb69Rzt8yKzceJBNyDcqtk5IxIi1B0f74uiMyE6MQYeZXoBPRQiI6SkQVRPSAl/VuIiImopLglRgcO2qU/eezc6W7ZbASTIQfXzVW0b7zRBM2H63XoCIhYpfPQCciI4BnAVwDoBjALURUrLJeKoAfAtgR7CIHq6vHgX1nVQI9T04oCoabpuehKFt57dEn1h+BXa5qJETY+LOHPhNABTNXMrMVwAoAi1XW+w2AJwBE3GVs9p7tQY9b//nwZANGpUv/eTCYjAb8YuF4RfuRs214a3e1BhUJEZv8CfRcAFUuv1c7275ERNMB5DPzGm8PRET3ElEZEZXV14fv67ha//nsPOk/D6avFw/H9IIMRfuTG46izaIcXSSECL5BHxQlIgOApwD81Ne6zPwCM5cwc0l2dvZgN+237SpX1Zkl/edBRUR48NoJivaG9m48u+m4BhUJEXv8CfQaAPkuv+c52/qkApgEYDMRnQQwG0BppBwYtfTY1fvPc6X/PNguLhyC6y4aqWh/eesJnGrs0KAiIWKLP4G+C8BYIhpNRGYANwMo7VvIzC3MnMXMhcxcCGA7gEXMXBaSigO05/R5WN36z4clGzBa5UINYvAeuGY8zKb+f1ZWuwO/XXNYo4qEiB0+A52ZbQCWAdgA4DCAlcxcTkSPEtGiUBc4WDsqmxRtMv48dPKHJOGeS0cr2v9+6Bw2Ha3ToCIhYodffejMvJaZxzFzETP/1tn2EDOXqqx7WaTsnQPAtspGRZuMPw+tH1x2AYalKru0Hn6vHJYeuwYVCREbdH2maKfVhr0qp6DLAdHQSo434cFrlcMYTzd14jm5/qgQIaPrQN95ogk99v4ntoxINqAo06RRRbHj+qm5mD1miKL9+Y8rUSlXNhIiJHQd6J9WNCja5hZI/3k4EBEeu34S4oz9X2ur3YEH3j4Ah5xBKkTQ6TrQt1Yo+8+/li/DFcPlgmGpuOfSMYr2nSea8PrO0xpUJIS+6bbvoaG9G4drWxXtcyXQg+/MXo+L7p/gQOkeI6pb+x8MfXxNORZk1CEvbYB/gjnTBnY/IXRMt3vonx1X7p2PG2LCsGQZfx5OiXEGPH65ckqAjh7GgxtbZIpdIYJIt4H+6Rfq/eci/L5WEI+bJyYp2j853Y03DnZqUJEQ+qTLQGdmbFU5ICr959r596+lYUSy8s/tsU9aUXnepkFFQuiPLgP9VGMnapq7+rUZScafaykt3oD/UOl66bIxfvz384rhpUKIwOky0D/5Qjk177QRZqSYdfl0o8bloxNUu14+P9eDp3e2aVCREPqiy4RTu/TZ3HzZO48Ev740TfXCIs/sasenVcp564UQ/tNdoFt67Pj0uLL//LJC5YWMRfglmw34769nwu18IzCAH64/j7oOmetFiIHSXaDvONEEi9v15oYkGDB5WJxGFQl300eacf/MVEV7Q5cD968/D5ucRSrEgOgu0DcdUU7ROn9UPIwGOd0/ktx/cQrm5Cm7wXbUWPGf26Q/XYiB0FWgMzM2qgT6ZYUyXDHSGA2EP1ydiawk5Z/g87vbUXqsS+VeQghvdBXoJxo6cLqp/4kqBgLmj5L+80g0LNmIp6/OhNqXp1981IyDdcpLBwohPNNVoG9SGd0yvSATGQm6epq6ckl+PH46W9mfbrEx7v1ADpIKEQhdJd1mlUucLRg/TINKRCB+UJKC6y5Qfos6027H995vQqfbQW4hhDrdBHp7t031+qGXXZitQTUiEESE31+VgfFZypkXD9T14P51MvJFCH/oJtA3H62D1d5/T25YajyKR6ZpVJEIRFKcAS9+YwiGqHSP/eNkNx7aLDMzCuGLbgJ9/cGzirYrJgyXqxNFkfw0E1765hDEq8xw/MbBTjz5mQxnFMIbXQS6pceuOv78mkkjNKhGDMb0kWb8YWEm1D6G/7S7Hc+VSagL4YkuAn3rFw3osPYfDZGWYMLsMUM1qkgMxsKiRDw0T72r7MnP2vCXfXKRaSHU6CLQ15cru1uunDAcZpMunl5M+u7UFCy7OEV12SNbWvHSJ5VhrkiIyBf1iddjd+Cjw+cU7QuluyXq/XR2Ku6ckqy67LE1h/G/Hx8Pc0VCRDa/Ap2IFhLRUSKqIKIHVJb/hIgOEdF+IvoHEY0KfqnqdlQ2obmzp19bYpwR88bJcMVoR0R4aF4abpqQqLr8d+uO4D83HJXRL0I4+Qx0IjICeBbANQCKAdxCRMVuq+0FUMLMkwGsBvBksAv1ZH15raJtwfhsJMTJxaD1wECEJ6/IwI3j1UP9mU0V+Pd3DsBml5OPhPBnD30mgApmrmRmK4AVABa7rsDMm5i5bxKV7QDyglumuh67A+sOKPvPr54o3S16YjQQfn9lBpYWq4f633ZW4b7/24NOq1ybVMQ2fwI9F0CVy+/VzjZPvgdg3WCK8teWY/Vo7Og/gZPZZMDlcrq/7hgNhCeuyMCtk5SXsAOAjw6fw9Lnt+FsiyXMlQkROYJ6UJSIbgdQAuD3HpbfS0RlRFRWX6+cSCtQb++pUbRdVTwcqQlyMQs9MhDhtwvSPY5+KT/TisXPbsXnVc1hrkyIyOBPoNcAyHf5Pc/Z1g8RXQnglwAWMbPqxSGZ+QVmLmHmkuzswR20bOnqwYcqo1tunObty4OIdkSEn81Jw8PfdD+M0+tcazeWPr8NK3aeDnNlQmjPn0DfBWAsEY0mIjOAmwGUuq5ARNMA/C96w1x5ymYIrDtQC6ut/4GwoclmGd0SI747dzSeu2064lXONbDaHXjg7QN44K39sPTI9LsidvgMdGa2AVgGYAOAwwBWMnM5ET1KRIucq/0eQAqAVUS0j4hKPTxc0Ly9V9nd8s0pOYgzRv3QeuGnay8aiZXfn4PsVPUrUq3YVYVFz2zF0bMyXYCIDX6lHzOvZeZxzFzEzL91tj3EzKXO21cy83Bmnur8WeT9EQenqqkTO08op8q9cbp0t8SaKfkZePdf5uKi3HTV5cfOtWPRM1vx189OwiFT8Aqdi8rd2VW7qxVtRdnJHt/UQt9yMxKx6r45WDpDfbRst82Bh0vLcfufd6DK7RKFQuhJ1AW61ebA31QOeN04PU+myo1hCXFGPLlkMv7jhos8zuHz2fFGLPyfLfjrZydhl711oUNRF+jry8+ivq3/IBqTgbDEw96ZiB1EhFtnFaB02VyMHaY+tLHDasfDpeW48U+fofxMS5grFCK0oi7QX9t2UtF29aQRGJ6mvCaliE3jR6ShdNnXcNusAo/rfF7VjG/+cSt+9e4BnHc7OU2IaBVVgX7oTCt2nTyvaL9jdtjmAhNRItFsxG9vuAiv3jUTI9PVP+wdDPzf9tO47D83489bT6DbJkMcRXSLqkB/bftJRdv4EamYOXpI+IsRUWHeuGys/9E83DIz3+M6LV09+M0Hh3DFf32Md/ZWS/+6iFpRE+gtnT14R2Xs+XfmjJKDocKr9MQ4/O7GyVh13xyPfesAUH2+Cz9+83N8/b8/xnv7aiTYRdSJmkBftbsKlp7+Z4amJphw/VQZey78c3HhEKz510vx4DXjkRJv8rje8foO/HDFPlz11MdYsfO0dMWIqBE1gf7uPuXe+ZIZeUj28sYUwp3ZZMD35xdh48/mY+mMPHj7clfZ0IEH3j6AS5/YhGc3VaBJDp6KCEdaXe2lpKSEy8rK/F6/o9uGd/bW4NVtJ3HsXO9Fgjf+dD7GZHv+Cv2lM3sHWKWIWDnTgkWY2c4AAA4ySURBVPIwh2tb8eT6I9h01Pfsn2aTAddPzcFts0Zhcl66dPUJTRDRbmYuUV0WLYHeh5mxvbIJu081YdnlY/27kwS6/gQp0PvsqGzE/3z0BbZVNvq1fvHINNwyMx/fnJKDjCRzUGsRwhtdBfqASKDrT5ADvc/2ykb8ceMX+LTCv2A3Gw24YsIwXD8tF/PHyaUPReh5C3TpgBbCxewxQzF7zFAcqG7B8x8fx7qDtfA22MVqd2DdwbNYd/AsUhNM+HrxCFx70QjMvSBLwl2Eneyhi+gUoj10d1VNnXht+yms2HkarRb/r1maEm/C/AuzceWEYbhs3DBkJku3jAgO6XKRQNefMAV6n06rDR/sr8UbO05jX4CXuDMQMDU/A/PGZWPeuGxMzk2HSebtFwMkgS6Brj9hDnRXh2tb8dbuary77wwa2lWvtuhVSrwJs0YPwZyioZg1eigmjEyVgBd+k0CXQNcfDQO9j83uwCcVDXh/3xlsKD+LDuvATkBKNhsxrSATM0ZlYvqoTEzJS5eRM8IjOSgqRAiYjAYsuHAYFlw4DJYeOz4+Vo/1B8/io8Pn0BZAf3uH1Y6tFQ3YWtHwZdvorN4LtlyUm46JuWkoHpkmIS98kj10IYLMamfsOmPFRycs2HjCglMtwZk6ICfFiPFZJlw4NA4XDjVh7NA4FGWakGDSyQlOEfCtKxpIl4sEutDQyWYbtpzqxpbT3dhR0402a/DecwQgP82IMZkmjM40YUyGCYUZRhSmmzAy1QiTIYrCXgLdL9LlIoSGCjNMKMww4Y4pybA5GAfqerCzxoodNd3YdcY6qIBnAKdb7TjdasfmU+5X8gJyU40oSDMhL82IvDQjclONyHH+DE82Is4YRYEvfJJAFyKMTAbCtBFmTBthxvdnpMDBjC+abNhda8XuWis+P9uD4+dtCMY+vM0BnGqxe+zyIQDZyQaMSDZieIoRI5INGJ5ixLAkI7KTDchOMiA7yYghiQYJ/ighgS6EhgxEzj7xONw6KRkA0NrtwMG6HpTX9+BAXQ8O1fegstnm9YzVgWAAdR0O1HU4gLoer+tmJBCGJhoxNNGAoYkGZCYaMCTRgMwE50+iAenxBmQkGJAeT0iLlw8BLUigCxFh0uINuCQ/Hpfkx3/ZZrExjjX24Gij7ct/j5+3oaYtPHO1N1sYzRYbjiuvAOlRchwhPd6AtHhCarwBaebef1PNhBSzASlmQqrZgGQzIdlMSOloQHK8ESnxJiTFm5BsNiLRbITZaJCZLf0kgS5EFEgwESYPN2Py8P5DFzusDpxotqGy2Y7K8zacarbhZIsNJ5vtOG9xeHi08OjoYXT02HGm3d977FBtNRkIiWYjksxGJMYZkWg2ITHOgCSzCQlxvaGfYDL0/hvXezs+zoh4kwEJcUZMGJmGGaMy1Tep1YCJEB0AlkAXIoolmw2YNMyMScOUy9q6HTjdakd1qw3VrXZUt9lxxuWnsUvbwPeXzcFos9gCGtvv6rtzCz0Hus74FehEtBDAHwAYAbzEzI+7LY8H8CqAGQAaAXybmU8Gt1QhRCBS4w2YmG3AxOw41eUWG6Ouw47adjvqOhw4125HXafd2a9uR0OnAw2dDjRpvKc/WLE066XPQCciI4BnAVwFoBrALiIqZeZDLqt9D8B5Zr6AiG4G8ASAb4eiYCFEcCSYCAXpJhSke48Bm4PR1NUb7uctDjR2OdDUZcf5LgeaLYwmiwPNLj8t3Q60dXNQRuoEQ4JJAt3VTAAVzFwJAES0AsBiAK6BvhjAcuft1QCeISJirc5aEkIEjclAGJZsxLBk/4PRwYw2K6O124EWC6PN6kBrtwNtVka7tTfw23t6b7dbGR1WB9opGR3ddnRYbejotqHTaken1Q77IIf3JMTFzsRn/gR6LoAql9+rAczytA4z24ioBcBQAA2uKxHRvQDuBYCCgoIBljwAcgaaEGFlAJDu/MkfxOMwM6x2B7qc4d5ptcPS89W/XT12dDlvW3rssNgcztu9/3bb7CjOSfO8AZ1lQ1gPijLzCwBeAHpP/Q/ntoUQ0YeIEG8yIt5kREaS1tVEPn++i9Sg/4dsnrNNdR0iMqH3g9m/izIKIYQICn8CfReAsUQ0mojMAG4GUOq2TimAf3LeXgJgo/SfCyFEePnscnH2iS8DsAG9wxZfZuZyInoUQBkzlwL4M4DXiKgCQBN6Q18IIUQY+dWHzsxrAax1a3vI5bYFwNLgliaEECIQsTOeJ4SWL18+qOV6FIvPWQitSaCr8Ceg+34A4JFHHvG6vq/l3rY32GBUu38ot9fH13MOhHw49IqU1yFS6hAqmFmTnxkzZnCk6n1ZvC/v+/F3fWbmhx9+OODt+XpsT4/p7f7e6vG1PX8F63GC/ViRwNf/maf1I+V1iJQ6YhV6j12q5mpMBHrfG8LfN1Igge7+47otb8u9bc91HW/L/K3VU5vrMl+h4b5dtecRyHMOhN4CJNDn4++OQ7hESh2xKuYD3Z83hK8w8rRcbT1v68+fP9/n9rxtw/05qD0nX7W6P6Y/QezPdgNZ7kuoPhwiQSCvTaS8DpFSh5BAD3gPx9P63oLRdR3XP3Jfwe5pO96WzZ8/32s4+/Oh4Ok5eHvu7r+r7bX7eq7+8vZhEo0CDcRIDlA9/H9obTD/jzEZ6P7uUavxFWruIegePr4+CFz3tPtuDySE/dljd69ZrR5Pj+vvnr77B5h7HQMJgMHeP5IF+nwC3SEJtUipI5oN5jWMyUB3FcgbIpDuCl8/nvak3Ze7d9e41jl//vwBB7o/z6XvPgP54HD/8HL/IHBtC5T744ZSuPd6BxrokbB3zhzcOiLlOYWbBPog+BPo/oSZpz1cX/dXu4+ndfr+dW/z1p3hbZm35Z6eh2sN/n54ua/rrVtoIP8HoTTQN9dA6wr0fnoOvcEEW7QJ1t93zAe6+x6wL2oB7v6v2n+Ia4i5btvfUFTbq/X0WL5qV6vT1zJP4ez+PAPpivFVrz/PY7B8/b8PdFuxFEahEquv4WCed8wHeqAC3TtVCzFXwezG8fXJ7im0H3YeLPW2nmu72rcGtW0Ntl5v/wfBovZYwdhbitUwGiytvolFEgn0MAkkfJm9h2EgXR2ujzWYYPTVbz/QDye1bbp/eHl6Tv6+UUP1Rvf15gnkzSVhFFyx+qE4mL8XCfQBcg3qQPeyfY03d1/PdXvu2x5ogKgFrNq2+j4E3Nv76vfGfa9fLdAHYrBv9ECCV6sahbyGAyGBPkDue6fuAd9321MYuj+Wr71etS6Rge7tum7TV6B72l6g+rp1XGsc6J5IMN/ovh4rEmqMVfLNJnAS6APkKXDVuhnc29UeK5DQVDuQG2iAqAWs++N7256WgllHqII3Ul4rEVsk0IPMfQ+073dfZ0QONjQlQAZGXjehJ94CnXqXh19JSQmXlZVpsm0hhIhWRLSbmUvUlsl86EIIoRMS6EIIoRMS6EIIoRMS6EIIoRMS6EIIoROajXIhonoApzTZ+FeyADRoXEOgpObwkJrDJxrr1rLmUcycrbZAs0CPBERU5mn4T6SSmsNDag6faKw7UmuWLhchhNAJCXQhhNCJWA/0F7QuYACk5vCQmsMnGuuOyJpjug9dCCH0JNb30IUQQjck0IUQQid0H+hEtJCIjhJRBRE9oLL8TiKqJ6J9zp+7tajTraaXiaiOiA56WE5E9LTzOe0nounhrlGlJl81X0ZELS6v80PhrlGlpnwi2kREh4ionIh+qLJORL3WftYcUa81ESUQ0U4i+txZ8yMq68QT0ZvO13kHERWGv1JFTf7UHVn54WleXT38ADACOA5gDAAzgM8BFLutcyeAZ7Su1a2meQCmAzjoYfm1ANYBIACzAeyIgpovA/CB1nW61TQSwHTn7VQAx1T+PiLqtfaz5oh6rZ2vXYrzdhyAHQBmu63zAwDPO2/fDODNKKk7ovJD73voMwFUMHMlM1sBrACwWOOafGLmLQCavKyyGMCr3Gs7gAwiGhme6tT5UXPEYeZaZt7jvN0G4DCAXLfVIuq19rPmiOJ87dqdv8Y5f9xHYywG8Ffn7dUAriAiClOJqvysO6LoPdBzAVS5/F4N9T/+m5xfp1cTUX54ShsUf59XpJnj/Pq6jogmal2MK+dX/Gno3QtzFbGvtZeagQh7rYnISET7ANQB+JCZPb7OzGwD0AJgaHirVPKjbiCC8kPvge6P9wEUMvNkAB/iq70EEVx70DsHxRQAfwTwrsb1fImIUgC8BeBHzNyqdT3+8FFzxL3WzGxn5qkA8gDMJKJJWtfkDz/qjqj80Hug1wBw/cTMc7Z9iZkbmbnb+etLAGaEqbbB8Pm8Ig0zt/Z9fWXmtQDiiChL47JARHHoDcbXmfltlVUi7rX2VXOkvtYAwMzNADYBWOi26MvXmYhMANIBNIa3Os881R1p+aH3QN8FYCwRjSYiM3oPtpS6ruDWH7oIvX2Ska4UwB3OERizAbQwc63WRXlDRCP6+kSJaCZ6//Y0fcM66/kzgMPM/JSH1SLqtfan5kh7rYkom4gynLcTAVwF4IjbaqUA/sl5ewmAjew86qgVf+qOtPwwabnxUGNmGxEtA7ABvSNeXmbmciJ6FL1Xzi4F8K9EtAiADb0H9e7UrGAnIvobekcqZBFRNYCH0XtABsz8PIC16B19UQGgE8B3tan0K37UvATAPxORDUAXgJu1fsMCmAvgOwAOOPtJAeDfARQAEfta+1NzpL3WIwH8lYiM6P1wWcnMH7i9D/8M4DUiqkDv+/Bm7cr9kj91R1R+yKn/QgihE3rvchFCiJghgS6EEDohgS6EEDohgS6EEDohgS6EEDohgS6EEDohgS6EEDrx/yzgz0DdpS6uAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "x5pvmhcZ2C01"
      },
      "source": [
        "And now we can use `NumericalInverseHermite` to use any `QMCEngine`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 299
        },
        "id": "s9ERuEQwuf5L",
        "outputId": "09b12cd3-0361-45ce-b24a-3dc262e27092"
      },
      "source": [
        "from scipy.stats import NumericalInverseHermite\n",
        "\n",
        "fni = NumericalInverseHermite(dist)\n",
        "sobol_engine = qmc.Sobol(d=1, scramble=True, seed=rng)\n",
        "sample = fni.qrvs(99, qmc_engine=sobol_engine)\n",
        "\n",
        "\n",
        "plot_pdf(dist, sample)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/scipy/stats/_qmc.py:1078: UserWarning: The balance properties of Sobol' points require n to be a power of 2.\n",
            "  warnings.warn(\"The balance properties of Sobol' points require\"\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxb1Z338c/PsuR9SWJnsZ2dhJAEyEZYyw4FShPWQgoFSluGmYFpp8tMO9NCyDBTeGam8zwttBQKQwtpIUApAQIhLYEkQEgcEiAJIXE2HGdzHO+bLOk8f9imsu6VLNmStf3er5dflo+upJ9k+6ujc889V4wxKKWUSn4Z8S5AKaVUdGigK6VUitBAV0qpFKGBrpRSKUIDXSmlUkRmvB64pKTETJgwIV4Pr5RSSWnTpk3HjDGldtfFLdAnTJhAZWVlvB5eKaWSkojsD3adDrkopVSK0EBXSqkUoYGulFIpQgNdKaVShAa6UkqliLjNcklnXp/hs+NteH2GiSV5ODIk3iUppVKABvoQ2rjvOH/8oIZV2w9zrMUNQEFWJnMnDOOKmWO4ek45Tod+aFJKDYwG+hDo8vr48Ytbebay2nJdc6eHtz6t5a1Pa1lWWc1DX53D6KLsOFSplEp22h2MsY4uL3/79CbbMA9Uub+eK36+lnW7jg1BZUqpn//855x00kncdNNNLF++nAceeCDotk8++SR33XVXVB73tttu4/nnnwdg7dq1zJgxg1mzZtHe3j6o+9Ueegy1dnr4+pMb2bD3eNi3Od7q5tb/3cDjt87j/BNHxrA6pdQvf/lL/vznP1NRUQHAggULhryGpUuX8qMf/Yibb7550PelgT5QBzf3u8k9q+rZsNf+HdflgCyH0Oy2njHK6zPcvbSSl24oZdKwgF9R2ewBlatUoprww1dj/hj7HviSpe3OO+9kz549XH755dx+++0MGzaMyspKHnroIZ577jnuu+8+HA4HRUVFrFmzps9tX331Ve6//35efvllSkpKPm9fvHgxu3fvpqqqimPHjvFP//RPfOtb38IYw913382qVasYO3YsLpcLgN/85jcsW7aMlStX8tprr7F06dJBPU8N9BhZva+DFz6xhrkAPzqnkK/OzCUnU3inupPv/7mBo62+Pts1uw3feuU4L36lhMIsHRlTKtoeeeQRXn/9dVavXk1JSQlPPvnk59ctWbKElStXUl5eTkNDQ5/bvfjii/zsZz9jxYoVDBs2zHK/H330EevXr6e1tZXZs2fzpS99ifXr1/Ppp5+yfft2jhw5wvTp07n99tv55je/ybp167jyyiu57rrrBv2c+k0KEXlCRI6KyNYg14uI/FxEqkTkIxGZM+iqklxzp49/fbPR0u7MgF9cNow75uST78rAkSGcOz6bFYtKOa3MZdl+d72H775Rj573VamhdfbZZ3Pbbbfx2GOP4fV6P29/8803efDBB3n11Vdtwxxg4cKF5OTkUFJSwgUXXMCGDRtYs2YNixYtwuFwUFZWxoUXXhiTusPp+j0JXBbi+suBKT1fdwC/GnxZye2Bd5s42OK1tC85v4grp+ZY2ktyHfzmyuFMKHJYrvvz3k6W7xzcjhKlVGQeeeQR7r//fqqrq5k7dy51dXUATJ48mebmZnbu3Bn0tiIS8udY6jfQjTFrgFB79RYCvzPd1gPFIjImWgUmm61H3Sz9uM3Sfs5YFzfOyA16u6LsDB67cjh5Tusv//61TTR1+mxupZSKhd27d3P66aezZMkSSktLqa7unqU2fvx4XnjhBW655Ra2bdtme9uXXnqJjo4O6urqeOuttzjttNM499xzefbZZ/F6vRw6dIjVq1fHpO5ojKGXA/5z8g70tB0K3FBE7qC7F8+4ceOi8NCJ5+GNLZa2XKfw0wuL+32nnjLCyX9fUsydK+r7tNe2+fjZ+mYWn1cU1VqVSgR2Oyzj7Qc/+AG7du3CGMNFF13EqaeeypYtWwCYNm0aS5cu5frrr+fll19m8uTJfW57yimncMEFF3Ds2DF+8pOfUFZWxtVXX82bb77J9OnTGTduHGeeeWZM6pZwxmdFZALwijFmps11rwAPGGPW9fz8F+CfjTEhz14xb948k9QnuLCZ5bK73sPFTx0l8BW959xCbp+VH/Zd3/HKcd7Y09GnLUNg+Q0lzJx1+kCqVUoNgcWLF5Ofn8/3v//9mD2GiGwyxsyzuy4a0ydqgLF+P1f0tKWdRze1WMJ8dF4GN83Mi+h+7j2vkJzMvr15n4Ela5p0B6lSKqhoBPpy4Jae2S5nAI3GGMtwS6o71OzljzusY+ffnJ1PVmZkO0XKCzL5zukFlvYNB928u7tuwDUqpWJr8eLFMe2d9yecaYt/AN4DThSRAyLyDRG5U0Tu7NlkBbAHqAIeA/4uZtUmsMe3tNAVsN+yKEtYNDP4jtBQbp+VxwmBBxUB/7Nqp/bSlVK2+t0paoxZ1M/1Bvj7qFWUhNq7fDyz1do7v/XUPPJcA/sQ5HQI/3hGAX//Wt8dpJX761m76xjnTrU96bdSKo3pIYhR8MaeDlq6+vaaczKF206NbOw80OUnZDNthPU992faS1dK2dBAj4LnbQ7x//LUHIbnWA8UikSGCN85wzqWvqW6gbd31g7qvpVSqUcDfZAOt3h5p7rT0n7tSdYjQgfii5OymVHqtLQ/vm5vVO5fKZU6NNAH6cUdbfgCRj/GFjps12YZCBHh2/Otc9jX7jrGjsNNUXkMpVRq0EAfBGOM7YqK10zLISOK6zdcPCmbScXW4ZvH12ovXSn1Vxrog/DR0S6q6j2W9mtPGthUxWAyRLh9trWX/tKWgxxt7rC5hVIqHWmgD8KfPrX2zueXuRhXFP1l5q+dlkNxdt9ev9vr4+n39kf9sZRSyUkDfYCMMbyx29o7jtbO0EA5TvslBJ5av5+OLutSvUqp9KOBPkDbj3moae4bpA6BL06OTaAD3HJKHk5H3156fVsXr289HLPHVEolDw30AXpjt3W45fRyF8XZsXtJR+U7uPKUMkv70+t12EUppYE+YIHL2wJcOjk75o978xnWdeQr99frFEallAb6QFQfb+OTY9bZLRdPjH2gzxk3jGmjrUePLl3/WcwfWymV2DTQB2DV9iOWthmlTioKoz+7JZCIcNPp1l76i5traO20vskopdKHBvoAvLHduhPy0kmx7533ump2ObmuvgcatXR6eGnLwSGrQSmVeDTQI1Tf6mbDXus5s4di/LxXQbaThbPKLe3PbtRhF6XSmQZ6hNZWHbNdu8VumdtYsht2+fBAo+4cVSqNaaBHaN0u67K1F07IRqK4dks4ZpYXMaOs0NL+XOWBIa1DKZU4NNAjYIxh7a5jlvZzx2fFoRq44bSxlrYXN9fg9vhstlZKpToN9Ajsrm3hUGPf+efODDijPDpL5UZqwalluDL7/gqPt7r5yyfWWThKqdSngR6BNTutvfM5Y1wDPm/oYBXnuvjijNGW9mWV1XGoRikVbxroEVhrM35+7rj4DLf0umGeddjl7Z21HGq0Lk2glEptGuhh6vR4Wb/HOl3xC3EO9LMmj6C8uO+CYD7TPZaulEovGuhh2rS/nvaAZWqHZWfYnu9zKGVkCNfNrbC0v7DpAMYYm1sopVKVBnqY1tnMbjl7rAtHxtBOV7Rz7RxroO+ubeXDA41xqEYpFS8a6GFaV2UzXXHc0B0dGsq4EbnMnzDc0v7CJp2TrlQ60UAPQ1NHF1trrL3dc+I8fu7v2rnWpQBe/uggnR49m5FS6SKsQBeRy0TkUxGpEpEf2lw/TkRWi8hmEflIRK6IfqnxU7nvuOVw/wlFDsoKHPY3iIMrTh5DtrPvr7OhrYvVO47GqSKl1FDrN9BFxAE8DFwOTAcWicj0gM1+DCwzxswGbgR+Ge1C48ludsvp5YnTO4fuBbvs5qQ/v0lnuyiVLsLpoc8Hqowxe4wxbuAZYGHANgboXVikCEipdVzX76mztJ1REZ+jQ0O5xmbn6FufHqWupTMO1Silhlo4gV4O+B96eKCnzd9i4GYROQCsAO62uyMRuUNEKkWksrbWepBOImoOMn6eaD10gHNOKGFUYd+6PD7DKx8dilNFSqmhFK2doouAJ40xFcAVwFMiYrlvY8yjxph5xph5paWlUXro2KrcV28ZPx83PDehxs97OTKEq2zWSf+jHmSkVFoIJ9BrAP/jyyt62vx9A1gGYIx5D8gGSqJRYLzZDrdMsk4RTBRXz7EG+ofVDeyubYlDNUqpoRROoG8EpojIRBFx0b3Tc3nANp8BFwGIyEl0B3pyjKn0Y73N2YnOmDQiDpWEZ9roQk4aY10n/cUPtJeuVKrrN9CNMR7gLmAl8Ands1m2icgSEVnQs9n3gG+JyIfAH4DbTAocdx50/DyBAx3gWpte+ouba/AFjh0ppVJKWOdNM8asoHtnp3/bPX6XtwNnR7e0+KvcX483IATHDs/pXgyrLU5FhWHBqWX8x4pP+oz91zS0s2Hf8YT+dKGUGhw9UjSEjXbDLRMTPxBHFmZzzhTrTmcddlEqtWmgh1C5r97SdtrExN0h6s9u2GXFx4fo6NKlAJRKVRroQXR6vGw50GBpP81mEaxEdOn00eS5+k6tbO708Gc9PZ1SKUsDPYitNY2Wky2X5LuYMCI3ThVFJsfl4LKZYyztOuyiVOrSQA/Cbrhl3vjhiMR//fNwXWMz7PLWzlqO6VIASqWksGa5pKONdoE+YVgcKglwcHPYm56RbRiTn8Ghlr9+0vD6DC+v2cDXZ+VH9rhlsyPbXik15LSHbsPnM2zab53hMi9Jxs97OTKEhSdah4he3KEnkFYqFWmg29hzrIX6tq4+bdnODGaUWY/ATHTXTMuxtH10tIuq4102WyulkpkGug278fPZY4fhdCTfyzV1hJOZNiey/qP20pVKOcmXUEMgYcfPB+iak6y99Bd3tOFL/tUZlFJ+NNBtVKbA+Lm/BVNzcARMzjnU4mP9AXd8ClJKxYQGeoCjzR3sr+u7UEuGwJxxxXGqaPBKch2cP956Qo4XdiTwgjRKqYhpoAf4YL/16NCpowooyLaOQyeTa06yznZ5raqDti6fzdZKqWSkgR5g82fW8fO545N3/LzXRROzKXD1HXdp6zKs3N0Rp4qUUtGmgR7ggxQN9OxM4cqp1p2jz3+iwy5KpQoNdD9uj48PD1hPaDFnXPIHOsC1NnPS3612U9PsiUM1Sqlo00D3s/1Qk2VBrhF5LsYnyYJc/Zk7xsWEor4rMBr0yFGlUoUGup9N+20OKBo3LKkW5ApFRLhuuvXN6fntbaTAGQOVSnsa6H5Sdfzc39XTcgh8e9rX6KXykM5JVyrZaaD7+cCmh57M88/tlBdkcvZYl6X9+e067KJUstNA73GwoZ1DjX2n8GVmCKdUpFagA1xnMyf91V3tOiddqSSngd7DbrhlelkhOQGncUsFX5ycTb6z78BLS5fhtSqdk65UMtNA72F3hGiqTFcMlOPMsJ2T/uw2nZOuVDLTQO9h10Ofk2I7RP1dbzPbZcNBN3sbdE66UslKAx3o6PKy7aDdAUWpN37ea85oJ5OHWc9AuEx76UolrbACXUQuE5FPRaRKRH4YZJuviMh2EdkmIr+Pbpmxte1gI13evvOwRxZkUV5sHZZIFSLCDTOsvfQXdrTh8emcdKWSUb+BLiIO4GHgcmA6sEhEpgdsMwX4EXC2MWYG8J0Y1Bozmz+zHz9PlQOKgrlmWg6ZAX8BR1t9vLWvMz4FKaUGJZwe+nygyhizxxjjBp4BFgZs8y3gYWNMPYAx5mh0y4wt+/Hz1B1u6VWS6+DiidmW9me367CLUskonEAvB6r9fj7Q0+ZvKjBVRN4RkfUiclm0ChwKdjNcZqfoDJdAdsMub+7t4EiLNw7VKKUGI1o7RTOBKcD5wCLgMRGxdHFF5A4RqRSRytra2ig99OAcbGjncJP1gKKTy4viVNHQOndcFqPz+v4ZeA0s0166UkknnECvAcb6/VzR0+bvALDcGNNljNkL7KQ74PswxjxqjJlnjJlXWlo60Jqjym78fEZZIdnO1DugyI4jQ/iKTS/9mW1teHXnqFJJJZxA3whMEZGJIuICbgSWB2zzJ7p754hICd1DMHuiWGfM2I2fp8twS68bZuSSEbD/t6bZy5rPdOeoUsmk30A3xniAu4CVwCfAMmPMNhFZIiILejZbCdSJyHZgNfADY0xdrIqOJvtAT/0dov7KCzJtTyL9+6067KJUMrEeWWLDGLMCWBHQdo/fZQN8t+craXR6vGyrabK0p+oh/6F8dWYebwZMV3xzbweHW7yMzk+P4Selkl1aHym67WATbm/fFQZLC7KoGJa6BxQFc/6ELMbkW3eOPqNHjiqVNNI60IOtf57qBxTZycwQbpiRZ2n//dZWy1G0SqnElN6BbndAURoOt/S6cUYujoD3sqOtPt7Yo8vqKpUM0jvQ7ZbMTeEVFvszOt/BpZOsR47+9sPWOFSjlIpU2gZ6uh9QFMzXTrEOu2w46GbHYevOY6VUYknbQLcbbkmnA4qCObPCxZTh1slPT723Pw7VKKUikb6Bnsbrt4QiItxi00t/cXMNje1dcahIKRWutA30TTY99LlpPH7u7+ppOZZzjra5vSzbWB3kFkqpRJCWgd7R5WW73RmKNNAByHdlcO1J1vVdnnx3H56AeftKqcSRloG+tcZ6hqJRhVmUFVlneKSr22blETgbv6ahnZXbjsSlHqVU/9Iy0IPNP0/HA4qCmVicyUUTreu7PL4uKdZcUyotpWeg280/1x2iFrfPyre0ffBZA5tt3hCVUvEX1uJcCe3g5og2N8awaa/1DHlz8mrBZlw9nZ1Z4WJaSSY7jnn6tP9m3V4e/qq+ASqVaNKuh17d5KW2re+OPVcGzCh1xqmixCUifMOml/7ax4fYX6dHjyqVaNIu0CsPui1tJ49ykp2p4+d2FkzNoSS375+Jz8Bja3UsXalEk36Bfsga6HPHuOJQSXLIyhRun2U90Oi5ygMca9EzGimVSNIu0DdpoEfsppPzyM/qu7ul0+PjyXf2xacgpZSttAr0xk4fO+s8lnYN9NCKsjK46fRxlvbfvbePlk7r66mUio+0CvTNh9wEnqphYrGDktz0XpArHLefMxGXo++fS1OHh6fX66JdSiWKtAp0HW4ZuFGF2Vw9u9zS/tiaPbS7vXGoSCkVKK0C3W6H6DwN9LDdef5kMgImA9W1uln6vvbSlUoEaRPoXV7DlsPW5V/nlWmgh2tiSR4LTi2ztP96zR46urSXrlS8pU2gf3Ksi3ZP3xH04mxh0rDkP1h2KN114QkELnlT29zJMxs+i09BSqnPpU2g284/H+0iQxfkisgJIwv40sljLO2/enu39tKVirO0CfSNNdZAn6Pj5wNy94VTLG1Hmjp1xotScZYWgW6MYYPNIf+nl2ugD8SJowu44uTRlvaHV1fR3KGnqVMqXsIKdBG5TEQ+FZEqEflhiO2uFREjIvOiV+Lg7a73UNfed0GuLAecPFIDfaD+8eKplhkv9W1dPLFuX1zqUUqFEegi4gAeBi4HpgOLRGS6zXYFwLeB96Nd5GDZ9c5nj3aRpQtyDdiUUQVcPbvC0v7Y2j3Ut1pfb6VU7IXTQ58PVBlj9hhj3MAzwEKb7f4NeBDoiGJ9UbHBZvx8vg63DNp3Lp6C09H3TbGl08NDq6viVJFS6S2cQC8H/E/3fqCn7XMiMgcYa4x5NdQdicgdIlIpIpW1tbURFztQG2166PPLrKdXU5EZOzyXRfPt13jR9dKVGnqD3ikqIhnAz4Dv9betMeZRY8w8Y8y80tLSwT50WA40eahp7judLjMD5ozRE1pEw10XnECOs+9aOF1ew4Ov74hTRUqlr3ACvQYY6/dzRU9brwJgJvCWiOwDzgCWJ8qOUbvhlpmlTnKdaTHBJ+ZGFmbzrXMnWdpXfHyYTfuPx6EipdJXOKm2EZgiIhNFxAXcCCzvvdIY02iMKTHGTDDGTADWAwuMMZUxqThCOl0x9v7m3EmUFliHsP7tlU8wJnB9S6VUrPQb6MYYD3AXsBL4BFhmjNkmIktEZEGsCxysDTXWs+rML9fx82jKy8rk+5dOtbRvqW7gjx/U2NxCKRULYY07GGNWGGOmGmMmG2P+vaftHmPMcpttz0+U3vnRVi97GvqOnwu6wmIsXDd3LNNGF1jaf/raDpr0YCOlhkRKDyS/d8DaOz+xJJOi7JR+2nHhyBB+cqXl8ASOtXTyf1ftikNFSqWflE629w5Yx8/PqtDhllg5+4QS2yUBfvvePnYcbhr6gpRKMym9duy71dYeugb6AB3cHNZm/zrPx+pPpM9SxV6f4cfPvs+y60ZEvrpl2ezItlcqjaVsD726ycNnTX3HzzNEjxCNtfKCTO46Ld/SXnnIze+3tsWhIqXSR8oG+nvV1uGWU0Y6KcxK2aecML45O5+JxdYTbz/4ThOHW3TNdKViJWXT7V2bHaJnjdXhlqGQlSn8+wXFlvZmt+HetxvjUJFS6SElA90YYx/oOn4+ZM4am8VXpuda2lfu7uCVne1xqEip1JeSgb673sPR1r7rn7syYK6u3zKk/uWcQkpyrX9iP36rgaOtOvSiVLSlZKDbTVecPcZFjq7fMqSKszNYfG6Rpb2hw/DDvzTosgBKRVlKJtw7Ol0xYXxpSjaXT862tL+5r5Nnt+msF6WiKeUCvctr7AN9rE5XjAcR4f4LiijJsf6pLVnbxO56TxyqUio1pVygbz7sptnd96N8gUuYNUoDPV5G5Dr46UXWWS9tXYa7XztOh0eHXpSKhpQL9Lf3W3vnZ4/NspwqTQ2tSyZlc/30HEv79mMeHnhHlwVQKhrSItDPG6/j54ng3nOLbA84evLDVl6v0qmMSg1WSgV6bZuXrbXWpVrP1UBPCPmuDH5x2TBcNn9131vVQNVxXWZXqcFIqUBfa9M7nzI8k/KClF6DLKnMHOniR+cUWtpbuwx3vFpPc6fP5lZKqXCkVKC//ZkOtySD207N44s2Uxn31Hv43qoGfDo/XakBSZlA9/oMa3T8PCmICP91cTGTbMbT39jTwX++2xyHqpRKfikT6B8f7aK+o+/H9ZxM4bQyDfREVJCVwaNXDifPaZ199KtNLSzTg46UiljKBPpf9nZY2s6scJGdqdMVE9UJw5389yXW+ekA/7K6wfYAMaVUcCkT6Ct3WwP9ggnWcVqVWC47IYfvnWE9ubTHB3/zynG21uhyu0qFKyUCfW+Dh53HrYeQXzJJAz0Z3HVaPtdMsx501NJluPWJDew91hqHqpRKPikR6Kv2WHvnp45yMjrfutNNJR4R4acXFjO/zLo8Q12rm689/j4HG/TAI6X6kxKB/obNcMul2jtPKlmZwqNXDmfqcOsxAwfq21n02HoON1p/z0qpv0r6QK9t87LpkHX9c7t5ziqxFWdn8LurRlBeYP1ktb+ujUWPredIk4a6UsEkfaD/ZU8HgYehTCp2MHmYHh2ajEbnO3jqqhGMsFlud++xVm749XscqNcpjUrZCSvQReQyEflURKpE5Ic2139XRLaLyEci8hcRGR/9Uu29YTN+fsnkHER0umKymjQsk6euGkFxtvV3uK+ujesfeY+qoy1xqEypxNZvoIuIA3gYuByYDiwSkekBm20G5hljTgGeB/5PtAu109jWxVqbw/11/Dz5TS918vRVJRTlWM8De6ixgxt+/R4fVjfEoTKlElc4PfT5QJUxZo8xxg08Ayz038AYs9oY0/s5eD1QEd0y7a3YeoiugLWcSnMzmD1aTwadCmaOdPL0N06nONf6+6xrdXPDo+/xxrbDcahMqcQUTqCXA9V+Px/oaQvmG8BrdleIyB0iUikilbW1teFXGcRLW2osbVdOzSFDh1tSxskVRTx7x5mUFliXcOjo8vE3T2/if9/ZqyecVooo7xQVkZuBecB/2l1vjHnUGDPPGDOvtLR0UI91qLGd9/cet7RfdaL1ABWV3E4cXcDzd55JxTDr79YYuO/l7fzzCx/R6fHGoTqlEkc4gV4DjPX7uaKnrQ8RuRj4V2CBMSbmi3C8/OFBAjtlE4sdnDJSh1tS0fgRebzwt2cxfYx1LXWAZZUHuPFRnauu0ls4gb4RmCIiE0XEBdwILPffQERmA7+mO8yPRr9Mqz9tPmhpW3hirs5uSWGjCrNZdueZnDfV/tPd5s8auOLna3l75+CH85RKRv0GujHGA9wFrAQ+AZYZY7aJyBIRWdCz2X8C+cBzIrJFRJYHubuo2HWkme2HrCcWXqjDLSkvPyuTx2+dx9fOsJ8Ze7zVza1PbODB13fg9ujZj1R6CevoG2PMCmBFQNs9fpcvjnJdIb3wgXVn6KmjnEws1oOJ0kGmI4N/u2om08sKueelrXR5rTtEf/XWbtbuquV/vjKLKaOsqzkqlYqS7kjRTo+X5yqrLe3aO08/i+aP4w/fOoORNjNgALbWNPGlX6zj0TW78Xi1t65SX9IF+utbD1PX2nftFmcGLJiqgZ6O5k0Yzopvf4EvTCmxvd7t8fEfK3Zw9S/fZftB6zCdUqkk6QJ96frPLG2Xn5BDSa4ulZuuSvKz+O3X5/ODL56II8N+p/jHNY18+aF1LHl5O00dXUNcoVJDI6kCfeeRZjbss849v/nk3DhUoxJJRobw9xecwIt/dxaTSvNst/H6DE+8s5cL/+ttnqusxuvTg5FUakmqQF+6fr+lberwTE6zOTGCSk+nVBTz6t1f4OtnTyDYDNZjLZ384PmPuPIX61i369jQFqhUDCVNoLd2evijzeyWm07O07nnqo8cl4N7vzyD5+88iykj84Nu98mhJm5+/H2++th6Nu23fvJTKtkkTaC//OFBmjv7njc0x+ngaptzUSoFMHf8MF75h3P4/qVTyXYG/1N/d3cd1/7qPW55YgMbbJaTUCpZJE2g/2GDdWfoVbPLKMxKmqeg4iAr08FdF07hz989j8tmjA657ZqdtXzl1+9x/SPvsmr7EXw6xq6SjMRrlbp58+aZysrKsLc/UN/GMxuqeWZjNcdaupeKeeXuc5gpe2JVokpB79d08h/rmvjwSP8zXSaW5HHbWRO4Zk45Bdm6Rs80KU8AAA42SURBVJBKDCKyyRgzz/a6ZAn0Xm6Pjze2H2bD3uMsWTgTDm6OQXUqlfmM4dVdHfzP+ib2NPS/QmOuy8HCWeV8df44ZpYX6j4bFVcpFegWGuhqgDw+w58+becXG5rZ3xje0rsnjSnk+rkVfPnUMts12pWKNQ10pULw+Ayv7mrnl5UtfFrn6f8GgCNDOPuEEr58yhgunT6aIpuzKikVCxroSoXBZwxr9nfy+JZW23PVBuN0CGdNLuGLM0Zz8UkjGVmo57RVsaOBrlSEdtV18fTHbbzwSRstXZH9j5w6ysn547M4f0I2p4x0Bl2OIKiy2ZFtr9KKBrpSA9Tq9vHyznaWbW/jg8ORrwFTlCWcVZHFWWOzOLPCxeRhmf3vVNVAVyGECnRdQFypEPJcGdw4M48bZ+ZRdbyLF3e086dP26lpDm8namOn4bXdHby2u/vUeKW5GcwvczG3zMW8MS5OKnHidOisGRUd2kNXKkI+Y/jgkJtXdnXwWlU7R1oHvtZ6lgNOHuli1mgnJ490cvJIFxNOmkdGpMM0Km3okItSMeIzhs2Hu1i1p4NVezrYXR/eLJlQ8lwOThpTyPSyQk4cXcC00QVMGVVAoR7cpNBAV2rIfNbo4a19nby1v4P1NW7aItyhGsqYomxOGJnP5NJ8Jo/MZ1JJHhNL8hhdmK09+jSiga5UHLi9hi2H3bxT3cn7NW4+OOzGHd7Qe0SyMjMYNzyX8SNyGTs8l7HDcqkYlkPFsFzKi3MozAljR6xKGrpTVKk4cDmE+eVZzC/vPqK0w2P4+KibjQfdVB50s+VwF8c7Bn+u006Pj11HW9h1tMX2+vysTMYUZTOmOIcxhdmMKspmdGE2IwuyGFWYzcjCLEbkuch06EJ3yU4DXakhkp0pnFaWxWll3QFvjOGzRi9bjrj5+GgXHx3pYnttV8Tz3vvT0ukJGfgAIjA810VJfhYj8l2MyO8O+ZJ8F8PzshiW62RYnothuS6G5TopynWSlamnfUw0GuhKxYmIML44k/HFmSw8sbvNZwzV2Sey7WATOw43s+NQE7uOtrC/rpVYruZrDNS1urtPwH4kvNvkOB0U5zopynFSmO2kMMdJYU5m9+XsTAqyneRnZ1KQnUl+Vvf3vKxM8lzdP+dlZeLK1E8F0aSBrlQCyRBhfOdOxo+AK0YAMxxAER2eQvbUe9jd87Wn3sO+Bg97Gzw0u+OzH6y9y0t7o5dDjR0Dvg9XBuRmO8lzZZLjcpDrcpDj7P6e68oku+dyjstBttNBtjODHOdfL2dndl+eVJrH+BH255INKp7732J08JgGulJJIDtTmF7qZHpp36mLxhjqO3zsb/TyWaOH6iYv1Y1eDjR7ONjspabZG5MdsdHi9oG7rYuGtsiPwvX37Yum8I+XTI1SVckrrEAXkcuA/wc4gN8YYx4IuD4L+B0wF6gDbjDG7ItuqUqpQCLC8BwHw3MczB5tPVm6zxjq2n0cavZysMXLkRYfh1u8HGn1crTVy5FWH8fafNRHYedsPGU7dTwfwgh0EXEADwOXAAeAjSKy3Biz3W+zbwD1xpgTRORG4EHghlgUrJQKX4YIpbkOSnMdnDIq+HZur6GuzUddu5faNh917T6Ot3d/r++5XN/hpb7D0NDuo7HThzeBztAX6pyx6SScHvp8oMoYswdARJ4BFgL+gb4QWNxz+XngIRERE69J7kqpiLgcwpgCB2MKwuvpGmNocRsaOn00dhiaOn3dX25DU4ePZrePZreh2e2jpbPnu9vQ2mVocfto7TK0uk3U3hRytIcOhBfo5UC1388HgNODbWOM8YhIIzACOOa/kYjcAdwBMG7cuAGWHEBXplNqyAlQ0PM1doD3YYyh0+Ojze2ltdNDm9tLm9tDu9vbfbnLS4fbS3tX988dXf5fPjo8f71cMSw38gJSMDuGdKeoMeZR4FHoPlJ0KB9bKZVYRKRntoqD4XnW8X8VuXAGnmro+yZc0dNmu42IZAJFdO8cVUopNUTCCfSNwBQRmSgiLuBGYHnANsuBW3suXwe8qePnSik1tPodcukZE78LWEn3tMUnjDHbRGQJUGmMWQ48DjwlIlXAcbpDXyml1BAKawzdGLMCWBHQdo/f5Q7g+uiWppRSKhI6edPP4sWLE+I+lFJqINIy0IOF7n333dfvNv0J5z7OP//8Ad23/32GW99g32ACb2/3+EP9JqZvmkrZS/4TXAyAiOD/vBcvXszixYv7tAduM5D79r/c+xiDuW//24Z7H8G286/H7rpe9913X5/b2z3+YJ7PQAz14ymVSEKd4AJjTFy+5s6da+Kl+2l3u/feew1g++W/TSjh3Eewy3b3Haqt97b+99FfXXb3CQR9XsGeQ7DHD6cW/5oGcp1dDYN5LKWSFd2TUWxzNS0C/d577w0auv6hFyyUwwmQcO4j3DeOwMc777zzQt6+9/n530+w7QKvD/ZcBvrVX2CHei1DvdH195xC3V86BXs6PddkNpjfU9oHul0vM1hAnHfeeZbbhAqOUD3ncILPrqdrV2+oWgNvE/jc+nvO9957b79vGoGvy0B66P29OYZ7n5E8XiTbp4J0eq7JbDC/Jw10m4C0+zmc3rp/iIYaaonlV2APO3BoZaD36/962L0h2H23ez0DhXoj6W+4yu6NNNTjDaQnH+q+orHNUNJATw4a6EHYfTTv/R4qRPz110u26+n6h5HdG4H/tv73FxjGod44QgVdfz3qYG8Awe4vWKAH+wQSOMQT7HcTScCGqj/wfvsTye8/mHD+6RIhQKP5JqZiJ1q/p5QO9MCAsfsHC/ZP13ubYOEYSWgGBmfv4waGe6hfaLDeaSQ1hLrvcO/vvPPOC2sYJtxhqMDH7e936X95MIEZ7D7Cvc9kCXR/8apH3zwiM8i/69QP9FABEKpH6C9Y4I4fP35Qodp73/61BAZ972X/XmSoN4FQYWwXuOEMJ4X6p4wkYIP9Tvp7jHA+XUUi2Jt8qOcQzmNHo/cfK/EK9ER7Y0t0Guh+wg25wICweWGCttl9j/ZXOM9joI/tP9bv/9xC1RFKuIEerbAL9YkrUpEO/fSK5M0rktvEUrzeUOL9vJPNYH5PKRfovSIdnw3nNoFB4n+73p+DhWGwMA6sI9Jed+Bt/O8nsNbA5x/qseyeX6jXOtR2kb7JhisWQRHufSZjoA+laH2SUpFJ2UDveXK2343p/w8u2D9ff73MYEEbToj198YSzldg7cHGugP/sfxvE26QD0So30mkYllfNB57oL3/VJNOb2TxltKBHhjOwf6J7P7gBtID630M/yEEuzeIYPXYhazd7YKFd6iecn9DJ3b3EwuBzz3R/tmH4rmnm3R93vGQ0oHeq79/0mDB3J9I/lDtgrk/gSHsH4IDCcRQ2w5Vj7G/N7FUlq7Blk6/43hLi0Dvz0D/4CK5XeDY9WAeI5r3pYaOvv4q1kIFelqutqiUUskq1GqLabkeulJKpSINdKWUShEa6EoplSI00JVSKkVooCulVIqI2ywXEakF9sflwbuVAMfi+PiR0npjS+uNLa03esYbY0rtrohboMebiFQGm/qTiLTe2NJ6Y0vrHRo65KKUUilCA10ppVJEOgf6o/EuIEJab2xpvbGl9Q6BtB1DV0qpVJPOPXSllEopGuhKKZUiUj7QReQyEflURKpE5Ic2198mIrUisqXn65vxqLOnlidE5KiIbA1yvYjIz3uey0ciMmeoawyop796zxeRRr/X9p6hrjGgnrEislpEtovINhH5ts02CfMah1lvwrzGIpItIhtE5MOeeu+z2SZLRJ7teX3fF5EJQ1/p57WEU2/C5ENYgq2rmwpfgAPYDUwCXMCHwPSAbW4DHop3rT21nAvMAbYGuf4K4DVAgDOA9xO83vOBV+L9uvrVMwaY03O5ANhp8/eQMK9xmPUmzGvc85rl91x2Au8DZwRs83fAIz2XbwSeTfB6EyYfwvlK9R76fKDKGLPHGOMGngEWxrmmoIwxa4DjITZZCPzOdFsPFIvImKGpziqMehOKMeaQMeaDnsvNwCdAecBmCfMah1lvwuh5zVp6fnT2fAXOulgI/Lbn8vPARSIiQ1RiH2HWm1RSPdDLgWq/nw9g/w9xbc/H6+dFZOzQlDYg4T6fRHJmz0fa10RkRryL6dXzUX823b0yfwn5GoeoFxLoNRYRh4hsAY4Cq4wxQV9fY4wHaARGDG2VfxVGvZA8+ZDygR6Ol4EJxphTgFX8tfegBu8DutedOBX4BfCnONcDgIjkAy8A3zHGNMW7nv70U29CvcbGGK8xZhZQAcwXkZnxrKc/YdSbVPmQ6oFeA/i/o1b0tH3OGFNnjOns+fE3wNwhqm0g+n0+icQY09T7kdYYswJwikhJPGsSESfd4bjUGPNHm00S6jXur95EfI17amkAVgOXBVz1+esrIplAEVA3tNVZBas3yfIh5QN9IzBFRCaKiIvunTDL/TcIGB9dQPc4ZaJaDtzSMxPjDKDRGHMo3kUFIyKje8dHRWQ+3X9vcfvn7anlceATY8zPgmyWMK9xOPUm0mssIqUiUtxzOQe4BNgRsNly4Naey9cBb5qevY9DLZx6kywfyIx3AbFkjPGIyF3ASrpnvDxhjNkmIkvoPnP2cuAfRGQB4KF7B99t8apXRP5A96yFEhE5ANxL944ajDGPACvonoVRBbQBX49Ppd3CqPc64G9FxAO0AzfG65+3x9nA14CPe8ZNAf4FGAcJ+RqHU28ivcZjgN+KiIPuN5ZlxphXAv7fHgeeEpEquv/fboxTrRBevQmTD+HQQ/+VUipFpPqQi1JKpQ0NdKWUShEa6EoplSI00JVSKkVooCulVIrQQFdKqRShga6UUini/wPNnX1hLbWGkAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nu2-pSO22NSY"
      },
      "source": [
        "There is also `NumericalInverseHermite.rvs` which uses MC."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-fOq96Ax6B7k"
      },
      "source": [
        "## Integration convergence\n",
        "\n",
        "We all like convergence plots don't we?\n",
        "\n",
        "In the following, I evaluate a squared sum in 3 dimension as an example."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 297
        },
        "id": "VDPWdjDqwFnl",
        "outputId": "30e1e0e5-3ee9-4c2f-d15b-e78f48663282"
      },
      "source": [
        "from collections import namedtuple\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.stats import qmc\n",
        "\n",
        "n_conv = 99999\n",
        "min_m, max_m = 4, 11\n",
        "ns_gen = 2 ** np.arange(min_m, max_m)\n",
        "\n",
        "\n",
        "def art_2(sample):\n",
        "    \"\"\"3-D squared sum.\n",
        "\n",
        "    True value 5/3 + 5*(5 - 1)/4\n",
        "\n",
        "    Art B. Owen. On dropping the first Sobol' point. arXiv 2008.08051, 2020.\n",
        "    \"\"\"\n",
        "    return np.sum(sample, axis=1) ** 2\n",
        "\n",
        "\n",
        "functions = namedtuple('functions', ['name', 'func', 'dim', 'ref'])\n",
        "case = functions('Art 2', art_2, 5, 5 / 3 + 5 * (5 - 1) / 4)\n",
        "\n",
        "\n",
        "def conv_method(sampler, func, n_samples, n_conv, ref):\n",
        "    samples = [sampler(n_samples) for _ in range(n_conv)]\n",
        "    samples = np.array(samples)\n",
        "\n",
        "    evals = [np.sum(func(sample)) / n_samples for sample in samples]\n",
        "    squared_errors = (ref - np.array(evals)) ** 2\n",
        "    rmse = (np.sum(squared_errors) / n_conv) ** 0.5\n",
        "\n",
        "    return rmse\n",
        "\n",
        "\n",
        "# Analysis\n",
        "sample_mc_rmse = []\n",
        "sample_sobol_s_rmse = []\n",
        "sample_sobol_rmse = []\n",
        "\n",
        "for ns in ns_gen:\n",
        "    # Monte Carlo\n",
        "    sampler_mc = lambda x: rng.random((x, case.dim))\n",
        "    conv_res = conv_method(sampler_mc, case.func, ns, n_conv, case.ref)\n",
        "    sample_mc_rmse.append(conv_res)\n",
        "\n",
        "    # Sobol'\n",
        "    engine = qmc.Sobol(d=case.dim, scramble=False)\n",
        "    conv_res = conv_method(engine.random, case.func, ns, 1, case.ref)\n",
        "    sample_sobol_rmse.append(conv_res)\n",
        "\n",
        "    engine = qmc.Sobol(d=case.dim, scramble=True)\n",
        "    conv_res = conv_method(engine.random, case.func, ns, n_conv, case.ref)\n",
        "    sample_sobol_s_rmse.append(conv_res)\n",
        "\n",
        "sample_mc_rmse = np.array(sample_mc_rmse)\n",
        "sample_sobol_rmse = np.array(sample_sobol_rmse)\n",
        "sample_sobol_s_rmse = np.array(sample_sobol_s_rmse)\n",
        "\n",
        "# Plot\n",
        "fig, ax = plt.subplots()\n",
        "\n",
        "\n",
        "# MC\n",
        "ratio = sample_mc_rmse[0] / ns_gen[0] ** (-1 / 2)\n",
        "ax.plot(ns_gen, ns_gen ** (-1 / 2) * ratio, ls='-', c='k')\n",
        "\n",
        "ax.scatter(ns_gen, sample_mc_rmse, label=\"MC\")\n",
        "\n",
        "# Sobol'\n",
        "ratio = sample_sobol_rmse[0] / ns_gen[0] ** (-2/2)\n",
        "ax.plot(ns_gen, ns_gen ** (-2/2) * ratio, ls='-', c='k')\n",
        "\n",
        "ratio = sample_sobol_s_rmse[0] / ns_gen[0] ** (-4/2)\n",
        "ax.plot(ns_gen, ns_gen ** (-4/2) * ratio, ls='-', c='k')\n",
        "\n",
        "ax.scatter(ns_gen, sample_sobol_rmse, label=\"Sobol' unscrambled\")\n",
        "ax.scatter(ns_gen, sample_sobol_s_rmse, label=\"Sobol' scrambled\")\n",
        "\n",
        "ax.set_xlabel(r'$N_s$')\n",
        "ax.set_xscale('log')\n",
        "ax.set_xticks(ns_gen)\n",
        "ax.set_xticklabels([fr'$2^{{{ns}}}$'\n",
        "                    for ns in np.arange(min_m, max_m)])\n",
        "\n",
        "ax.set_ylabel(r'$\\log (\\epsilon)$')\n",
        "ax.set_yscale('log')\n",
        "\n",
        "ax.legend(loc='lower left')\n",
        "fig.tight_layout()\n",
        "plt.show()\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deVxV1frH8c8CEQRxApxBcAxFTMQBsrLBodQ0K8vSokFLc0hFb976JXUtK3FAs4xyrmuDeb0q3swG08TZnAJHRHBKUVEgkGn9/mCIeTxwBp7363VedTbn7P1sQb+svddZj9JaI4QQQpgaK2MXIIQQQhRFAkoIIYRJkoASQghhkiSghBBCmCQJKCGEECaplrELqA7Ozs7a3d3d2GUIIYQowoEDB+K01i4Ft9eIgHJ3d2f//v3GLkMIIUQRlFLnitoul/iEEEKYJAkoIYQQJsnsLvEppRyAj4FUYJvW+ksjlySEEKIKmMQISim1TCl1RSl1rMD2AUqpE0qp00qp17M3DwPWaq1HA49Ue7FCCCGqhUkEFLACGJB3g1LKGlgMPAR0BEYopToCLYHY7JdlVGONQgghqpFJBJTWejtwvcDmHsBprXWU1joV+AoYApwnK6SghPqVUmOUUvuVUvuvXr1aFWULIYSoQqZ8D6oFf4+UICuYegILgY+UUgOBjcW9WWsdCoQC+Pr6VmjJ9vW/X2DOlhNcjE+meYM6TOvfgaFdW1RkV0IIIcrJlAOqSFrrJOD5srxWKTUYGNy2bdtyH2f97xeYse4of6WkgJU1F+KTmbHuKICElBBCVANTDqgLgGue5y2zt5WZ1nojsNHX13d0eQ8+Z8sJktMySDj0PTd+WUqtuk5Y13Nh9KYm7BnQHVdX13yPRo0aoZQq72GEEEIUw5QDah/QTinlQVYwPQU8XZ4dVGYEdTE+GYDaTdpSz3co6QlXybgVR/zZY8ydu520tLR8r69Tp06h0Cr4qFevXrnrEEKImsokAkoptQboAzgrpc4DM7XWS5VS44EtgDWwTGv9R3n2W5kRVPMGdeh2ayvTW39D8zZxXNTOfJg+nAP1+rJjeh/+/PNPYmNji3xs3bqVS5cukZmZmW+f9erVKzXE6tSpU95ShRDCIqma0PLd19dXl3ctvn0bPsXrwJvUUam525J1bY51m0X3R14u9f1paWlcunSp2BCLjY3lypUrhd7n5ORUYoC1aNGC2rVrl+tchBDClCmlDmitfQttt+SAynOJb/SpU6fK9+b5XnAzli2n0/nyaBr+rtb4u1rTqW0rrKeWayBXrJSUFC5cuJAvtGJiYvI9j4+PL3hONGnSpMQQa9asGdbW1iUeW2YoCiFMRY0MqBwVGUER1ADQLP89lRk/3ebPpKw/J8fa0Ovevvj7++Pv70/Pnj2pX7++4YvOlpiYWOIoLDY2lqSkpHzvsba2pnnz5sUG2LH4Wry/7RIp6X9fgqxjY83sYZ0lpIQQ1a5GBpQhRlAAWmui4zXhsRmEX7EjPKkVR44cITMzE6UUXl5euYHl7+9PmzZtqm1Gn9aa+Pj4EgPs/Pnz3L59O/8brWtRy9EZa0dnatVzwdrRGecmzVk4um9ukDVo0EBmJgohqlyNDKgcFRpBHfkGNk6EtOS/t9nUgcELwXs4CQkJ7N27l/DwcMLDw9m1axc3b94EwMXFJV9gdevWzaiTH7TWXL16NTewXlz8Pem3rpKeEEfGrbisGYoJ10Dnn9Th4OBQ6qSOunXrGumshBCWQgKqIg0Lj3wDP70DN89D/ZbwwFvgPbzIl2ZmZhIZGZkbWOHh4Zw8eRIAGxsbfHx88oVW8+bNK3NKlXLX+z9zIT453zadmYGLdTIfDWlV7Ejs8uXLFPx5adCgQYkB1rJlS+zs7Krz9IQQZkYCyggdda9evcru3btzA2vv3r2kpKQA0KpVq3yB5e3tTa1a1TPrP2eVjOS0v9faLcs9qNTUVC5evFji5cS4uLhC73NxcSkxxJo3b46NjU2VnKsQwvTVyICq1D2oKpCamsrhw4fZuXMn4eHh7Ny5k4sXLwJgb29Pz549cwOrV69eNGrUqMpqqapZfMnJyZw/f77EEMu5FJrDysqKpk2blhhiTZs2xcqqbGsbywxFIcxLjQyoHMYaQZVGa01sbGy+y4KHDh0iIyNrZOPp6ZlvlNWhQweLmLSQkJBQ6szEv/76K997atWqRYsWLUoMMWdnZ/576GKFRodCCOORgDLBgCpKUlIS+/btyxdaN27cAKBRo0b4+fnlBlb37t1xcHCo2IHKcX+tummtuXHjRqkzE1NTU/O9z87ODhyc0A5O1KrnjLWjC7XqOVPL0ZlmLVry69tPVOlHAoQQFSMBZSYBVVBmZiYnT57MF1iRkZFA1ued7rzzznyjLDc3t9J3WsoMRXOQmZmZb2ZizuPjTXvISIgj/VYcGYmFZyY6OjqWOjPR3t7eSGclRM1UIwPK1O5BGcr169fzTb7Ys2dP7iWxFi1a5AusO++8s/DSSHk+45VPfVeYfKwazqDq5J2hqDMzyEi8TvqtOOpn3mS0T/1Cgfbnn38W2kejRo1KXW7K1ta2uk9NCItVIwMqhzmPoMoiPT2dI0eO5BtlnTt3Dsi67NW9e/fcwPLz88NlcTugqO+7gqD4Irabj/LOULx9+3ah5aYKPq5fL9jsmTItN1VdszKFMHcSUBYcUEU5f/48u3btyg2sgwcPkp6eDkA7l9r4N9e56wt2dLHCSimLGEGB4WfxJSUllTozMSEhId97rKysSlxuytXVlcaNG5d5ZmJVnJcQpkICqoYFVEHJycns378/K7C+X0v43gPE/ZX1va9vC35utfHvOxT/YWPo0aMHjo6ORq7YvNy8ebPUmYk5n4HLYWNjQ8uWLUsMsZxGmBX97JoQ5kACqoYHVEH68Nec+upNdkVeIPxPO8Lj6vLHmfNorbGyssLb2zvfvSx3d3eLmOJuLFprrl27VurMxJxRbo6cRph/ZtQlvU6j7NmJf6+f6Orqyp6gwUY6KyEMo0YGlKVOkqgq8fHx7NmzJ/ey4O7du0lMTASgadOm+QLLx8dHJgoYWGZmZpGNMGNiYtiw8ygZCVfJSLxBwfuH0ghTmLsaGVA5ZARVMRkZGRw7dizf5IuoqCgAateuja+vb77QatKkiZErtlw5sxN1RnrWzMSEq2TcisMhLZ7BbW3zBdrVq1cLvV8aYQpTJgElAWUQly9fzjf5Yv/+/bkfmG3Tpk2+wOrUqVOpjRNF2ZTnHlRKSkqpkzqqqhGmEBUhASUBVSVu377NwYMHcwNr586duZ8tcnR0pFevXtXW3NHSGXIWn6EaYbq5ueV77uLiIvcqRblJQElAVQutNdHR0fkuC5pKc0dRdhVthGlra1vqzMSiGmHKFPqaTQJKAspoqr25owmvM2hJCjbCLOpx4cKF3MWPcxRshJlUqz6/XtBkOjSilqML1vWccXCoK1PoaxAJKAkok5GZmUlERES+UVbOLMtKN3e0gHUGLUlGRgaXL18uMcQuXb4MBf4dsrJ1wK5hY/r4eEojzBqgRgaUTDM3H1evXs03+WLfvn0Va+5owesMWir3aetJT7xO+q2rfy/0m3CV9Ftx3OGYKo0wa4AaGVA5ZARlflJTUzl06FC+yRc5zR0dHBwKNXds2LBh1huDGmCp6wxaqrwL/ObVokEddr5+P2AajTBF1ZGAkoAya6U1d+zYsWNWYCVuxt/pJu2drPLfiJcRlMky1DJOVdkIUybyVC0JKAkoi1NSc0enOgo/V2v8W1rj72FP99Eh2PccZeSKRXGqYxZfZRphljYzsX79+sWGmMxQLJ0ElASUxctt7vh1CDs3f8WuqFtExmU1LKxVq1ah5o6urq5GrliYmryNMGNiYooMsYsXL5KZmb8RZt26dYsMrpjbdVh1JJH0Oo2wqp01qUMW+S1MAkoCqkYqqbljy5YtCzV3lJvqojTp6elcunSpxJFYUY0wrezq5i70W8+5KRMf6VVoZmJNXd9SAkoCSlByc8c6deoUau7o7Oxs5IqFOcpphHnXm9+QlhBHxq2rpOf7bxyZKQmF3te4ceNSZyZaYiNMCSgJKFGMCxcuFGrumJaWBkD79u3zjbI8PT1l1pcos5JmKP4woWeFGmE2a9asxBBr0qSJ2f2MSkBJQIkySk5O5sCBA/lGWTkrhDdo0KDQ+oJ169Y1csXCVFV2hmJFG2GWNjPRycnJpGYmWkxAKaVaA28A9bXWj5flPRJQojK01pw5c4adO3fmBtYff/yR29yxS5cu+UZZrVq1Mqm//MK4qnIWX2UaYZZlZmJ1nZdJBJRSahkwCLiitfbKs30AEAJYA59rrd8vw77WSkAJYympuWOzZs3yBVbXrl1r7M1vYXzFNcLMt9zUpUuFZiY6OjoWG14nEmqzcM91bvP3pKLKzE40lYC6B0gEVuUElFLKGjgJ9AXOA/uAEWSF1ewCu3hBa30l+31VHlBhUWGEHAzhctJlmjo0ZZLPJAa2Hlju/QjLV1JzR1tb23zNHf38/CrW3FEWwRVVJC0trdSZiVeuXCn0Pqs69ahVz4Wmz81HKat8q3+Uh0kEVHYh7sCmPAHlBwRprftnP58BoLUuGE4F91NiQCmlxgBjANzc3LrlzNQqq7CoMILCg0i8kUjq1VTs29hjZ21HkH+QhJQoE4M2d5RFcIWRpaSkcOHCBWJjY3kieCPpCXGk37qKTkvBedBUABRw9v3y//toygH1ODBAa/1S9vNRQE+t9fhi3u8EvEvWiOvz0oIMKjaC6re2H5eSLnFl/RWurL+CfXt7nB9ypr1/e7YO31qufQkBlWzuKIvgChNSlvUTy6O4gDK7CfVa62vAK2V5bZ7VzMt9nMtJlwFw6u+Elb0V17ZcIyYkhstfXyY0PpRRo0ZVvm+RqFFsbW3x8/PDz8+PqVOnFtnc8V//+lfRzR3PnqNNQ1V48sXN88Y5GVGjTevfocjZidP6dzDocUxhBFWhS3zlUZkRVA6dobm1/xbxW+JJiErAxcWF8ePHM27cOPkwpzCYEps72iv8Xa1zH92aWVPH2U1GUMIoLG4WX3Yh7uQPqFpkTZJ4ALhA1iSJp7XWfxjgWBXuB5VzDyol4+/PGNhZ2zHTbyaO5x0JDg5m06ZN1KlTh4CAACZPnky7du0qW7IQ+WRmZhIZGUn4NwsJ37iK8JhUTl7Lmm1lYw0+Hdvi/+DgijV3FMJEmERAKaXWAH0AZ+BPYKbWeqlS6mFgAVkz95Zprd815HGrahZfZGQk8+bNY9WqVaSlpTF06FACAwPx9/c3ZPlCZMmexXf1Ygy7bzQkXHcl/PQN9u7dW7HmjkKYCJMIKGOp6s9BXb58mcWLF/Pxxx9z/fp1/Pz8CAwMZMiQISXPzBLCAFJTUzl8+HDuB4kLNnfs0aNHvinuuc0dhTARNTKgqrvle1JSEitWrGDevHlERUXRpk0bpkyZQkBAAPb29lV+fCGgHM0dsx/t27eXlS+EUdXIgMpR3StJZGRksH79eubMmcOePXtwcnJi3LhxvPrqqxX7gKYQlVRSc8dGjRrlC6zu3bvLL1SiWklAGWGpI601O3fuJDg4mA0bNlC7dm2effZZpkyZwh133FHt9QiRI7e5Y57AioyMBKS5o6h+NTKgqvsSX0lOnDjB/PnzWblyJSkpKQwePJjAwEDuvvtuubwiTII0dxTGUiMDKocpLRZ75coVPv74YxYvXkxcXBzdu3dn2rRpPProozLbSpiUam3uKOsM1mgSUCYSUDn++usvVq1axbx58zh16hQeHh5MnjyZ559/XvoLCZNVJc0dZZ3BGq9GBpQpXeIrTkZGBhs3biQ4OJidO3fSsGFDxo4dy/jx42nWrJmxyxOiRCU1d6xfvz5+fn65gdWjRw8cHR0L70TWGazxamRA5TDFEVRRdu3axdy5c1m3bh02NjaMHDmSKVOm0KlTJ2OXJkSZ5DR3zBtYx44dy23u6O3tnW+U5e7ujnq7IVDUv0MKguKr+xSEEUhAmUFA5Th9+jQLFixg2bJlJCcn8/DDDxMYGEifPn1kQoUwOzdv3mTPnj25HyTO29yxadOm+Lsk4d80FX9Xa3yaWWNbK/tnXEZQNYYElBkFVI5r167xySefsGjRIq5cuYKPjw+BgYE8/vjjMoNKmK1CzR23bSXqfFbbEVtr8G1ujX8rW/wfG4ffiED57GANUCMDyhzuQZVFSkoKq1evZu7cuZw4cQI3Nzdee+01XnrppaKv6QthZi7/HEr46lnsOvEn4Zes2X8hjdS0dKACzR2F2amRAZXDXEdQBWVmZhIWFkZwcDDbt2+nfv36vPzyy0ycOJEWLcq3zL20sxemrFLNHYXZkYCygIDKa+/evcydO5e1a9diZWXF008/zdSpU/H29i71vcW1EpF29sJUFdXc8ciRI0U3d/T3p02bNnK/1oxIQFlYQOU4e/YsCxYsYOnSpSQlJdG/f38CAwN54IEHiv0LWrAZY45mDs344fEfqrpkIQyixOaOLi75Aqtbt27SAduESUBZaEDluH79Op9++ikLFy7k8uXLdOnShcDAQJ588slCEyq8V3qji5jWq1Acee5IdZUshEFlZmYSERGRb5SVc+/ZxsYGHx+ffKElzR1NR40MKEuZJFEet2/f5t///jfBwcFERETQokULXnvtNUaPHp17nV5GUKKmuHr1ar6VL/bt2yfNHU1QjQyoHDVhBFVQZmYmW7ZsITg4mJ9//hlHR0fGjBnDpEmTOJJ2RO5BiRopNTWVQ4cO5Zt8kdPc0d7enp49e+YGVq9evWjUqFHlDyrrDJZKAqqGBVReBw8eZO7cuXz99dcopXjyySfp8WQPNt3eJLP4RI1WWnNHT0/PfKOsDh06lG/yhawzWCYSUDU4oHLExMQQEhJCaGgoiYmJPPDAAwQGBtK/f3+Z8SREttKaO+ZdX7B79+44ODgUvzNZZ7BMJKAkoHLFx8cTGhpKSEgIFy9exMvLi6lTpzJixAhsbW2NXZ4QJqWk5o7W1tZ07dq1+OaOQQ2QdQZLJwElAVVIamoqX331FcHBwRw9epRmzZoxceJEXn75ZRo2bGjs8oQwWdeuXcvX3HHv3r1FN3eMXsidda9iY13gCoWMoPKRgJKAKpbWmq1btxIcHMzWrVtxcHDgpZde4rXXXsPd3d3Y5Qlh8tLS0go1d4yJiQGgjg10b26Nf0tr/F2t8fOoi/NTH8k9qDxqZEDVxGnmlXX48GHmzp3LmjVryMzM5IknniAwMBBf30I/O0KIEpw/fz5rivuGFez89Rd+v5BMembW1yrc3NFC1ciAyiEjqPI7f/48Cxcu5NNPP+XWrVv06dOHwMBAHnrooRr9F0mIikpOTmb//v35RllxcXFAOZo7WigJKAmoCrl16xaff/45CxYsIDY2Fk9PT6ZOncozzzyDnZ2dscsTwmxprTl9+nS+wPrjjz9Kbu5oobNtJaAkoColLS2Nb7/9ljlz5nDo0CGaNGnChAkTeOWVV3BycjJ2eUJYhPj4ePbs2ZMbWIWaO+YJLB8fH4uZdSsBJQFlEFprfvnlF+bMmcP333+Pvb09L7zwApMnT6Z169bGLk8Ii1KouWN4OFFRUQDY2tri6+ubG1h+fn5m29xRAkoCyuCOHTvG3Llz+fLLL8nIyGDYsGEEBgbSs2dPY5cmhMW6fPlyvvUF9+/fT2pqKmC+zR0loCSgqszFixdZtGgRS5YsIT4+nt69exMYGMjgwYONMqFCmjGKmsRozR0NuMagBJQEVJVLSEhg2bJlzJ8/n3PnztG+fXumTp3KqFGjqq0XjzRjFDVdtTR3NPAagxJQElDVJj09ne+++445c+Zw4MABXFxcGD9+POPGjcPZ2blKjy2tRIQo7NatW4WaO966dQuoYHNHA68xaFEBpZQaCgwE6gFLtdYl/ssjAWUcWmu2b99OcHAwmzZtok6dOgQEBDB58mTatWtXJceUZoxClK7SzR0NvMagyQSUUmoZMAi4orX2yrN9ABACWAOfa63fL8O+GgLBWusXS3qdBJTxRUZGMm/ePFatWkVaWhpDhw4lMDAQf39/gx5HRlBCVEy5mjv+PIpaiecL78TcR1BKqXuARGBVTkAppayBk0Bf4DywDxhBVljNLrCLF7TWV7LfNxf4Umt9sKRjSkCZjsuXL7N48WI+/vhjrl+/jp+fH4GBgQwZMsQgs43kHpQQhlFSc0eHOrb0bKbxbwH+rtb0bVOLWrb2lnEPSinlDmzKE1B+QJDWun/28xkAWuuC4ZTzfgW8D2zVWv9Y2vEkoExPUlISK1asYN68eURFRdGmTRumTJlCQEAA9vb2ldq3zOITwvAKNXf8KYxDx6OwrwU33uuIdd+ZljGLr4iAehwYoLV+Kfv5KKCn1np8Me+fCDxH1kjrkNZ6SRGvGQOMAXBzc+t27ty5KjgTUVkZGRmsX7+eOXPmsGfPHpycnBg3bhyvvvqq2X7oUIiaIikpiVOnTnHnnXdWaj/FBZRZrvqptV6ote6mtX6lqHDKfk0o8DZwsHbt2tVboCgza2trHnvsMXbt2sVvv/3G3XffzaxZs2jVqhVjxozhxIkTxi5RCFEMBweHSodTSUwloC4AedpQ0jJ7W6VorTdqrccY7INposoopbjrrrv4z3/+w/HjxwkICGD16tXccccdDBkyhB07dmCOM06FEBVnKgG1D2inlPJQStUGngI2GLkmYSTt27dnyZIlnDt3jpkzZxIeHs4999xDr169+Pbbb0lPTzd2iUKIalDtAaWUWgPsAjoopc4rpV7UWqcD44EtQCTwjdb6DwMca7BSKvTmzZuV3ZUwgsaNGxMUFMS5c+f45JNPuHHjBsOHD6d9+/YsWrQod5VnIYRlMssP6paXzOKzDBkZGWzcuJHg4GB27txJw4YNGTt2LOPHj6dZs2bGLk8IUUEWNUmirGQEZVmsra0ZOnQov/32G+Hh4dx///3Mnj0bd3d3XnzxRSIiIoxdohDCgCw6oGSShOXy8/Nj7dq1nDx5ktGjR7NmzRo6derEoEGD2LZtm0yoEMICWHRACcvXtm1bPvroI2JjY/nXv/7Fvn37uO++++jevTtr1qyRCRVCmDGLDii5xFdzODk58eabb3Lu3Dk+++wzEhMTefrpp2nTpg0LFiwgISHB2CUKIcrJogNKLvHVPHZ2drz00ktERESwYcMG3N3dmTx5Mq6urrz++utcuFDpj9cJIaqJRQeUqLmsrKwYPHgwv/76K3v27KF///7MmTMHDw8PAgICOHr0qLFLFEKUotwBpZRyyF593OTJJT4B0KNHD77++mtOnz7N2LFjWbt2Ld7e3gwYMIAff/xRJlQIYaJKDSillJVS6mmlVJhS6gpwHLiklIpQSs1RSrWt+jIrRi7xibw8PDwICQkhJiaG9957j8OHD9O3b1+6du3KF198QVpamrFLFELkUZYR1C9AG2AG0FRr7aq1bgz0BnYDHyilRlZhjUIYVKNGjZgxYwbR0dEsW7aMtLQ0Ro0aRevWrQkODkZG3EKYhlJXklBK2WitS/zVsiyvMSZZSUKURGvN999/T3BwMD///DOOjo6MGTOGSZMm4erqWvoOhBCVUuGVJHKCRym1q8AOHZVSXfO+xtTIPShRFkopHnroIX766ScOHDjA4MGDWbBgAa1bt2bkyJH8/vvvxi5RiBqpPJMkbAGUUvMAtNYJwMdVUZShyD0oUV4+Pj58+eWXREVFMWHCBP773//i4+PDgw8+yJYtW2RChRDVqDwBpZRSTYCR2S3XAepUQU1CGJ2bmxvz5s0jNjaWDz74gMjISAYMGIC3tzcrV64kNTXV2CUKYfHKE1AzgB3Av4H5Sqlx5Xy/EGanQYMGTJ8+nbNnz7Jy5UqUUgQEBODh4cEHH3xAfHy8sUsUwmKVOWC01t9rrdtrrV8DvgbaAi9WWWVCmJDatWvz7LPPcvjwYbZs2UKnTp14/fXXcXV1ZfLkyURHRxu7RCEsTllm8SldyovK8hpjUEoNBga3bdt29KlTp4xdjrAwhw8fZu7cuaxZswatNU888QRTp07F17fQZCQhRAkq0w/qF6XUBKWUW4Ed1lZK3a+UWgk8Z6hCDUkmSYiq1KVLF1atWsXZs2eZMmUKmzdvpnv37tx3332EhYWRmZlp7BKFMGtlCagBQAawRil1MXsFiSjgFDACWKC1XlGFNQph0lq2bMmHH35IbGwsc+fO5cyZMwwaNAgvLy+WLl1KSkqKsUsUwiyVq+W7UsoGcAaStdZmc3dYPqgrqlNaWhrffPMNwcHBHDp0iCZNmjBhwgReeeUVnJycKrzfsKgwQg6GcDnpMk0dmjLJZxIDWw80YOVCGIdBWr5rrdO01pfMKZyEqG42NjY888wzHDx4kB9//JGuXbvy5ptv4ubmxoQJE4iKiir3PsOiwggKD+JS0iU0mktJlwgKDyIsKqwKzkAI01DmEZRSakoRm28CB7TWhwxalYHJCEoY27Fjx5g7dy5ffvklGRkZDBs2jMDAQHr27Fmm9/db249LSZcKbW/m0IwfHv/B0OVWOxkd1myGGEH5Aq8ALbIfL5N1f+ozpdR0g1QphIXy8vJi+fLlREdHM336dLZu3UqvXr2455572LBhQ6kTKi4nXS7XdnMio0NRnPIEVEvAR2s9VWs9FegGNAbuAQKqoLZKk7X4hKlp3rw5s2fPJjY2lgULFhATE8OQIUPw9PQkNDSU5OTkIt/X1KFpubabk5CDIaRk5J9IkpKRQsjBECNVJExFeQKqMXA7z/M0oInWOrnAdpMh08yFqXJ0dGTSpEmcPn2ar776CkdHR15++WVatWrFO++8Q1xcXL7XT/KZhJ21Xb5tdtZ2TPKZVJ1lVwlLHh2KyilPQH0J7FFKzVRKBQE7gX8rpRyAiKooTghLV6tWLZ588kn27dvHtm3b6NmzJzNnzsTV1ZVx48aR8wHzga0HEuQfRDOHZigUzRyaEeQfZBH3aSx5dCgqp7zTzH2Bu7Kf7tRam8XMA5kkIcxJREQE8+bNY/Xq1aSlpTF06FACAwPx9/c3dmlVIuceVN7LfHbWdhYTwKJ0BplmTtZlvd75yWkAACAASURBVEyyPrhrkj2ghDB3HTt25PPPP+fcuXP885//ZNu2bdx11134+/uzbt06MjIyjF2iQVny6FBUTnmmmU8CRgPfAQp4FAjVWi+quvIMQ0ZQwpwlJSWxfPly5s2bx9mzZ2nTpg1TpkwhICAAe3t7Y5cnRKUVN4IqT0AdAfy01knZzx2AXVprb4NWWgUkoIQlyMjIYN26dQQHB7N3716cnJwYN24cr776Kk2aNDF2eUJUmCEu8SmyLu3lyMjeJoSoBtbW1jzxxBPs3r2bHTt20Lt3b2bNmkWrVq0YM2YMx48fN3aJQhhUeQJqOVmz+IKUUm8De4BlVVOWEKI4Sil69+7N+vXriYyMJCAggFWrVuHp6ckjjzzC9u3bpTW9sAjlaVg4D3geuAbEAc9predXVWHFUUp5KqWWKKXWKqXGVvfxhTAlHTp0YMmSJcTExPDWW28RHh7OvffeS8+ePfn2229JT083dolCVFipAaWUSlBK3VJK3QK2Ae9lP3ZkbyszpdQypdQVpdSxAtsHKKVOKKVOK6VeL2kfWutIrfUrwHD+nvIuRI3WuHFj3n77bWJiYvj444+5ceMGw4cPp3379ixatIjExERjlyhEuZUaUFprR611vTwPxzyPeuU83gqy1u/LpZSyBhYDDwEdgRFKqY5Kqc5KqU0FHo2z3/MIEAZsLufxhbBo9vb2jB07luPHj7Nu3TqaNWvGxIkTcXNz44033uDSpcILzgphqsr7OahK0VpvB64X2NwDOK21jtJapwJfAUO01ke11oMKPK5k72eD1voh4JnqrF8Ic2Ftbc2jjz7Kzp072blzJ/fddx+zZ8/G3d2dF198kT/++MPYJQpRqmoNqGK0AGLzPD+fva1ISqk+SqmFSqlPKWEEpZQao5Tar5Taf/XqVcNVK4SZ8ff357vvvuPkyZO89NJLrFmzBi8vLwYOHMgvv/wiEyqEyTKFgCoXrfU2rfVErfXLWuvFJbwuVGvtq7X2dXFxqc4ShTBJbdu2ZfHixcTExPDOO++wb98+7r//fnx9fVmzZg1pabI4jDAtphBQFwDXPM9bZm+rNGm3IURhzs7O/N///R/nzp0jNDSUpKQknn76adq2bcv8+fNJSEgwdolCAKYRUPuAdkopD6VUbeApYIMhdiztNoQoXp06dRg9ejQRERFs2LABd3d3pkyZgqurK//4xz+4cMEgvycKUWHVGlBKqTXALqCDUuq8UupFrXU6MB7YAkQC32itDXIHV0ZQQpTOysqKwYMH8+uvv7Jnzx769+9PcHAwHh4eBAQEcPToUWOXKGqocrXbMFeyFp8Q5XP27FkWLFjA559/zl9//UX//v0JDAzkgQceQClZ4UwYlqHabZgVGUEJUTEeHh6EhIQQGxvLe++9x+HDh+nbty9du3bliy++kAkVolpYdEDJPSghKqdRo0bMmDGD6Oholi5dSmpqKqNGjcLDw4Pg4GDklz9RlSw6oIQQhmFra8sLL7zAsWPHCAsLo3379kybNg1XV1cCAwOJjY0tfSdClJNFB5Rc4hPCsKysrHj44Yf5+eef2b9/P4MGDWLBggW0bt2akSNH8vvvvxu7RGFBLDqg5BKfEFWnW7du/Pvf/+bMmTNMmDCB//73v/j4+PDggw/y/fffywoVotIsOqCEEFWvVatWzJs3j9jYWD744AMiIyN56KGH8Pb2ZsWKFdy+fdvYJQozZdEBJZf4hKg+DRo0YPr06Zw9e5aVK1eilOL555/Hw8OD999/nxs3bhi7RGFmLDqg5BKfENWvdu3aPPvssxw+fJgtW7bg5eXFjBkzcHV15bXXXiM6OtrYJQozYdEBJYQwHqUU/fr144cffuDQoUM8+uijLF68mDZt2vDUU08hH54XpZGAEkJUuS5durB69WqioqKYMmUKmzdvpnv37tx3332EhYWRmZlp7BKrTFhUGP3W9sN7pTf91vYjLCrM2CWZDYsOKLkHJYRpcXV1Zc6cOcTGxhIcHMzp06cZNGgQXl5eLF26lJSUFGOXaFBhUWEEhQdxKekSGs2lpEsEhQdJSJWRRQeU3IMSwjTVr1+fqVOnEhUVxRdffIGtrS0vvfQS7u7uvPvuu1y7ds3YJRpEyMEQUjLyh25KRgohB0OMVJF5seiAEkKYNhsbG5555hkOHjzI1q1b6dq1K2+++SZubm5MmDCBqKgoY5dYKZeTLpdru8hPAkoIYXRKKR588EH+97//ceTIEYYPH86nn35Ku3bteOKJJ9izZ4+xS6yQpg5Ny7Vd5CcBJYQwKZ07d2b58uVER0czffp0tm7dSq9evbj77rv573//a1YTKib5TMLO2i7fNjtrOyb5TDJSRebFogNKJkkIYb6aN2/O7NmziY2NZcGCBcTGxjJ06FA8PT0JDQ0lOTnZ2CWWamDrgQT5B9HMoRkKRTOHZgT5BzGw9UBjl2YWpGGhEMIspKen89133zFnzhwOHDiAi4sL48ePZ9y4cTg7Oxu7PFEJNbJhoRDCctSqVYsnn3ySffv2sW3bNnr27MnMmTNxdXVl3LhxnDp1ytglCgOTgBJCmBWlFPfeey8bN24kIiKCkSNHsnTpUjp06MCwYcMIDw83donCQCSghBBmy9PTk88++4xz587xz3/+k23btnHXXXfh7+/PunXryMjIMHaJohIkoIQQZq9p06bMmjWL2NhYFi1axOXLl3nsscfo0KEDH3/8MX/99ZexSxQVIAElhLAYDg4OjB8/nlOnTvHtt9/i5OTEq6++ipubG2+99RZ//vmnsUsU5WDRASXTzIWomaytrXn88cfZvXs3O3bsoHfv3syaNYtWrVoxZswYTpw4YewSRRlYdEDJWnxC1GxKKXr37s369euJjIwkICCA1atXc8cddzBkyBB27NghrelNmEUHlBBC5OjQoQNLlizh3LlzzJw5k507d3LPPffQq1cvvv32W9LT041doihAAkoIUaM0btyYoKAgYmJi+OSTT7h+/TrDhw+nffv2LFq0iMTERGOXKLJJQAkhaiR7e3teeeUVjh8/zn/+8x+aN2/OxIkTcXNz44033uDSpUvGLrHGk4ASQtRo1tbWDB06lN9++43w8HDuv/9+Zs+ejbu7Oy+++CIRERHGLrHGkoASQohsfn5+rF27lpMnTzJ69GjWrFlDp06dGDRoENu2bZMJFdVMAkoIIQpo27YtH330ETExMbzzzjvs3buX++67j+7du/PVV1/JhIpqIgElhBDFcHZ25v/+7/+IiYkhNDSUxMRERowYQZs2bViwYAEJCQnGLtGiSUAJIUQp7OzsGD16NBEREWzYsAF3d3cmT56Mq6srr7/+OhcuXDB2iRbJLANKKeWglNqvlBpk7FqEEDWHlZUVgwcP5tdff2XPnj3079+fOXPm4OHhQUBAAEePHjV2iRalWgNKKbVMKXVFKXWswPYBSqkTSqnTSqnXy7CrfwDfVE2VQghRuh49evD1119z+vRpxo4dy9q1a/H29mbAgAH8+OOPMqHCAKp7BLUCGJB3g1LKGlgMPAR0BEYopToqpTorpTYVeDRWSvUFIoAr1Vy7EEIU4uHhQUhICDExMbz33nscPnyYvn370rVrV7744gvS0tKMXaLZqtaA0lpvB64X2NwDOK21jtJapwJfAUO01ke11oMKPK4AfYBewNPAaKVUkeeglBqTfRlw/9WrV6vupIQQAmjUqBEzZswgOjqaZcuWkZaWxqhRo2jdujXBwcHIotXlZwr3oFoAsXmen8/eViSt9Rta69eAfwOfaa0zi3ldqNbaV2vt6+LiYtCChRCiOLa2tjz//PMcO3aMzZs30759e6ZNm4arqyvTpk0jNja29J0IwDQCqkK01iu01ptKeo202xBCGItSioceeoiffvqJAwcOMHjwYObPn0/r1q0ZNWoUhw4dMnaJJk9V9408pZQ7sElr7ZX93A8I0lr3z34+A0BrPdtQx/T19dX79+/Pty0tLY3z58+TkpJiqMMIM2dnZ0fLli2xsbExdinCQsXExBASEpL7maoHH3yQwMBA+vXrh1LK2OUZjVLqgNbat9B2EwioWsBJ4AHgArAPeFpr/YcBjjUYGNy2bdvRp06dyve1s2fP4ujoiJOTU43+wRBZtNZcu3aNhIQEPDw8jF2OsHDx8fGEhoYSEhLCxYsX8fLyIjAwkBEjRlC7dm1jl1ftiguo6p5mvgbYBXRQSp1XSr2otU4HxgNbgEjgG0OEE5TcsDAlJUXCSeRSSuHk5CQjalEtGjRowPTp0zl79iwrV65EKUVAQAAeHh588MEHxMfHG7tEk1Dds/hGaK2baa1ttNYttdZLs7dv1lq311q30Vq/a6jjlXYPSsJJ5CU/D6K61a5dm2effZbDhw+zZcsWOnXqxOuvv46rqyuTJ0/m3Llzxi7RqMx2kkRZSMt3IYQ5UErRr18/fvjhBw4dOsSjjz7KRx99RJs2bRgxYgQHDhwwdolGYdEBZeqUUowcOTL3eXp6Oi4uLgwa9PcKTv/73//w9fWlY8eOdO3alalTpxqjVCFENenSpQurVq3i7NmzTJkyhc2bN+Pr68t9991HWFgYmZlFfrLGIll0QJn6NHMHBweOHTtGcnIyAFu3bqVFi78/Anbs2DHGjx/PF198QUREBPv376dt27bGKlcIUY1atmzJhx9+SGxsLHPnzuXMmTMMGjQILy8vli5dWiPul1p0QBnyEt/63y9w1/s/4/F6GHe9/zPrfzfM6sUPP/wwYWFhAKxZs4YRI0bkfu3DDz/kjTfe4I477gCyOn+OHTvWIMcVQpiHevXqMWXKFM6cOcOXX36Jra0tL730Eu7u7rz77rtcv15wcR7LYdEBZSjrf7/AjHVHuRCfjAYuxCczY91Rg4TUU089xVdffUVKSgpHjhyhZ8+euV87duwY3bp1q/QxhBDmz8bGhqeffpqDBw/y448/0rVrV958801cXV2ZOHEiUVFRxi7R4Cw6oAx1iW/OlhMkp2Xk25aclsGcLScqtV8Ab29voqOjWbNmDQ8//HCl9yeEsGxKKR544AH+97//cfToUYYPH86SJUto164dw4cPZ8+ePcYu0WAsOqAMdYnvYnxyubaX1yOPPJL7Ib28OnXqVGNn7wghSufl5cXy5cuJjo5m+vTp/PDDD/Tq1Yt77rmHDRs2mP2ECosOKENp3qBOubaX1wsvvMDMmTPp3Llzvu3Tpk3jvffe4+TJkwBkZmayZMkSgxxTCGE5mjdvzuzZs4mNjWXBggXExMQwZMgQPD09CQ0NzZ2IZW4koMpgWv8O1LGxzretjo010/p3MMj+W7ZsycSJEwtt9/b2ZsGCBYwYMQJPT0+8vLws8jqzEMIwHB0dmTRpEqdPn+arr77C0dGRl19+mVatWvHOO+8QFxdn7BLLpdrX4qtOJa3FFxkZiaenZ5n3tf73C8zZcoKL8ck0b1CHaf07MLRrsV1BhJkq78+FEKZMa8327dsJDg5m06ZN1KlTh4CAACZPnky7du2MXV4uk1ks1hiKWs1c/iESRZGfC2GpIiMjmTdvHqtWrSItLY2hQ4cSGBiIv7+/sUszjcVihRBCGIenpyefffYZ586d44033uDXX3/lrrvuwt/fn3Xr1pGRkVH6TqqZBJQQQtQgTZs25V//+hcxMTEsWrSIy5cv89hjj9GhQwc+/vhj/vrrL2OXmEsCSgghaiAHBwfGjx/PqVOn+Pbbb3FycuLVV1/Fzc2Nt956iz///NPYJVp2QJn6WnxCCGFs1tbWPP744+zevZsdO3bQu3dvZs2aRatWrRgzZgzHjx83Wm0WHVDSbkMIIcpGKUXv3r1Zv349kZGRBAQEsGrVKjw9PXnkkUfYvn071T2pzqIDytS9++67dOrUCW9vb+68885Slyhxd3cv1+cYgoKCCA4OBiAgIIBt27ZVplyjW7FiBePHjy/ya3Xr1i3XvvL+2Qgh8uvQoQNLliwhJiaGmTNnEh4ezr333kvPnj355ptvSE9Pr5Y6JKCMZNeuXWzatImDBw9y5MgRfvzxR1xdXY1dVpWprh9oIYThNG7cmKCgIGJiYvjkk0+4ceMGTz75JO3atWPhwoUkJiZW6fEloMrqyDcw3wuCGmT998g3ldrdpUuXcHZ2xtbWFgBnZ2eaN28OwE8//UTXrl3p3LkzL7zwArdv385934cffkjnzp3p0aMHp0+fBiA6Opr7778fb29vHnjgAWJiYgodr379+tSuXbvQ9j59+pDzGbG4uDjc3d2BrNHKsGHDGDBgAO3atWP69OkAZGRkEBAQgJeXF507d2b+/PkAnD59mgcffJAuXbrg4+PDmTNn2LZtG3fffTePPPIIHTt2BGDo0KF069aNTp06ERoamltH3bp1mTZtGp06deLBBx9k79699OnTh9atW7Nhw4bc18XGxtKnTx/atWvH22+/XeSf7Zw5c+jevTve3t7MnDkzd/u7775L+/bt6d27NydOVH6hXyFqCnt7e1555RWOHz/Of/7zH5o3b86kSZNwc3Pjn//8J0lJSVVzYK21xT+6deumC4qIiCi0rViHv9Z6VhOtZ9b7+zGrSdb2CkpISNBdunTR7dq102PHjtXbtm3TWmudnJysW7ZsqU+cOKG11nrUqFF6/vz5WmutW7VqpWfNmqW11nrlypV64MCBWmutBw0apFesWKG11nrp0qV6yJAhWmutZ86cqefMmVNiHffee6/et2+f1lrrq1ev6latWmmttV6+fLn28PDQ8fHxOjk5Wbu5uemYmBi9f/9+/eCDD+a+/8aNG1prrXv06KHXrVuXew5JSUn6l19+0fb29joqKir39deuXdNaa/3XX3/pTp066bi4OK211oDevHmz1lrroUOH6r59++rU1FR96NAh3aVLl9yamjZtquPi4nLfn1O7g4OD1lrrLVu26NGjR+vMzEydkZGhBw4cqH/99Ve9f/9+7eXlpZOSkvTNmzd1mzZtivyzKdfPhRA12M6dO/WwYcN0mzZtdFpaWqX2BezXRfzbLSOosvjpHUgrsNhiWnLW9gqqW7cuBw4cIDQ0FBcXF5588klWrFjBiRMn8PDwoH379gA899xzbN++Pfd9OSuejxgxgl27dgFZlwuffvppAEaNGsVvv/1W4bryeuCBB6hfvz52dnZ07NiRc+fO0bp1a6KiopgwYQLff/899erVIyEhgQsXLvDoo48CYGdnh729PQA9evTAw8Mjd58LFy6kS5cu9OrVi9jYWHKWoKpduzYDBgwAoHPnztx7773Y2NjQuXNnoqOjc9/ft29fnJycqFOnDsOGDSt0rj/88AM//PADXbt2xcfHh+PHj3Pq1Cl27NjBo48+ir29PfXq1eORRx4xyJ+REDWVv78/3333HUeOHKFWrVpVcgwJqLK4eb5828vI2tqaPn368Pbbb/PRRx/x3XfflfoepVSR/19RtWrVyl2Sv2AL6ZzLjzm1pqen07BhQw4fPkyfPn1YsmQJL730Uon7d3BwyP3/bdu28eOPP7Jr1y4OHz5M165dc49pY2OTez5WVla5x7ayssp3/6rgORd8rrVmxowZHDp0iEOHDnH69GlefPHFMv1ZCCHKL+eX0apg0QFlsM9B1W9Zvu1lcOLECfIuYHvo0CFatWpFhw4diI6Ozr2/tHr1au69997c13399de5//Xz8wOyfpP56quvAPjyyy+5++67y1yHu7t7bs+ptWvXlvr6uLg4MjMzeeyxx5g1axYHDx7E0dGRli1bsn79egBu375d5KfRb968ScOGDbG3t+f48ePs3r27zHXm2Lp1K9evXyc5OZn169dz11135ft6//79WbZsWe7N2wsXLnDlyhXuuece1q9fT3JyMgkJCWzcuLHcxxZCVK+qGZeZCK31RmCjr6/v6Ert6IG3YOPE/Jf5bOpkba+gxMREJkyYQHx8PLVq1aJt27aEhoZiZ2fH8uXLeeKJJ0hPT6d79+688sorue+7ceMG3t7e2NrasmbNGgAWLVrE888/z5w5c3BxcWH58uVlriMwMJDhw4cTGhrKwIEDS339hQsXeP7553NHXbNnzwaygvTll1/mrbfewsbGhm+//bbQewcMGMCSJUvw9PSkQ4cO9OrVq8x15ujRowePPfYY58+fZ+TIkfj65l9fsl+/fkRGRuaGd926dfniiy/w8fHhySefpEuXLjRu3Jju3buX+9hCiOolq5mX1ZFvsu453TyfNXJ64C3wHm7gSoWxyWrmQlS/4lYzt+gRlEF5D5dAEkKIamTR96CEEEKYLwkoIYQQJkkCSgghhEmSgBJCCGGSJKCEEEKYJLMLKKVUH6XUDqXUEqVUH2PXUxk1ud1GdHQ0Xl5eRX4t7wK2ZbFt2zYGDRpkqNKEECaiWqeZK6WWAYOAK1prrzzbBwAhgDXwudb6/RJ2o4FEwA6o3FpDRpS33YatrS1xcXGkpqYau6wyy8jIwNra2thlCCEsWHWPoFYAA/JuUEpZA4uBh4COwAilVEelVGel1KYCj8bADq31Q8A/gKL7LVSBsKgw+q3th/dKb/qt7UdYVFil9mcq7TYWLlxIx44d8fb25qmnngKyVrl4/vnn6dy5M97e3rlrBNatW5epU6fSpUsXdu3axTvvvEP37t3x8vJizJgxud02+/Tpw+TJk/H19cXT05N9+/YxbNgw2rVrx5tvvpl77PT0dJ555hk8PT15/PHHi1we6YcffsDPzw8fHx+eeOKJ3CWMvv/+e+644w58fHxYt25d+b8BQgjTV9QS51X5ANyBY3me+wFb8jyfAcwow35qA2tL+PoYYD+w383NrdDy7uVpq7DpzCbtu9pXe63wyn34rvbVm85sKvM+CjKVdhvNmjXTKSkpWuu/W2dMnz5dT5o0Kfc1169f11pntcT4+uu/W4zktM7QWuuRI0fqDRs2aK2zWnhMnz5da631ggULdLNmzfTFixd1SkqKbtGihY6Li9Nnz57VgP7tt9+01lo///zzubXmtAC5evWqvvvuu3ViYqLWWuv3339fv/3227l/RidPntSZmZn6iSeeyP2zqCxptyFE9cOE2220AGLzPD+fva1ISqlhSqlPgdXAR8W9TmsdqrX21Vr7uri4VKrAkIMhpGTkX+k7JSOFkIMhFd6nqbTb8Pb25plnnuGLL77IXTL/xx9/5NVXX819TcOGDYGsFc0fe+yx3O2//PILPXv2pHPnzvz888/88ccfuV/LaWfRuXNnOnXqRLNmzbC1taV169bExmZ9u11dXXMXex05cmShunfv3k1ERAR33XUXd955JytXruTcuXMcP34cDw8P2rVrh1KKkSNHlvl8hRDmw+yWOtJarwPKdE1HKTUYGNy2bdtKHfNy0uVybS+rnHYbffr0oXPnzqxcuZKuXbuW+B5Dt9sICwtj+/btbNy4kXfffZejR48W+1o7O7vc+04pKSmMGzeO/fv34+rqSlBQUL52HXnbZeRt25G3fUZZWmf07ds3d1HcHIcOHarAmQohDCksKoyQgyFcTrpMU4emTPKZxMDWpS84XR6mMIK6ALjmed4ye1ulaa03aq3H1K9fv1L7aerQtFzby8IU2m1kZmYSGxvLfffdxwcffMDNmzdJTEykb9++LF68OPd1N27cKPTenDBydnYmMTGxTK06CoqJickdBf773/+md+/e+b7eq1cvdu7cmftnkZSUxMmTJ7njjjuIjo7mzJkzAIUCTAhRtcKiwggKD+JS0iU0mktJlwgKD6r0vfmCTCGg9gHtlFIeSqnawFPABkPs2FD9oCb5TMLO2i7fNjtrOyb5TKrwPhMTE3nuuedyJyhEREQQFBSUr91G586dsbKyKrLdRkhICPPnzwey2m0sX74cb29vVq9eTUhI2S49ZmRkMHLkSDp37kzXrl2ZOHEiDRo04M033+TGjRt4eXnRpUsXfvnll0LvbdCgAaNHj8bLy4v+/ftXqH1Fhw4dWLx4MZ6enty4cYOxY8fm+7qLiwsrVqxgxIgReHt74+fnx/Hjx7Gzs8ttD+Lj40Pjxo3LfWwhRMVVxW2PolRruw2l1BqgD+AM/AnM1FovVUo9DCwga5r5Mq31u4Y8riHabVTHcFYYn7TbEKJ03iu90RTODoXiyHNHyr0/k2i3obUeUcz2zcBmQx/PUPegAAa2HiiBJIQQZN3euJR0qcjthmQKl/iqjKHuQQkhhPhbVdz2KIrZzeIzJK21QWbCCctQnZe7hTBnOVeTqvq2h0UHVEmX+Ozs7Lh27RpOTk4SUgKtNdeuXcPOzq70FwshquW2R7VOkjCWoiZJpKWlcf78+Xyf3RE1m52dHS1btsTGxsbYpQhRo5jEJAlTYmNjg4eHh7HLEEIIUQyLniRhqM9BCSGEqH4WHVAyi08IIcyXRQeUEEII81UjJkkopa4C5yqxC2eg7K1si1cfMOb1xqKOb0rnVpl9lPfcynqssryuur+vhvqeGYKhz72851aVf/aG2nfOfqr63MrzekP+XBvi57GV1rpw24mienDIo1BvqSJ7lVRgP6FGPo9Cxzelc6vMPsp7bmU9VlleV93fV0N9z0zl+16Zc6vKP3tD7TtnP1V9buV5vSF/rqvy51Eu8VWvjRZ8fEPsuzL7KO97y/r6srzO2N9XYzL2uZv6z3Rl9lNVP9Nlfa2xv7c14xJfZSml9usi5uhbAjk382Op5wVybuaoKs9LRlBlE2rsAqqQnJv5sdTzAjk3c1Rl5yUjKCGEECZJRlBCCCFMkgSUEEIIkyQBJYQQwiRJQAkhhDBJElDloJRyUErtV0oNMnYthqSU6qOU2qGUWqKU6mPsegxFKWWllHpXKbVIKfWcsesxJKXU3dnfr8+VUuHGrseQlFJuSqn1SqllSqnXjV2PoSilOiqlvlFKfaKUetzY9RiCUqq1UmqpUmpt9nMHpdRKpdRnSqlnKrt/CagClFKu06lspQAAA9FJREFUSqlflFIRSqk/lFJ5exj/A/jGWLVVVgnnpoFEwA44b7wKK6aE8xoCtATSMMPzguLPTWu9Q2v9CrAJWGncKiumhO9bZ2Ct1voFoKsRS6yQEs7rIWCR1nos8KwRSyy3En4Oo7TWL+Z56TCyvnejgUcqe9wa2w+qBOnAVK31QaWUI3BAKbUVaAFEkPWPuLkq7tx2aK1/VUo1AeYBlf7Np5oVd14dgHCt9afZv+H9ZNQqK6bIc9NaR2R//WngxeLfbtKK+77tBtYqpV4AVhu1woop7rxWAzOVUo8ATkatsPxK+znM0RI4mv3/GZU9qARUAVrrS8Cl7P9PUEpFkhVOfQAHoCOQrJTarLXONFqhFVDcueX5IbsB2Bqrvooq4Xt2HkjNflml/7IYQwnnFqGUcgNuaq0TjFljRZVwbg8DM7XW27N/sVhuxDLLrZS/Z68qpayBdcassbxK+jks8NLzZIXUIQxwhU4CqgRKKXeyLjHs0Vpvzd4WAMSZWzgVlPfclFLDgP5AA+AjI5ZVaXnPi6zf+hYppe4GthuxLIMocG6QNXIyq3+8i1Pg3C4BQUqpp4Fo41VVeQX+nrkD/yTrF905RiuqkgqckxPwLtBVKTUDWAh8pJQaiAHW8pOVJIqhlKoL/Aq8q7U2q992SmOp52ap5wVybubIEs+rus9JJkkUQSllA3wHfGkpP1g5LPXcLPW8QM7NHFnieRnjnGQEVYBSSpE1K+q61vo1Y9djSJZ6bpZ6XiDnZo4s8byMdU4SUAUopXoDO8iaiZJzn+mfWuvNxqvKMCz13Cz1vEDOzRxZ4nkZ65wkoIQQQpgkuQclhBDCJElACSGEMEkSUEIIIUySBJQQQgiTJAElhBDCJElACSGEMEkSUEIIIUySBJQQQgiTJAElhIlTSr2slNJKKc882yKVUh7GrEuIqiYBJYTp60xWf52BAEopO6AJZt6KQojSSEAJYfq8gQ/IDiiymmYe17JOmbBwElBCmL6OwH+Bxkqp+mSNqI4YtyQhqp4ElBAmTCnlClzTWicDW8nqfOxN1qrSQlg0CSghTFtn/g6jzWRd5usMHFFKuSqlPlNKBSulHjRahUJUkVrGLkAIUaK8o6VfgU+BOtnbugOpwEKtdYxxyhOi6sgISgjTljuC0lrfJuveU6rWOl5rvRVYBHyklGphxBqFqBLSsFAIM6WU+gCwBuyAyVrrNCOXJIRBSUAJIYQwSXKJTwghhEmSgBJCCGGSJKCEEEKYJAkoIYQQJkkCSgghhEmSgBJCCGGSJKCEEEKYJAkoIYQQJun/AdwTsSBhyF2cAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_0SHgi4g7LGx"
      },
      "source": [
        "## Random number generators and seed\n",
        "\n",
        "I have a final note on RNG with NumPy. You may have seen or used lot of things like `np.random.seed(123)` or `np.random.RandomState(123)`. Please don't.\n",
        "\n",
        "There are 3 problems:\n",
        "\n",
        "* Using a fix seed is in general dangerous in a sense that you can forget it and then instead of having a RNG, well it's not anymore random and hence just a MC draw. Also small values, such as commonly used 0, 123, etc. tend to produce bad entropy for the generator.\n",
        "* Using `np.random.seed(...)` fix the global seed. But new code are not relying on this at all. Hence you might only fix a portion of the code.\n",
        "* `np.random.RandomState` should not be used for new code because it's slower and has statistical issues. New code should use `np.random.Generator` which is better in every way.\n",
        "\n",
        "Have a look at the documentation for more details:\n",
        "https://scipy.github.io/devdocs/tutorial/stats.html#random-number-generation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OwCwpG1Rhhg-"
      },
      "source": [
        "# PyTorch\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VoGIgD4IfDru"
      },
      "source": [
        "### Basics\n",
        "\n",
        "The implementation of Sobol' sequences in PyTorch is essentially a C++ (PyTorch ATen) port of the `scipy` version, so their behavior is equivalent. There are some slight differences in the API (largely argument names), but otherwise using the PyTorch [`SobolEngine`](https://pytorch.org/docs/stable/generated/torch.quasirandom.SobolEngine.html) is very similar to using `scipy`'s `Sobol`. \n",
        "\n",
        "*Note:* \n",
        "Currently `SobolEngine` is limited to being run on the CPU (of course the generated sample tensors can easily be moved to GPU). [`cuRAND`](https://docs.nvidia.com/cuda/curand/host-api-overview.html#generator-types) supports generating Sobol' sequences on the GPU, but this is currently not exposed in PyTorch (always happy to accept PRs...)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uFDYUHTJ7N_o",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "outputId": "385b6ffc-767e-4552-cd57-c5efaf5f4693"
      },
      "source": [
        "import torch\n",
        "\n",
        "torch.__version__  # make sure we're good (1.9+)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            },
            "text/plain": [
              "'1.9.0+cu102'"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 53
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bMURUs6s1EHy"
      },
      "source": [
        "t_sobol_engine = torch.quasirandom.SobolEngine(\n",
        "    dimension=2, scramble=True, seed=29211402,\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "usxanizZ1Hr4",
        "outputId": "78c31ea0-583e-4528-f3c4-cc5bdd30d7eb"
      },
      "source": [
        "t_sobol_engine.draw(4)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "tensor([[0.9967, 0.5888],\n",
              "        [0.0412, 0.4169],\n",
              "        [0.2833, 0.9479],\n",
              "        [0.7389, 0.0573]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 55
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7gRvyyg1xpTB",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "4f387b93-703c-4009-b8a8-a7553e0bcc23"
      },
      "source": [
        "t_sobol_engine.reset()\n",
        "t_sobol_engine.draw_base2(2)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "tensor([[0.9967, 0.5888],\n",
              "        [0.0412, 0.4169],\n",
              "        [0.2833, 0.9479],\n",
              "        [0.7389, 0.0573]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 56
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_CLwohIlYqnO"
      },
      "source": [
        "### Illustrative example: Auto-differentiating through quasi-random samples.\n",
        "\n",
        "We consider the case where we draw random samples $\\xi$ from some distributon $\\mathcal{D}$, where $\\mathcal{D} = \\mathcal{D}(x)$ depends on some parameter $x$. We'd like to compute gradients of $\\xi$ w.r.t. $x$. We don't want to do that by hand. In various cases we can do this by combining the \"reparameterization trick\" [1] with auto-differentiation capabilities of modern Deep Learning frameworks. \n",
        "\n",
        "[1] *D. P. Kingma and M. Welling. Auto-Encoding Variational Bayes. arXiv e-prints, 2013.*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-uTOD3zMYzTb"
      },
      "source": [
        " from torch import Tensor\n",
        "\n",
        "# create a tensor, let torch know we want to compute gradients\n",
        "x = torch.tensor([1., 2., 3.], requires_grad=True)\n",
        "\n",
        "# do some kind of transformation\n",
        "def dummy_transform(z: Tensor) -> Tensor:\n",
        "  return torch.sin(z) ** 2\n",
        "\n",
        "# create a random SPD matrix and add the transformed b to its diagonal\n",
        "a = torch.rand(3, 3)\n",
        "M = a @ a.t() + torch.diag_embed(dummy_transform(x))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zo_J_eHyajaP"
      },
      "source": [
        "Let's draw some samples from a Multivariate normal.\n",
        "\n",
        "Note that if $\\epsilon_i \\sim \\mathcal{N}(0, 1)$ iid and $LL^T = M$, then $\\xi := L \\epsilon \\sim \\mathcal{N}(0, M)$. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kHN_JZWFZDBR",
        "outputId": "db977bdb-3515-4924-8c7f-6345b055287e"
      },
      "source": [
        "# compute the Cholesky factor of M\n",
        "L = torch.linalg.cholesky(M) \n",
        "\n",
        "# draw quasirandom samples epsilon\n",
        "t_sobol_engine = torch.quasirandom.SobolEngine(\n",
        "    dimension=3, scramble=True, seed=29211402,\n",
        ")\n",
        "sobol_samples = t_sobol_engine.draw(2)\n",
        "epsilon = torch.distributions.Normal(0, 1).icdf(sobol_samples)  # 2 x 3\n",
        "\n",
        "# correlate samples \n",
        "xi = epsilon @ L.t()  # 2 x 3\n",
        "xi"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "tensor([[ 2.8662,  1.9570,  2.1361],\n",
              "        [-1.4842, -1.0919, -1.4085]], grad_fn=<MmBackward>)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 58
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4qEFutgScM7B"
      },
      "source": [
        "We can now compute gradients of these samples w.r.t. our input using PyTorch's auto-differentiation feature"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "fL5sQzh_cd0y",
        "outputId": "42267098-7e0c-4360-94cb-51e89a6701c1"
      },
      "source": [
        "# compute gradient of the first element of the first sample w.r.t x\n",
        "torch.autograd.grad(xi[0, 0], x)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(tensor([1.1680, -0.0000, -0.0000]),)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 59
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FUia92aphAfK"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}