scoremixin_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fisher Score Mixin Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates how to use the ``qinfer.ScoreMixin`` class to develop models that use numerical differentiation to calculate the Fisher information. We test the mixin class with two examples where the Fisher information is known analytically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preamble"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import division, print_function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.\n",
      "  warnings.warn(self.msg_depr % (key, alt_key))\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/qinfer/metrics.py:51: UserWarning: Could not import scikit-learn. Some features may not work.\n",
      "  warnings.warn(\"Could not import scikit-learn. Some features may not work.\")\n",
      "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/IPython/parallel.py:13: ShimWarning: The `IPython.parallel` package has been deprecated. You should import from ipyparallel instead.\n",
      "  \"You should import from ipyparallel instead.\", ShimWarning)\n",
      "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/qinfer/parallel.py:53: UserWarning: Could not import IPython parallel. Parallelization support will be disabled.\n",
      "  \"Could not import IPython parallel. \"\n"
     ]
    }
   ],
   "source": [
    "from qinfer import ScoreMixin, SimplePrecessionModel, RandomizedBenchmarkingModel\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "try:\n",
    "    plt.style.use('ggplot')\n",
    "except:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Precession Model Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create two models, one that uses ScoreMixin's numerical ``score`` method, and one that uses ``SimplePrecessionModel``'s analytic ``score`` method. To make the first model, we declare a class that does nothing but inherits from both the ``ScoreMixin`` class and ``SimplePrecessionModel``; note that ``ScoreMixin`` is first, such that its implementation of ``Model.score()`` overrides that of ``SimplePrecessionModel``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class NumericalSimplePrecessionModel(ScoreMixin, SimplePrecessionModel):\n",
    "    pass\n",
    "\n",
    "analytic_model = SimplePrecessionModel()\n",
    "numerical_model = NumericalSimplePrecessionModel()\n",
    "\n",
    "expparams = np.linspace(1, 10, 50)\n",
    "modelparams = np.linspace(.1,1,50)[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We verify that both models compute the same score by plotting the score for a range of experiment and model parameters. Since this is a single-parameter model, the score is a scalar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(50, 50)\n",
      "(50, 50)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fc9bf549470>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAC8CAYAAACOhp56AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnV1slNfd4H/nnGdmjG3wxNQkBMJL3rJZisu2CWE3CtlA\nwkWk3MSqFEtpFCVqL6pKbYnVLpBGiF70IkFJAIWIq6qNmptyg6VI1V6sBNsWKRJ+gW5qym7ZzZsP\n0gRjg8EYz8zznLMXz8x4/IU94/l6nvn/JH/M53nmmd/85/zPp3LOOQRBEIRIoht9AIIgCELlSBAX\nBEGIMBLEBUEQIowEcUEQhAgjQVwQBCHCSBAXBEGIMMsK4hcuXODVV19lz549DA4OLvlxw8PDyyl2\nWTSqbHnN0ULclnKbtezZVBzErbX85je/4fXXX+ftt9/mzJkzXLlyZUmPbcWTL685OojbUm4zlz2b\nioP45cuXWbt2LT09PXiex44dOzh79mw1j00QGoK4LUSJioP42NgYq1evLl7u7u5mbGysKgclCI1E\n3BaihKp02v1HH33EX//6V370ox8B8Kc//YnLly/zgx/8oKoHKAj1RtwWooRX6QO7u7u5du1a8fLY\n2Bjd3d1z7jc8PDyj/ai/v5+/uAxnXKbSoivHwQ6V4i8uCw5A1a3oJ3SSv9hs3cprhrIbUe5+r5MT\nJ04UL/f29tLb21vWc4jb5dFKfjWy7IXcrjiIb9q0ia+++oqRkRHuuecezpw5w549e+bcb74P0RmX\n4ZC7VWnRZeMc4MA5xT6jeMOfyN9SP9H308kbwcTid4xR2Y0od7/XSX9//7KeQ9wuj1byq5FlL+R2\nxUFca80Pf/hDfv3rX+Oc4+mnn2b9+vXLOsha4lA4p/K1FEFYGHFbiBIVB3GA7373uxw9erRax1JV\nSlv6ZwquqGctRYgm4rYQFZYVxJuefJrpkFqKEDPEbSFPrIN4Mc0UhJghbgsFYhXE755mCkJ0EbeF\nhYhVEAfCNFM6eoQ4Im4L8xC7IO5QOCu1EyF+iNvCfEQ+iEuaKcQVcVtYCpEP4kCxpx6Qzh4hXojb\nwiLEIohLT70QV8RtYTFiEcQlzRRii7gtLIJszyYIghBhJIgLgiBEmEg2p8zptReEmCBuC+USySAO\nFCc+lPbeC0IsELeFMohsEJeJD0JcEbeFcohsEJdeeyG2iNtCGUjHpiAIQoSRIC4IghBhItOcIr32\nQlwRt4XlEJkgDkivvRBfxG2hQiIVxGUtZSGuiNtCpUQqiOPyP9JrL8QNcVuoEOnYFARBiDCL1sSP\nHz/OuXPn6Orq4q233gJgYmKCI0eOMDIywpo1axgYGKC9vb3mBysI1UTcFuLAojXxp556itdff33G\ndYODg2zdupWjR4/S29vLyZMna3aAglArxG0hDiwaxDdv3kxHR8eM64aGhti5cycAu3bt4uzZs7U5\nOkGoIeK2EAcqahMfHx8nnU4DkE6nGR8fr+pBCUKjELeFqFGV0SlK1aZHvXQShAy9am1cg95/cVuo\nNct1u6Ignk6nuXHjRvFvV1fXgvcdHh5meHi4eLm/v58dKsXepRSkmJ4Ekf9/+obKeEIn2U9nxY+P\nWrmNLLua5bo5/8y6vWSCzIkTJ4r/9/b20tvbu+RyxG0pt95lL9ftJQVx5xyu5Oti27ZtnD59mr6+\nPk6fPs2jjz664GPn+xCdcRkOuVtLKDd8AXMnQVQu+n46eSOYqPjxUSu3kWVXs1yXH0dd9GHGbaro\nyi+THfT395fxvOL2cmm1cqtd9nLdXjSIHz16lIsXL3Lr1i1+/OMf09/fT19fH4cPH+bUqVP09PQw\nMDBQlRczB1fyI5MgYsuS0kkH1imc0yjraGOKNjIkyTBFG1OujSlSZZUrbgu1ph5uLxrE9+zZM+/1\nBw4cWMLRCcLSmXfNEFeYkg7OaqzVJGyO9Vzhm+7/8QBX+MT8C5+aDXxqHiirPHFbqBe1dDta0+6F\n2FJoVnCltVI3M520VmMDTdI6HuAK/8UN8QgX+Ehvx2nF14k1jXsBgrAAtXZbgrhQN+6aWhakDoBs\nANkAFzhIGEgaMAacAqfQztLlxtngPudbXOIKa/k/ehMJL1evlyIIM2ik2xLEhbozJ7V0YK3GWYWd\n8OHTm7hPx1G3MpgHOzEPdqLXtmNdmHIaF/BPdy9DPMItOvlbYgtfmvvxRWehwTTCbbFeqCvzp5bh\nxsDWauxti7t8E/fRl6ivb+I9sYZE2pLYYLBoAqfRzvIV9zGE4f/yIKN6Ndf0agJM416Y0PI0ym0J\n4kJVmC+ddLOvzwtdkF0ph1Lkh/iF17uMxX49ifvfY6jPrqMeTOLd6SDhdWJRGDQWzQiruUrPnHZG\nGeghVJtmd1uCuFBVZqSTbmbtxAUKGxicb1A4EoksXiKLNpbAGQJjCLosbnMHLnMvjHZivtuN6+kg\nwGBROHQ4HAsFqJmigwRxoWY0q9sSxIWqMTuddEzXTpwDGxhsNoHNehgVoHWGZDKH5+XwncHHw3Y5\n7OZ2XPe9cMdH39eG62kjwOBQxfG0xTIlagt1oJndbrog3qg1MoSlseD7UzrrzOZnIWpwOR3epBQu\n0FjfYLMeWoNKQULlSJgMGg+NJVhpcCtX4L7ZUfKBgYDSGW3RDNzidnMTVbebLogXmLeTQGga5uuF\ndy5fm7ijceMa22GwVz10l4+XzkHSJ1AW6/loZVGpgMAYIEGACdNKpxdOJ2OCuN3cRM3tpgziRcll\n1++mZKHJC7bQuTOpsVc9bNpg/93D+xcfb6WP7rBYLyCwBgVozxJojSOBRYfDrEpqI3EMcuJ2cxNF\nt5sviJekLlFNm+OGK/5i5vtjHco5lLVYNKBBKeyUxl7zcBjsFx4qDYn1PibhYxN+KHV+KftCLSXq\nTSVLQtxuOuLgdvMFcaH5KJHcOl0c99oxMsKqK1+w8ssr+J3tTKy/n4n1a8mu7CRYZzErNYlNU+ie\nAJdSRcEtWmqiQnMQA7cliAtLophmFiYuWE3q6nW+8deLrPu3s0zd183Vxx7G3tsOKxNYfLxUMuzY\nWeXjUhQFd8xdclMQGkXU3ZYgLhS5a++8ne6dt0EoetvIGD3/a5gH//v/4Pam+7H3djLxyIMEq1Zi\nV2mMayexOpN/CjX9I80JQp2Js9sSxIW5FDp2XL6Z8GuD/VLjvjSotMWstSTW+gRrV3L94Yf4zO4i\nc183NzY9SHZFR9iJQ6EnvuRppfYtNJoYui1BXJhB6eiJ4kSGfxrsOQ97zsP71xzedp/EfZlQ9G2b\n8deuItfZzu0H1pJtb59uF1Qz5ZYgLjSSuLotQVyYpqR33hZnoymCrzzs+QTBH5OYR8PaSurRLMH9\nnVy/7z8y9sjmcMKDUqD1jNRSArfQFMTYbQniLcbddlm3vin+aBNgTIDyfMx9DvufLEHG4v2HLPoB\ni0sqnDJYMz2karoMVaytSNu3UC9a1W0J4i3K3A16Q9GDTIIgkyCRyqBX+CS9LHZdgN3uE9yXw3wj\nQG0IsMmwA2d6Ftqs5xaEBtFqbksQb0Hm32U9L/pUAn8yiad8TJsl6eVw63IE9xnswxo04IEzs3vk\nBaHxtKLbEsRbgNK1j4udOjb8MToo/ugE6BUOrRxeWw7l2bA9UMPsNenntgs2v+xC/ChtQmlVtxcN\n4qOjoxw7dozx8XGUUuzevZtnn32WiYkJjhw5wsjICGvWrGFgYID29vZ6HLNQIaWSF7aM8rSPp3Kk\nTBaTsuRMgEkGmISPSlisUvllMqORWpaDuB0fWtntRYO4MYaXX36ZjRs3MjU1xb59+/jOd77DqVOn\n2Lp1K8899xyDg4OcPHmSF198sR7HLFTAbMltEIqujCWhc7R5dzAmwKR8cs4PKx/KFcfERiW1LAdx\nOx60utt6sTuk02k2btwIQFtbG+vWrWN0dJShoSF27twJwK5duzh79mxND1QoD+emf3Dgshp72yO4\nlYBMuJefl/DRxoJW4VKZxaFUDlSYp86fWs73Ez3E7WhSdJtwqnyru11Wm/jVq1f59NNPeeihhxgf\nHyedTgPhh2F8fLwmBygsj+J42IyHve1hpwxeRw7T4eMlsxhtsVqRI0mQXzIzyqllpYjbESS/3kmr\nu73kID41NcU777zDK6+8Qltb25zblYr+yYgbpT31dsoQ3EoQ3E6EqWVHQDKZDWeeockxPZ04DmKX\ng7gdPYqLTVnV8m4vKYgHQcDbb7/Nk08+yfbt24GwhnLjxo3i366urnkfOzw8zPDwcPFyf38/O1SK\nvQuU5XT+JDuoRSrzhE6yn86qP2+jynULXijsH+LYoT0OdoPtUFjfYlIepq0DrVL5e4X3Dp9CVe20\n3+19riUnTpwo/t/b20tvb++C9xW3m7fcxdx+QnuYthS227S028q5xXf+O3bsGCtXruTll18uXvfB\nBx/Q2dlJX18fg4OD3L59e8mdP2/amxxyt4qX3aydo6ffsOqLvt908kYwUfXnbVS5pUMHS9eFUFiS\nJktKZ/lvJs2bkxNkMylyuSSmzUe35TBtQfjYguQzainLP/d71coZ73M9GDXryrq/uN285S7m9i8T\nHRwNrpHNpFra7UVr4pcuXeLPf/4zGzZsYO/evSileOGFF+jr6+Pw4cOcOnWKnp4eBgYGKj+6GTue\nCOUyvaCPwuWX1dQojAlI6QwJk2VF+230igBl/bA7W4WL4Lcy4nbzcze3jQpoT9zGGL+l3V40iG/e\nvJk//OEP89524MCBqh9QFHqDG82cCQ4lC/ooHFpbtArQ2hFWYhTOKMIWQlvz7aKigrjdfJTjNgqs\n0i3vtszYjCpu5hApG4Q1D20snvHxtI8yjkB7+HjkSBA4E9vOHSFGLNFtqzQZUvh4Le22BPGIUpS8\nZEspgITJ4RmfhJfDofCVR4Ahl991uxUlF6LFUt22hEG81fdtlSAeEeZLMwvTiwG0tijl0MYW/y+I\nbgtjZFsszRSiQaVuO1RddpNvdiSIR4lCmmkVQRCujewcGC/AGB/PCxf7QRHKnV8XYvYuJILQdFTg\ndrNtztAoJIhHiIK01mqsr/F9g3LgeT6e55NI5MKAnZ/k0KjdtwWhXCpxWyonIRLEm5gZaaZVuMLi\nPoEGp9DKoXSYYmptUdrlx8XOWg+5iXYhEQSootvitQTxpqeYZmqCnIfNeuFYWc8nmciivQDjBSVj\nY9UCExwEoclYrtuiNyBBvOkpppmBxmY9/KkkChf20idyeMkcYYVEFUUXhCggblcHCeJNxpw0MwjT\nTOubcI1k5dA631NvArQXlIw8aZ1ZakL0ELdrgwTxZqSQZgaaIOthMx7OKrRxmLZMmGYmc6DJSz53\niU1BaErE7aojQbwJmZFmZjz8yTDNNB1ZEm1ZTCqHKkw7ljRTiBDidvWRIN4ELNRTX0wzAaVdvpbi\nY5J+yfhYSTOF5kXcrj0SxJuF+XrqC2lmRz7NTOVAO5ngIEQLcbumSBBvEubtqVdhO2GiLYtJ+GBc\nOF5WJBcihLhdWySIN4iF1ouwwfRCPmGaaTGp3Kw0U0QXmhdxu75IEG8kpWlmEE43xukwvcwPswp7\n6iXNFCKGuF03JIg3kOk0U2F9g+8btHJ4CZ9EwkebAKUdStaIECKGuF0/JIg3CjedaoY1kfA6pV24\nclsqhzGBpJlC9BC364oE8QZRWHazsOC98QI846ONxXgBCkkzhWgibtcXCeKNIl9bsYFG63DbKeP5\nxb0DZZlNIbKI23WlYUF8Rg92i6RTC73mcNeSgEQinK1mncY6g6SZ0UTcFrfrSWNr4vke7MI3d6tQ\n2Lnb6ACTDMJ1k43NL3bvcE5mqkUecVvcrhOLBvFcLsfBgwfxfZ8gCHjsscd4/vnnmZiY4MiRI4yM\njLBmzRoGBgZob28vq/Diwu5u8fvGhUJboFaFBe8DVCHFzG9d0iq1t0YjblcXcbsxLBrEE4kEBw8e\nJJVKYa3lwIEDPPzww3z00Uds3bqV5557jsHBQU6ePMmLL7649JJdyU+M39hCmunyv1T+U62VxZhw\n/8DCYj+SZi4NV6XAKG4vD3G7+lTi9pJym1QqBYQ1lyAIABgaGmLnzp0A7Nq1i7Nnz5ZfeivhQGHx\nlE+bmSKps8WNXwvbTgll4qZrf5WOdhC3q4C4XX3KcHtJbeLWWvbv38/XX3/NM888w6ZNmxgfHyed\nTgOQTqcZHx+vzsHHkMKGxRqHp3MkdRaUIlCGAF3cM1AoAzd9XpeDuL08xO0aUKbbSwriWmsOHTrE\n5OQkb731Fp9//vmc+yg1f4HDw8MMDw8XL/f397NDpdinpzt+6vlN/YROsp/OupXn8m/If1UJjLcS\nT/t42scBFkOACc9DsdZSfXaoFHtr8swNKrd0iNo8NZQTJ04U/+/t7aW3t3fBpxK3K6fRbjfK65qW\nXYHbZY1OaW9vZ8uWLVy4cIF0Os2NGzeKf7u6uuZ9zHwfojMuw5v2VkNmbO2nkzeCibqVF+QMNmd4\nLelxKJOhre0OqdQd0Aofj8AZLLokXar+udgLHHK3qv689Sy32FY4Z9THzPO13+ukv7+/7OcXt8un\n0W43yutql71ctxdtE7958yaTk5MAZLNZPv74Y9atW8e2bds4ffo0AKdPn+bRRx9d1guJKy5Q2KyH\nzRn8TIpcLkXOJReQXLgrs9oJl4u4vTzE7SqyDLcXrYnfuHGD9957D2stzjkef/xxHnnkER566CEO\nHz7MqVOn6OnpYWBgoOLjjxszdzPR2JzBGY2fTWBSPoYECjs9DE3aDOdldk99sX21SoFB3C4fcbs6\nVNPtRYP4hg0bePPNN+dc39nZyYEDB8ousFUofKMq7fJrJidIuAw64eeHZMnaEUumJM2s5rhrcbsy\nxO0qUgW3Ze2UGjC9EL5CaYtJhYvfJ0wG7dnifYQlUNpT30ITZ5oVcbuKVMltCeJVYuZuJmBduIqb\n8QJMwsckfBLJjCy/OQ9uzj+lt1W3CUUoH3G7curhtgTxauNAK4fRPuhw30CtwndQlt9cgNm987Nv\nkxp4cyBul08d3JYgXkUKa2UYFYTTjpUPSmGVLhk2JMzLQhMcJIA3BeL2Mqix2xLEq0X+W7XQ4eNp\nn6SXxaLxnVfyBors8zF3NEOp4XLOGoq4vSxq7bYE8SphXbibtw00yjmMCooz1mxxOL5IXmC+CQ4z\nkXPVLIjb5VFvtyWIVwlnCxvCemjl8E2AcQFOydoRCyIjTyKBuF0BdXRbgniVcFYRBAY/56F1uJdg\ngB/eJqLPy3SaKTQz4nb51NNtCeIV4oq/wr9KOzwvh1YBXiJAm2DmOgji+qyhaqU1FDk5zYS4XT6N\ndFuCeKUUFsTPv2FaW0wyAOfQJhx+JbXMeZAmlOZH3K6MBrktQXw5lPTaGxOEu3obP+zwcRq3tD03\nWgppQokI4nbZNMptCeLLwE5o7C2DvWVgVQ7d5dBdAcWNBYFWzzWlCSWaiNuL0yxuSxBfBu6mwX6Z\nwL+SQN0PGoteqUEraSkoRZpQIoe4vUSawG0J4svA3dQEX3j4f0+hA4dZGWDv91FaNC9FmlCih7i9\nNJrBbQniZTBjDWAHaoXFdPuwPoNZ7aPabf7G1g5YxfOkpAklKojbS6MZ3ZYgXgGFN0+tsngPZPFW\n+eh0gEoHS9grqUUoDD9zas4C+ELzIm4vgSZzW4J4ubjp8bF6VYBeZTFYnAYM+S/k1q6twPR042aQ\nXFgi4vaSaDa3JYiXgcbSbcdY7V/nnuA61xNpbiTTXE+mafXdTObtqZcPfWQQtxem2d2WIF4GGst6\n/wpbsxf5VuYSwyu+xUX9La57aZl+DCUL/jRHmiksHXF7EZrYbQniZaCxPBBc4T9nzvL0nf9Ju55k\n3Oti2H2rWb6UG4oj3LpLiB7i9t1pZrcliC/C7FRqUrUz5nXzZXItY949TOoV+Vub8w2uJYtPdmi9\ncxIlxO2FiZLbSw7i1lpee+01uru72bdvHxMTExw5coSRkRHWrFnDwMAA7e3ttTzWhuKcInCGL7z7\nGVKPcD2R5nLim/zTW9va6aYr7eiJ3nloda9B3F6QiLi95EFDf/zjH1m3bl3x8uDgIFu3buXo0aP0\n9vZy8uTJmhxgU5BfRyIUfR1nV2zjw1XP8m8rHuZKorVFL+x+3syS342W9hrE7bsQFbeXFMRHR0c5\nf/48u3fvLl43NDTEzp07Adi1axdnz56tzRE2AdbX2CkPfzLBRHYlV10Pn3vruWZWM6namymzqjnO\nlf6UvnBV8hMNWt1rELdLiarbSwri77//Pi+99BJKTb+I8fFx0uk0AOl0mvHx8docYaNx4LKGYMIj\nuJ4iuJ3AZU1r7+7tptsJo3wOWtprELfnI4JuLxrEz507R1dXFxs3bsTdZWxN6QchVjhwGU0wkcDP\ni27zorcqUUkz70bLew3i9jxE0e1FOzYvXbrE0NAQ58+fJ5vNcufOHd59913S6TQ3btwo/u3q6pr3\n8cPDwwwPDxcv9/f3s0Ol2KvzJ8lBPdOUJ3SS/XTe9T4zPtIa7EqDSxpsGnRbAtOm0CpZVrk7VIq9\n5R9uVahm2Q4Vzt5bwvu2lHNdC06cOFH8v7e3l97e3jn3Wa7XIG4XaJTb1S43qm4rd7dqyCwuXrzI\nhx9+yL59+/jggw/o7Oykr6+PwcFBbt++zYsvvrik53kjuMWbwa2Sb7v6ib7fdPJGMHHX+xTOiHMK\nLLiMwWU0LqsxK3x0u49ZEZRV7l61kkPuVqWHvSyWW7YrbBAwp6f+7u/bUs51tbmRuq/sx1TLaxC3\n60k1yo2D2xUvadPX18fHH3/Mnj17+Nvf/kZfX9/SHxyF9qZCuxgKnQpIrMqQ+sYk3qosOlme5HEg\nimlmJSzLaxC3I0jU3S5rss+WLVvYsmULAJ2dnRw4cKDCYpu3p7dA2EOtw29qz6JMuMv3dKdPcx//\ncpk92WF6qFn8Xnf1vAZxu/mJm9syY3MBnFNYG35DK+3Q2rZer31JminbucQHcZtYuS1BfAFC0TXW\napR1WKfRzi7+wLhQ2HaqlT7YLYK4HS+3JYgvgFIOpRxaWaZ3FYzHm74QM9JMop9mCvMjbsfLbQni\nC6C1RauwdqJ0KH1LEKM0U5gfcZtYuS1BfB6UCtsJlbJo5Uq+uaP/rX1XYpZmCnMRt+P3OiWI5ylN\ntxQUO3yMtlinsU7HbjGg2TME4pZmCiHidrzdliA+Gxe2ElqnUc5hXUGAGBPTNFOYhbgdSySIl1Lc\n/DQcfmWVnpYgrsQ4zRRKELdjiwTxAg5cVmOzCgIFKSAFKtXoA6stUU8zm22/w6ZE3CbObksQL+Ar\n7JgmuObhbhu8nhxqTQ5SrTJ+NnqSzyDutcrlIG43+gCWxyJuSxDP4wKwYwb770mCUYPKOXR7AN2N\nPjJhUVokba4UcTvCLMFtCeIFLDClcDcVblTDvQqyjT6o6uOKv6LL7DRTxb97bnmI25GhErcliOdR\nHpjVPupfwftGDnN/Dr0yhummm/47s9feEbm0M/9aPOWTUDkSKgesbeghNSPidrzdliCeR3kO8w0f\n0x5AFlSnRXVE/Gt9IQopWkH0KFKy5KunfFboO7SbyQYfVHMibkeMMt1u6SA+YxKEcegui+62KOVi\ntyznvGtHRKgNeb6e+kKq6SmflM7QZqbqf2BNSqu4veCknhZyu6WDeJHCpIf84vCVb5XR5MRh4kP+\nuI0K8JSPp3yMDsipBJO0N/bYmhFxOzpU6LYE8ULq4sApVfw/dsQozSytoSRNFqcUvvLIkWj0ETYX\n4nZ0WIbbLR3EC72+xR5gN31L3IhLmgnh+2VUQFJnadNTTNGGj0fGpeL41lWEuN3cVNPtlg7iWlna\n1SQruIOHz5RqY0q1kSXmU9mixqw006gAT/tYpZmijRwJAkxY04zO57imiNsRoQput3QQNwR0qAm6\n1Rgr1B2u08117hHRm4n50kydI1CaQJui5IEzMtmnBHE7AlTJ7ZYP4u16km49RqeawHcet11How+r\nKsxI1yIU2+6WZnrKJ2UypEyGKZdPM0nFarRFtWgJtyP2dtfK7ZYO4gGGSToYYzWTtHOTVWRJNvqw\nqktJ51akOn3cdPtgaU+9VZqMS81MM4U5iNtNTJXdbukgbtHcpgOHIkGOO6wgE6d0s2ThnEgFu5I0\nM6FytOkpkjpLoAy+MmRJShPKIrSE21EccVMDt1s6iIe1lXamaENRulVVPCiMDyZiTQ3F41aQUDlW\nmElSJsMd2smSIOOkCWUxWsHt8J9ovaZauB3Xof+CIAgtgXJOltUXBEGIKg2piZ84caIRxTa0bHnN\nrYG8z/Evt9Flz0aaUwRBECKMBHFBEIQIY371q1/9qhEFr1mzphHFNrRsec2tgbzP8S+30WWXIh2b\ngiAIEUaaUwRBECKMBHFBEIQIU/cZmxcuXOB3v/sdzjmeeuop+vr6albW8ePHOXfuHF1dXbz11lsA\nTExMcOTIEUZGRlizZg0DAwO0t1d3R5jR0VGOHTvG+Pg4Sil2797Ns88+W/Oyc7kcBw8exPd9giDg\nscce4/nnn6/Lawaw1vLaa6/R3d3Nvn376lZus1Avt1vNaxC374qrI0EQuJ/85Cfu6tWrLpfLuV/8\n4hfuiy++qFl5f//7390nn3zifv7znxev+/3vf+8GBwedc86dPHnSffDBB1Uv9/r16+6TTz5xzjl3\n584d97Of/cx98cUXdSl7amrKORee61/+8pfuH//4R13Kdc65Dz/80B09etS98cYbzrn6nOtmoZ5u\nt6LXzonbC1HX5pTLly+zdu1aenp68DyPHTt2cPbs2ZqVt3nzZjo6Zi6/OTQ0xM6dOwHYtWtXTcpP\np9Ns3LgRgLa2NtatW8fo6Ghdyk6lwkWOcrkcQRAA9XnNo6OjnD9/nt27dxevq0e5zUI93W5Fr0Hc\nXoi6NqeMjY2xevXq4uXu7m4uX75cz0NgfHycdDoNhFKOj4/XtLyrV6/y6aef8tBDD9WlbGst+/fv\n5+uvv+aZZ55h06ZNdSn3/fff56WXXmJycrJ4Xb3PdSNptNtx9xrE7YVo+Y5NpWq3CtrU1BTvvPMO\nr7zyCm1JVnmDAAAB1ElEQVRtbXUpW2vNoUOHOH78OJcvX+bzzz+vebmF9tmNGzfi7jJitZbnWphJ\n3LwGcXsh6loT7+7u5tq1a8XLY2NjdHd31/MQSKfT3Lhxo/i3q6urJuUEQcDbb7/Nk08+yfbt2+ta\nNkB7eztbtmzhwoULNS/30qVLDA0Ncf78ebLZLHfu3OHdd9+t6+ttNI12u1W8BnF7NnWtiW/atImv\nvvqKkZERfN/nzJkzPProozUt0zk34xt027ZtnD59GoDTp0/XrPzjx4+zfv16nn322bqVffPmzWLK\nl81m+fjjj1m3bl3Ny/3+97/P8ePHOXbsGK+++irf/va3+elPf1q3c90M1NvtVvIaxO27UfcZmxcu\nXOC3v/0tzjmefvrpmg4xPHr0KBcvXuTWrVt0dXXR39/P9u3bOXz4MNeuXaOnp4eBgYE5nUTL5dKl\nSxw8eJANGzaglEIpxQsvvMCmTZtqWvZnn33Ge++9h7UW5xyPP/443/ve95iYmKj5ay5w8eJFPvzw\nw+IwrHqV2wzUy+1W8xrE7bsh0+4FQRAiTMt3bAqCIEQZCeKCIAgRRoK4IAhChJEgLgiCEGEkiAuC\nIEQYCeKCIAgRRoK4IAhChJEgLgiCEGH+P/WF9sPNoHx1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc9b3ccf630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "analytic_score = analytic_model.score(np.array([0],dtype=int),modelparams, expparams)[0,0,...]\n",
    "print(analytic_score.shape)\n",
    "numerical_score = numerical_model.score(np.array([0],dtype=int),modelparams, expparams)[0,0,...]\n",
    "print(numerical_score.shape)\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(analytic_score)\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(numerical_score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we verify that both models give the same Fisher information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fc9bf450d30>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAC8CAYAAACOhp56AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFtVJREFUeJzt3VtoXNe9x/HvvszsseygqYIUErtGFNcYT02b2j6EBOJJ\n/BDwS0QhA6kJDe1D6aGtI1rqpMHYgj60wY4t4qCn0ib1S/ViBUNfPWkxBCRs97hyfU5MLnUKiWWp\nGus+t30e5qKt8R6P4mhGWtLvA3/W3muvPWtp9l//TEYz25bv+z4iImIke7UXICIiD05FXETEYCri\nIiIGUxEXETGYiriIiMFUxEVEDPaVivjVq1d55ZVXOHLkCENDQ8s+b3R09KtM+5Ws1tz6mc2i3Na8\na3XuWg9cxIvFIr///e95/fXXOXXqFJcuXeLf//73ss7diE++fmZzKLc171qeu9YDF/GbN2/y6KOP\n0tnZieu6PPXUUwwPD6/k2kRWhXJbTPLARXxiYoKHH364ut/R0cHExMSKLEpkNSm3xSTWg37t/oMP\nPuDvf/87P/7xjwH461//ys2bN/nhD3+4ogsUaTXltpjEfdATOzo6uHPnTnV/YmKCjo6Oe8aNjo4u\nef8olUqRTn9COv0JYAWCkH27vG2HHLdqjldap9wG+xYjmXyIdHouMC44xqnpL0dweU7I0NrtkL7k\nHkj/b6DfofTshz2OG4iwsa5fjiKW42M7BRy3gG0Xca08jl3AsQq45HEosJ+H+B/u4FAoRx63Ztsl\nXw0nsF9/u9RGKvuFIm6+iFMoYud8yAGxJPwnXdquRL7cZlnaH4zyMT8LfrmvkINsARaKpcgCC+XI\nBtrkiRMMDg5W8y2RSJBIJPgylNv32V5Duf1fbOHvjAfyemPm9gMX8R07dvD5558zNjbG1772NS5d\nusSRI0fuGRf2S5ROf0Jf3/vcm6jBpKtkVSWsmn675nglGyKEZ0nlmE1f30RgXHBMtByRpfMFlxYp\nRzSwHbZf22dD33vlPq+mDY7zgFg5vJCxHuD5ECtCrIDlFYhEs0S8HJFIDs9eIGpliVpZYszjsYDL\nY5zjQ6Jky7GAV20XiFXb+eo5lf17t+fL40tjN1XOKeTwsnnshTzM+TALtAO3+krbcyxtZ2raYMws\nhj8L/gzkZ2EmB1N5uJuHKUpxN9BOU0r0VCp1Tx5+GcptjMhth0d4l48Ceb0xc/uBi7ht2/zoRz/i\nN7/5Db7v8+yzz7Jt27YHfTiRNUO5LSZ54CIO8J3vfIf+/v6VWovImqHcFlPoG5siIgZTEa/Lajxk\nLfqSy77/cEOfA2nA0Ouq3A6lIi4iYjAVcRERg6mI12XoPz36JZd9/+GGPgfSgKHXVbkdSkVcRMRg\nKuIiIgZTERcRMZiKuIiIwVTERUQMpiIuImIwFXEREYOpiIuIGExFvC5D762g+0tIQ4ZeV+V2KBVx\nERGDqYiLiBhMRbwuQ++toPtLSEOGXlfldigVcRERg6mIi4gYTEVcRMRgKuIiIgZTERcRMZjbaMDA\nwACXL1+mvb2dkydPAjA9Pc2ZM2cYGxujq6uL3t5e2tramr5YkZWk3Jb1oOEr8WeeeYbXX399Sd/Q\n0BB79uyhv7+fRCLB+fPnm7ZAkWZRbst60LCI79q1i82bNy/pGxkZ4cCBAwAkk0mGh4ebszqRJlJu\ny3rwQO+JZzIZ4vE4APF4nEwms6KLWhsMvbeC7i/xlSi31zDldqgV+cOmZa2fJ0QkSLkta13DP2yG\nicfjTE5OVtv29va6Y0dHRxkdHa3up1Ipksnu8p4ViLB9u7xthxy3ao5XWqfcBvsWI5ncXDMuOMap\n6bcWp620TsjQ2u2QvmSC0rPtBMKt8zhuIMLGuha4NrhgOTa2Y+O4EWy7iGtvwrEKOBRwyeNQ4HEe\nwuGbOBTKkcdl6RiXfDUW9ytja/vz1f5IZd8p4kaLWE4RIj5sAmJJ+DqQC0S+3GZr+oORXWytXCnc\nHGwugFuEtiJ0AAvlyAZagMHBQSoSiQSJRILlUm5jTG7vZQsu3wjk9cbM7WUVcd/38f3Few3s3buX\ndDpNT08P6XSaffv21T037Jconf6Evr73uTdRg0lXyapKWDX9ds3xSjZECM+SyjGbvr6JwLjgmGg5\nIkvnCy4tUo5oYDtsv7bPhr73yn1eTRsc5wGxcnghYz3A8yFWhFgByysQiWaJeDkikRyev0DUyRK1\nssSYx2MBl8c4x4dEyZZjAa/aLhCrtvPVcyr7927Pl8eXxm6qnFPI4WXz2At5mPNhFmgHbvWVtudY\n2s7UtMGYWQx/FvwZyM/CTA6m8nA3D1OU4m6gnQaSJ06QSqVCczGMchtjc9vhEd7lo0Beb8zcbljE\n+/v7uX79OlNTU/zkJz8hlUrR09PD6dOnuXjxIp2dnfT29jZ6GJE1R7kt60HDIn7kyJHQ/mPHjq34\nYkRaSbkt64G+sSkiYjAVcRERg6mIi4gYTEVcRMRgKuIiIgZTERcRMZiKeF2Gft1a95eQhgy9rsrt\nUCriIiIGUxEXETGYinhdfuMha9GXXPb9hxv6HEgDhl5X5XYoFXEREYOpiIuIGExFXETEYCriIiIG\nUxEXETGYiriIiMFUxEVEDKYiXpehX8vVV5OlIUOvq3I7lIq4iIjBVMRFRAymIl6XoV/L1VeTpSFD\nr6tyO1TDf+1+fHycs2fPkslksCyLgwcPcujQIaanpzlz5gxjY2N0dXXR29tLW1tbK9YssiKU27Ie\nNCzijuPwgx/8gO7ububn5zl69Cjf/va3uXjxInv27OH5559naGiI8+fPc/jw4VasWWRFKLdlPWj4\ndko8Hqe7uxuAWCzG1q1bGR8fZ2RkhAMHDgCQTCYZHh5u6kJFVppyW9aDL/We+O3bt/n000/ZuXMn\nmUyGeDwOlH4ZMplMUxYo0grKbTHVsov4/Pw8b775Ji+//DKxWOye45a1fj53KRuLcltM1vA9cYBC\nocCpU6d4+umn2b9/P1B6hTI5OVlt29vbQ88dHR1ldHS0up9KpUgmu8t7ViDC9u3yth1y3Ko5Xmmd\nchvsW4xkcnPNuOAYp6bfWpy20johQ2u3Q/qSCUrPthMIt87juIEIG+ta4NrgguXY2I6N40aw7SKu\ntQnHKuBQwCWPQ4HHeQiHb+JQKEcet2bbJV8NJ7Bff7vURir7ThE3WsRyihDxYRMQS8LXgVwg8uU2\nW9MfjOxia+VK4eZgcwHcIrQVoQNYKEc20AIMDg5SkUgkSCQS1KPcDkxvWG7vZQsu3wjk9cbM7WUV\n8YGBAbZt28ahQ4eqfXv37iWdTtPT00M6nWbfvn2h54b9EqXTn9DX9z73Jmow6SpZVQmrpt+uOV7J\nhgjhWVI5ZtPXNxEYFxwTLUdk6XzBpUXKEQ1sh+3X9tnQ9165z6tpg+M8IFYOL2SsB3g+xIoQK2B5\nBSLRLBEvRySSw7MXiFpZolaWGPN4LODyGOf4kCjZcizgVdsFYtV2vnpOZf/e7fny+NLYTZVzCjm8\nbB57IQ9zPswC7cCtvtL2HEvbmZo2GDOL4c+CPwP5WZjJwVQe7uZhilLcDbTTQPLECVKpVGguhlFu\nY2xuOzzCu3wUyOuNmdsNi/iNGzf429/+xvbt2/nVr36FZVm8+OKL9PT0cPr0aS5evEhnZye9vb2N\nHkpkTVFuy3rQsIjv2rWLP//5z6HHjh07tuILWjsMfR9U95dYNuW2YZTbofSNTRERg6mIi4gYTEW8\nLkPvraD7S0hDhl5X5XYoFXEREYOpiIuIGExFXETEYCriIiIGUxEXETGYiriIiMFUxEVEDKYiLiJi\nMBXxugy9t4LuLyENGXpdlduhVMRFRAymIi4iYjAV8boMvbeC7i8hDRl6XZXboVTERUQMpiIuImIw\nFXEREYOpiIuIGExFXETEYCriIiIGUxEXETGY22hALpfj+PHj5PN5CoUCTzzxBC+88ALT09OcOXOG\nsbExurq66O3tpa2trRVrFlkRym1ZDxoW8UgkwvHjx/E8j2KxyLFjx3j88cf54IMP2LNnD88//zxD\nQ0OcP3+ew4cPt2LNLWLovRV0f4llU24bRrkdallvp3ieB5ReuRQKBQBGRkY4cOAAAMlkkuHh4SYt\nUaR5lNtiuoavxAGKxSKvvvoqX3zxBc899xw7duwgk8kQj8cBiMfjZDKZpi5UpBmU22K6ZRVx27Z5\n4403mJ2d5eTJk9y6deueMZYV/r8no6OjjI6OVvdTqRTJZHflrECE7dvlbTvkuFVzvNI65TbYtxjJ\n5OaaccExTk2/tThtpXVChtZuh/QlE5SebScQbp3HcQMRNta1wLXBBcuxsR0bx41g20VcaxOOVcCh\ngEsehwKP8xAO38ShUI48bs22S74aTmC//napjVT2nSJutIjlFCHiwyYgloSvA7lA5MtttqY/GNnF\n1sqVws3B5gK4RWgrQgewUI5soAUYHBykIpFIkEgkqEe5HZjesNzeyxZcvhHI642Z28sq4hVtbW3s\n3r2bq1evEo/HmZycrLbt7e2h54T9EqXTn9DX9z73Jmow6SpZVQmrpt+uOV7JhgjhWVI5ZtPXNxEY\nFxwTLUdk6XzBpUXKEQ1sh+3X9tnQ9165z6tpg+M8IFYOL2SsB3g+xIoQK2B5BSLRLBEvRySSw7MX\niFpZolaWGPN4LODyGOf4kCjZcizgVdsFYtV2vnpOZf/e7fny+NLYTZVzCjm8bB57IQ9zPswC7cCt\nvtL2HEvbmZo2GDOL4c+CPwP5WZjJwVQe7uZhilLcDbTTQPLECVKpVGgu3o9yG+Ny2+ER3uWjQF5v\nzNxu+J743bt3mZ2dBSCbzXLt2jW2bt3K3r17SafTAKTTafbt29fooUTWFOW2rAcNX4lPTk7y9ttv\nUywW8X2fJ598ku9+97vs3LmT06dPc/HiRTo7O+nt7W3FekVWjHJb1oOGRXz79u387ne/u6d/y5Yt\nHDt2rCmLEmkF5basB/rGpoiIwVTERUQMpiIuIrKSaj9V2mQq4iIiK632Y/9NLOYq4nUZem8F3V9C\nGjL0upqU2y0q4KAiLiKy8vyaaCIV8bqa/Mw3y5dc9v2HG/ocSAOGXleTcru2iDdxOhVxEZGV5Afa\nFvz3UkVcRMRgKuIiIgZTERcRMZiKuIiIwVTERUQMpiIuImIwFXEREYOpiNelryYb+xxIA4ZeV+V2\nKBVxERGDqYiLiBhMRbwu3V/C2OdAGjD0uiq3Q6mIi4gYTEVcRMRgDf+1+4pischrr71GR0cHR48e\nZXp6mjNnzjA2NkZXVxe9vb20tbU1c60iK055LaZb9ivxv/zlL2zdurW6PzQ0xJ49e+jv7yeRSHD+\n/PmmLFCkmZTXYrplFfHx8XGuXLnCwYMHq30jIyMcOHAAgGQyyfDwcHNWKNIkymtZD5ZVxN955x1e\neuklLGvxA/KZTIZ4PA5APB4nk8k0Z4UiTaK8lvWgYRG/fPky7e3tdHd34/v1P5YT/EUQWeuU17Je\nNPzD5o0bNxgZGeHKlStks1nm5uZ46623iMfjTE5OVtv29vbQ80dHRxkdHa3up1Ipksnu8l7tPwdd\nu2+Xt+2Q41bN8UrrlNtg32Ikk5trxgXHODX91uK0ldYJGVq7HdKXTFB6tp1AuHUexw1E2FjXAtcG\nFyzHxnZsHDeCbRdxrU04VgGHAi55HAo8zkM4fBOHQjnyuDXbLvlqOIH9+tulNlLZd4q40SKWU4SI\nD5uAWBK+DuQCkS+32Zr+YGQXWytXCjcHmwvgFqGtCB3AQjmygRZgcHCQikQiQSKRoNZXzWtQbq92\nbu9lCy7fCOT1xsxty7/fy5Aa169f58KFCxw9epRz586xZcsWenp6GBoaYmZmhsOHDy/rcU6cSNPX\n9z73Jmow6SpZVQmrpt+uOV7JhgjhWVI6dvz4Y/T1TQTGBcdEyxFZurbg0iLliAa2w/Zr+o4fhr73\nyn1eTRs81wNi5fBCxnqAV4RYEWIFLK9AJJol4uWIRHJ4zgJRO0vUzhJjHo8F/pvHOMcNomTLsYBX\nbmMs4JUjxnz1HK96bL56zCvvx6rtPJsq52RzeAt5vIU8zpwPs0D7cfior7Q9x9J2pqYNxkwp/Bko\nzpba3CxM5WAqD3fzMEUp7gbaaeDE8tO5aqXyGpTbrcztn/EI7/J/gbzemLn9wJ8T7+np4dq1axw5\ncoR//OMf9PT0POhDiawZymsxzbI/Jw6we/dudu/eDcCWLVs4duxYUxYl0krKazGZvrFZl6H3VtD9\nJaQhQ6+rcjuUiriIiMFUxEVEDKYiLiJiMBVxERGDqYiLiBhMRVxExGAq4iIiBlMRFxExmIp4XYbe\nve5LLvv+ww19DqQBQ6+rcjuUiriIiMFUxEVEDKYiXpeh91bQ/SWkIUOvq3I7lIq4iIjBVMRFRAym\nIi4iYjAVcRERg6mIi4gYTEVcRMRgKuIiIgZTERcRMZiKeF2G3ltB95eQhgy9rsrtUCriIiIGUxEX\nETGYinhdht5bQfeXkIYMva7K7VCW7/vr56cREdlgVuWV+ODg4GpMu6pz62feGHSd1/+8qz13Lb2d\nIiJiMBVxERGDOSdOnDixGhN3dXWtxrSrOrd+5o1B13n9z7vacwfpD5siIgbT2ykiIgZTERcRMZjb\n6gmvXr3KH//4R3zf55lnnqGnp6dpcw0MDHD58mXa29s5efIkANPT05w5c4axsTG6urro7e2lra1t\nRecdHx/n7NmzZDIZLMvi4MGDHDp0qOlz53I5jh8/Tj6fp1Ao8MQTT/DCCy+05GcGKBaLvPbaa3R0\ndHD06NGWzbtWtCq3N1peg3L7vvwWKhQK/k9/+lP/9u3bfi6X83/5y1/6n332WdPm++c//+l//PHH\n/i9+8Ytq35/+9Cd/aGjI933fP3/+vH/u3LkVn/c///mP//HHH/u+7/tzc3P+z3/+c/+zzz5rydzz\n8/O+75ee61//+tf+hx9+2JJ5fd/3L1y44Pf39/u//e1vfd9vzXO9VrQytzdiXvu+cruelr6dcvPm\nTR599FE6OztxXZennnqK4eHhps23a9cuNm/evKRvZGSEAwcOAJBMJpsyfzwep7u7G4BYLMbWrVsZ\nHx9vydye5wGlVy6FQgFozc88Pj7OlStXOHjwYLWvFfOuFa3M7Y2Y16Dcrqelb6dMTEzw8MMPV/c7\nOjq4efNmK5dAJpMhHo8DpaTMZDJNne/27dt8+umn7Ny5syVzF4tFXn31Vb744guee+45duzY0ZJ5\n33nnHV566SVmZ2erfa1+rlfTauf2es9rUG7Xs+H/sGlZzbuv8Pz8PG+++SYvv/wysVisJXPbts0b\nb7zBwMAAN2/e5NatW02ft/L+bHd3N/59PrHazOdallpveQ3K7Xpa+kq8o6ODO3fuVPcnJibo6Oho\n5RKIx+NMTk5W2/b29qbMUygUOHXqFE8//TT79+9v6dwAbW1t7N69m6tXrzZ93hs3bjAyMsKVK1fI\nZrPMzc3x1ltvtfTnXW2rndsbJa9BuV2rpa/Ed+zYweeff87Y2Bj5fJ5Lly6xb9++ps7p+/6S/4Lu\n3buXdDoNQDqdbtr8AwMDbNu2jUOHDrVs7rt371b/ly+bzXLt2jW2bt3a9Hm///3vMzAwwNmzZ3nl\nlVf41re+xc9+9rOWPddrQatzeyPlNSi376fl39i8evUqf/jDH/B9n2effbapHzHs7+/n+vXrTE1N\n0d7eTiqVYv/+/Zw+fZo7d+7Q2dlJb2/vPX8k+qpu3LjB8ePH2b59O5ZlYVkWL774Ijt27Gjq3P/6\n1794++23KRaL+L7Pk08+yfe+9z2mp6eb/jNXXL9+nQsXLlQ/htWqedeCVuX2RstrUG7fj752LyJi\nsA3/h00REZOpiIuIGExFXETEYCriIiIGUxEXETGYiriIiMFUxEVEDKYiLiJisP8HMo92mK9nGnYA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc9b3cc43c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "analytic_fisher_info = analytic_model.fisher_information(modelparams, expparams)[0,0,...]\n",
    "numerical_fisher_info = numerical_model.fisher_information(modelparams, expparams)[0,0,...]\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(analytic_fisher_info)\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(numerical_fisher_info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Randomized Benchmarking Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test that we get multiparameter Fisher information calculations correct as well, we compare to the zeroth-order non-interlaced randomized benchmarking model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class NumericalRandomizedBenchmarkingModel(ScoreMixin, RandomizedBenchmarkingModel):\n",
    "    pass\n",
    "\n",
    "analytic_model = RandomizedBenchmarkingModel()\n",
    "numerical_model = NumericalRandomizedBenchmarkingModel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now make experiment and parameters to test with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "expparams = np.empty((150,), dtype=analytic_model.expparams_dtype)\n",
    "expparams['m'] = np.arange(1, 151)\n",
    "\n",
    "modelparams = np.empty((500, 3))\n",
    "modelparams[:, 0] = np.linspace(0.1, 0.999, 500)\n",
    "modelparams[:, 1] = 0.5\n",
    "modelparams[:, 2] = 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make sure that the returned Fisher information has the right shape. Note that the Fisher information is a four-index tensor here, with the two indices for the information matrix itself, plus two indices that vary over the input model parameters and experiment parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "afi = analytic_model.fisher_information(modelparams, expparams)\n",
    "assert afi.shape == (3, 3, modelparams.shape[0], expparams.shape[0])\n",
    "nfi = numerical_model.fisher_information(modelparams, expparams)\n",
    "assert nfi.shape == (3, 3, modelparams.shape[0], expparams.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We check that each Fisher information matrix has errors that are small compared to the analytic FI alone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6626103445273564e-07"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(afi - nfi) / np.linalg.norm(afi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we plot the trace-inverse of each to check that we get the same Cramer-Rao bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def tr_inv(arr):\n",
    "    try:\n",
    "        return np.trace(np.linalg.inv(arr.reshape(3, 3)))\n",
    "    except LinAlgError:\n",
    "        return float('inf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def crb(fi):\n",
    "    return np.apply_along_axis(tr_inv, 0, np.sum(fi.reshape((9, modelparams.shape[0], expparams.shape[0])), axis=-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5gAAAGICAYAAAAkgU8PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl0VEXaBvDnrWyQFRKC7LIEjYQtGGSLQFgFRINKRAXX\nQUHQKIuIQQQFFQQV9QNcUFQYJTMsKiqKKC5BRhQR6LhFXFFkGyXsSaq+PyIxMECS7tt9+3Y/v3Pm\nnEknfe9zMs28qVtVb4kxxoCIiIiIiIjIQ8ruAERERERERBQYOMAkIiIiIiIiS3CASURERERERJbg\nAJOIiIiIiIgswQEmERERERERWYIDTCIiIiIiIrIEB5hERERERERkCQ4wiYiIiIiIyBJBMcA8cuQI\nJk6ciI0bN9odhYiIyG+wPhIRkdWCYoD5yiuvoFOnTpX+eZfL5cU01mFO6zklK3Nayyk5AedkZU5n\nCNT6CDgnK3Nayyk5AedkZU7rOSWruzkdN8CcN28ehg8fjnHjxh33+qZNm3DbbbchOzsbK1asKHt9\n8+bNaNCgAWJjYyt9j0D/H93XnJITcE5W5rSWU3ICzsnKnL7H+ng8p2RlTms5JSfgnKzMaT2nZA2a\nAWZGRgZycnKOe01rjQULFiAnJwezZ89GXl4etm/fDgDIz8/Ht99+i7y8PKxZs8aOyERERF7H+khE\nRP4g1O4AVZWcnIxdu3Yd91pBQQHq1q2LxMREAECXLl2wYcMG1K9fH0OGDAEAvP/++4iJifF5XiIi\nIl9gfSQiIn8gxhhjd4iq2rVrF2bMmIFZs2YBANavX48vvvgCN910EwDggw8+QEFBAa6//vpKXc/l\nch03BZyVlWV9aCIi8ku5ubll/z0lJQUpKSk2pvEM6yMREVnJnRrpuBlMbzjZL+vXX3+1KU3lxcTE\noLCw0O4YFXJKTsA5WZnTWk7JCTgnq1Ny1qtXj4Om03BqfQSc8xlkTms5JSfgnKzMaT2nZHW3Rjpu\nD+bJxMfHY/fu3WVf7927F/Hx8TYmIiIish/rIxER+ZojB5jGGJRf2ZuUlIQdO3Zg165dKC4uRl5e\nHtLS0qp8XZfLddw0MBERBb7c3FzHdPSrCOsjERFZyZ0a6bg9mHPmzEF+fj4KCwsRFxeHrKwsZGRk\n4PPPP8fChQthjEGPHj2QmZnp0X2csATIKdPrTskJOCcrc1rLKTkB52R1Ss569erZHcEyrI/Hc8pn\nkDmt5ZScgHOyMqf1nJLV3RrpuD2Y2dnZJ309NTUVqampHl37WDMD7schIgoeubm5jm/uA7A+EhGR\n9dypkY4bYHpTIPyBQUREVcNBU8VYH4mIglPQNvkhIiIiIiIi+3GASURERERERJbgALMcdskjIgo+\ngdRF1ltYH4mIgpM7NZJ7MMvhHhMiouDDPZgVY30kIgpO3INJREREREREtuEAsxwuASIiCj5cIlsx\n1kciouDEJbIe4hIgIqLgwyWyFWN9JCIKTlwiS0RERERERLbhAJOIiIiIiIgswQEmERERERERWYID\nzHLYxICIKPiwyU/FWB+JiIITm/x4iE0MiIiCD5v8VIz1kYgoOLHJDxEREREREdmGA0wiIiIiIiKy\nBAeYREREREREZAkOMMthEwMiouDDJj8VY30kIgpObPLjITYxICIKPmzyUzHWRyKi4MQmP0RERERE\nRGQbDjCJiIiIiIjIEhxgEhERERERkSU4wCQiIiIiIiJLcIBJREREREREluAAsxy2YSciCj48pqRi\nrI9ERMGJx5R4iG3YiYiCD48pqRjrIxFRcOIxJURERERERGQbDjCJiIio0kxJid0RiIjIj3GAeQrG\nGLsjEBER+R2zcZ3dEYiIyI9xgHkqX222OwEREZHfMauW8SEsERGdEgeYp6DfWmZ3BCIiIv9TdJQP\nYYmI6JQ4wDyV7T/C/PK93SmIiIj8ivQdBL2KD2GJiOjkOMA8BekxEOat5XbHICIi8ityXjfg1x9h\nftpmdxQiIvJDHGCegnTrC7PlM5g9O+2OQkRE5DckLAzScyAMt5IQEdFJcIBZjsvlQm5uLgBAIqMh\n6b1YQImIAlxubi5cLpfdMfxa+foIANL1Apj8z2F2/mpjKiIi8jZ3amSol7I4UkpKClJSUsq+lj6Z\n0HePguk3GFIzwcZkRETkLVlZWXZH8Hv/Ux8joyDdB8C88S/Itdk2JiMiIm9yp0ZyBvM0JLYmpEtP\nzmISERGdQHpdBLPpE5hdO+yOQkREfoQDzApIn0EwH78H88deu6MQERH5DYmKhnTrB/Pmv+2OQkRE\nfoQDzApIjXhIpwyYt9lRloiIqDzpfRHMxo/ZEI+IiMpwgFkJ0vcSmLw1MPv+sDsKERGR35DoWEjX\nPpzFJCKiMhxgVoLUTIB06Arz5lK7oxAREfkV6Z0J82ke92ISEREADjArTfpnwXz8LsyeXXZHISIi\n8hsSEwfJGADz2kt2RyEiIj/AAWYlSY14SLcLWECJiIhOIH0yYbZuhNn+o91RiIjIZhxgVoH0HQSz\neQPMbz/bHYWIiMhvSPVISL/LoJe/aHcUIiKyGQeYVSCR0ZC+g6BXLLI7ChERkV+R7v2An7+HKfjS\n7ihERGSjULsDeNv27dvxxhtvoLCwEC1btkSfPn08up5kDIB55zWY77+BNDnLopRERES+ZXl9DAuH\nXHQF9LLnocY/ABGxKCkRETlJwM9g1q9fH8OHD8ftt9+Ob775xuPrSXgEZOAQ6KXPwxhjQUIiIiLf\ns7o+AoB0zAD2FwJbN1pyPSIich7HzWDOmzcPGzduRFxcHGbNmlX2+qZNm7Bw4UIYY5CRkYHMzMyy\n73366adYvXo1unbtakkG6dIL5u0VgGsj0PJcS65JRETkCb+ojyEhUJlDoZcuhEppC1EhllyXiIic\nw3EzmBkZGcjJyTnuNa01FixYgJycHMyePRt5eXnYvn172ffT0tIwceJEfPjhh5ZkkJAQqMuuhc59\nFqa42JJrEhERecIf6iMAILUjEBkF89Fq665JRESO4bgZzOTkZOzadfxZlAUFBahbty4SExMBAF26\ndMGGDRtQv3595Ofn4z//+Q+KiorQrl0764K0OQ94dyXM+6sgPS+07rpERERu8Jf6KCJQQ4ZDz5kK\nk3Y+JDLKsmsTEZH/c9wA82T27t2LhISEsq/j4+NRUFAAAGjRogVatGhh+T1FBOryf0DPngTToSsk\nOtbyexAREXnCjvoIANKoGaR1e5jXl0AGX++VexARkX8KiAGmp1wuF1wuV9nXWVlZiImJqfiNyS1x\nsFN34K2liLz2Vu8FPIXw8PDK5bSZU3ICzsnKnNZySk7AOVmdkhMAcnNzy/57SkoKUlJSbEzjX9yu\njwD00BEoHH8dIvtdipC6DbwV8ZSc8hlkTms5JSfgnKzMaT0nZXWnRgbEADM+Ph67d+8u+3rv3r2I\nj4+v9PtP9ssqLCys1HvNBYOhJ9+M4k49IfUaVfqeVoiJial0Tjs5JSfgnKzMaS2n5ASck9VJObOy\nsuyO4TV21keEhAG9B6Fw4eMIGT2p0ve0ipM+g8xpHafkBJyTlTmt55Ss7tZIxzX5AQBjzHFHhCQl\nJWHHjh3YtWsXiouLkZeXh7S0tCpf1+VyHTdKrwyJiYUMGAydu4DHlhAROVBubu5xs3RO5k/1EQCk\n10Bg+48w+Zuq/F4iIrKfOzVSjMNGRXPmzEF+fj4KCwsRFxeHrKwsZGRk4PPPPy9rw96jR4/j2rC7\n49dff630z5riYuh7s6Eyh0LadfLovlXhpKcfTsgJOCcrc1rLKTkB52R1Ss569erZHcEy/lgfAcB8\nvh562QtQ98yBhIZ5dO+qcMpnkDmt5ZScgHOyMqf1nJLV3RrpuCWy2dnZJ309NTUVqampPk5TSkJD\noa4aAf3sI1At2kKqVbclBxERBS9/rI8AgLYdgA/fhnlrOWRA4C5HJiKiUo5cIust7i4BAgA5uxXk\n7FYwr71kcSoiIvKmQFoi6y0e1UcRqCtuhFn9CsyuHRYnIyIib3KnRjpuBtObPO0eKJddBz3lFphO\nGZAGTSxMRkRE3hLITX6s4nF9TKwD6ZMJ/dJTULfcDRGxMB0REXlL0DT58VcSWwNy8VXQi+bBaG13\nHCIiIr8hfTKB3b8Dn39sdxQiIvIiDjDL8WQJ0DFyfh/AGJi8dyxKRURE3sQlshWzpD6GhkFdNRL6\n5WdgDh+0KBkREXkTl8h6yIoDtkUpqKE3Qz8yGaZVGqRG5c8bIyIi3+MS2YpZUR8BQM5uCTmnDczy\nRZArbrQgGREReROXyPoJadgEcn5f6MXzeTYmERFROZJ1PczGdTDf5tsdhYiIvIADTC+RCy8Hft8O\n82me3VGIiIj8hkTFQF05AnrhYzBHj9gdh4iILMYBZjlW7DE5RsLCoK69Feblp2AK/7TkmkREZD3u\nwayYlfURACS1I+TMZjCv/tOyaxIRkfW4B9NDVu0xOUaang3p2B3mpacgN4637LpERGQd7sGsmNX1\nEQDkihtLj/Y6twukyVmWXpuIiKzBPZh+SC6+CubH72A+X293FCIiIr8hMXGQy/8B/dwcmKIiu+MQ\nEZFFOMD0MgmPgLr2VujF82D2/WF3HCIiIr8h7c8H6jaAWbHI7ihERGQRDjB9QJq3gHTuAf3CE+wq\nS0RE9BcRgRo6CuaT92G+3mJ3HCIisgAHmOVY3cSgPLnoSuC/u2E+fNsr1yciIvewyU/FvFofY2Kh\nrr4F+tlHYQ7u98o9iIjIPWzy4yFvNDE4RkLDoP4xFnrmRJizW0HOqOeV+xARUdWwyU/FvFkfAUBa\nnQtp3R5m8ZOQ4WO9dh8iIqoaNvnxc1K3IeTCIdALHoYpKbE7DhERkd+Qy66D+akA+pMP7I5CREQe\n4ADTxySjPxAZBfPaS3ZHISIi8hsSEQF1wxiYl5+G2f273XGIiMhNHGD6mCgFdf1tMHnvwORvsjsO\nERGR35DGzSF9L4F+6iGYYh5dQkTkRBxgluPNJgblSWxNqOtvL21o8Od/vX4/IiI6NTb5qZiv6iMA\nSO+LgZg4mOUv+uR+RER0au7USA4wy0lJSfFZswc5pw2kax/oZ2bDaO7HJCLyNfP9NwBKGxh4s4FN\nIPBpfVQK6rpsmE/zYL74xCf3JCKiv5kjh2G2fgbAvRrJAaaN5MLLAWNgVvrmqTAREf1Nr1pqdwQ6\nBYmOhRo+Dvr5x2H27LI7DhFRcNn3B/SieW6/nQNMG4kKgfrHWJgP3oL5arPdcYiIgos2dieg05Ck\ncyC9M6GffgimuNjuOEREwcNoQLk/TOQA02ZSIx7q+uzSo0v+2Gt3HCKi4GG03QmoAtJ3EFA9CmbZ\n83ZHISIKHloDwgGmo0mLVEi3ftDzH2TXPCIiX9EcYPo7UQrqhtthPl/P8zGJiHzFGM5gBgLpPxiI\nqQHz0tN2RyEiCg6cwXQEiY6FuvkumJeegvn5e7vjEBEFPq0BEbffzgFmOb5sw36isvMxv9kK/cFb\ntmQgIgoqf+3B5DElFbOzPgKANGwCGTIceu79MPv32ZaDiCgolNuD6U6NDPVGJqdKSUmxtVW9VI+E\nGnUX9MyJMPXPhDRLti0LEVHA+2sG01fHbziZ3fURAFSHbtA/fQf99Cyo7HsgKsTWPEREAavcHkx3\naiRnMP2M1GkAdc0t0PNnsOkPEZE3cQ+m48gl15Qe77XsRbujEBEFLs0usgFH2pwH6dq3dCnQ0SN2\nxyEiCkzcg+k4EhICNXw8zGd50OvfszsOEVFgYpOfwCQXXg6pdQbMwsdg+JSdiMh6/P9WR5KYWKjR\nd8PkPgvzbb7dcYiIAg+b/AQmEYFclw2zdxfMq/+0Ow4RUeAxxu4E5Cap3wjq+ttLj/fa+ZvdcYiI\nAovhEtmAJWHhpa3Z//M+9Lp37Y5DRBRYOIPpaNKyHWTgEOjH74U5sN/uOEREgaNckx93cIDp5yS2\nBtQtd8P8+zmYb7baHYeIKHBwgOl4qnt/SMtzoec9AFNcZHccIqLAwD2YgU/qNYL6xxjoJ2fC7Nhu\ndxwiosDAJj8BQQZfB0RUg1k0F4bLnomIPMc9mMFBWqRCBg2DnjOFx5cQEVlBczASCESFQA0fB/PL\njzCvLLY7DhGR83EPpnVcLhdyc3PtjnFKKr03JL039JypMAcP2B2HiMjZ/prBzM3NhcvlsjmMf/P3\n+ijVqkPdOhlmw0fQ771hdxwiImcrtwfTnRoZ6o1MTpWSkoKUlBS7Y5yW9B8M/Lm39IzMSbPsjkNE\n5Fx/7cHMysqyOYj/c0R9jK0BddsU6Jl3wsTGAd0vsDsSEZEzlduD6U6N5Aymw4gIZMhwIDoGB//v\nfhhdYnckIiJn4h7MgCOJdaBuuRt68XwUuTbZHYeIyJm4BzP4iAqBumEMzL4/YV56mk0NiIjcwT2Y\nAUkaNYMaPg4H50yF+ek7u+MQETkP92AGJwkLR9S4+2C+/wZm2QscZBIRVRVnMAOWnNMG1W+4Dfqx\ne2F+/cnuOEREzsJzMIOXREZD3TYFZsunMK8vsTsOEZGz8BzMgBbeoRvksmuhH7kHZuevdschInIO\nnoMZ3CQ6FmrMvTDr34d+e7ndcYiInIMzmAFPdcyADLwc+uHJMHt22R2HiMgRjNYQ7sEMbhJbE2rM\nfTDvvQG9lu3ZiYgqhXswg4LqegGk10DohyfxHGkiosrQ3INJACS+Vukg841/Q69bY3ccIiL/xxnM\noKF6XQzp3BP64bthCv+0Ow4RkX8z3INJf5HEOlC33wuz7EXoj9+zOw4RkX/jHsygogZkQVI7cZBJ\nRFQRzmBSeVK3AdTY+2CWPQ+dx5lMIqJT4gxm0JHMqyCtz4OelQOz7792xyEi8k8eNvkJtTCKX9qw\nYQM2btyIQ4cOoUePHmjdurXdkbxO6jaEGjsN+uHJ0LoE6vw+dkciIvI/QT6DGZT1UQTIvAoICYGe\nNQlqzH2QGvF2xyIi8i9aAx40+Qn4AWb79u3Rvn17HDhwAC+++GJQFFAAkDoNoMZNg559N3RJMVT3\n/nZHIiLyL0F+fnDQ1kcRyEVXQCsFPSsHauw0SM0Eu2MREfkP49kSWccNMOfNm4eNGzciLi4Os2bN\nKnt906ZNWLhwIYwxyMjIQGZm5nHvW7p0KS644AJfx7WV1K4HNW469OxJ0CUaqueFdkciIvIfATaD\nyfpYNerCy6FDQqFn3VU6yIxPtDsSEZF/0EHW5CcjIwM5OTnHvaa1xoIFC5CTk4PZs2cjLy8P27dv\nL/v+4sWL0a5dOzRu3NjHae0niXWgxk2HeecV6LeW2R2HiMh/BNgeTNbHqlP9LoV0uwD6obtgdv5m\ndxwiIv/g4R5Mxw0wk5OTERUVddxrBQUFqFu3LhITExEaGoouXbpgw4YNAIA333wTW7duxfr16/HO\nO+/YEdl2UusMqDsehMlbA73seZggXxZGRAQg4GYwWR/do/oMgvS9BPqhiTC//GB3HCIi+3EPJrB3\n714kJPy9fyI+Ph4FBQUAgH79+qFfv352RfMbUjMBavwD0I9NBQ7MBa4aAVEhdsciIrJPEDxsY32s\nHNW9H3RkFPTDd0ONyoE0S7Y7EhGRfYJtD6Y3uFwuuFyusq+zsrIQExNjY6LKCQ8Pr1rOmBiYex7F\ngVmTIM89isjROZDQMO8F/EuVc9rIKVmZ01pOyQk4J6sTcv7x1wxmbm5u2WspKSlISUmxK5LfcWp9\nBNz4DPYcgKL4Wjg4935UH3UXwtq09164cpzwbwVgTm9wSlbmtJ6/Zz0cFg4TEQHAvRoZEAPM+Ph4\n7N69u+zrvXv3Ij6+8m3HT/bLKiwstCyft8TExLiV04zKgX56Fv68fwLUzRMhEdW8kO5v7ua0g1Oy\nMqe1nJITcE5Wf89pjCmbwczKyrI5jfcEa30E3PwMJrWAjLwTB56YDnXlTZC0dO+EK8ff/60cw5zW\nc0pW5rSev2fVhw8BxSUA3KuRjtuDCZT+YVB+H2FSUhJ27NiBXbt2obi4GHl5eUhLS6vydV0u13Gj\n9EAlYeFQN02A1IyHnj0JpvBPuyMREfmW+btDXm5u7nGzdE7G+ug5SWoBNeZe6CXPQL/3ht1xiIh8\nr9weTHdqpBiHdXyZM2cO8vPzUVhYiLi4OGRlZSEjIwOff/55WRv2Hj16/E8b9qr69ddfLUrsPZ4+\n/TDGwKxYDPPph1C33gM5o56F6f7m709pynNKVua0llNyAs7J6u85TXEx9OjBaPjqf+yOYhnWx+N5\nXCN37YB+bCqkTQfIJVdDPNiPdDr+/m/lGOa0nlOyMqf1/D2rfj0XOHIYDUbf6db73V4iu23bNmze\nvBk//PADdu/ejcLCQuzfvx8AEBsbi+joaCQmJqJp06Zo06YNGjZs6O6tjpOdnX3S11NTU5GammrJ\nPYKFiEAGDYVOqAX90ESokRPZ2ICIgoPx7IyvithRI1kfrSWJdaAmzID+v+nAM7OB67IhYeF2xyIi\n8j4Pa2SVBphFRUVYtWoV1q1bh9q1a6NJkybo2bMnEhISEBUVhaioKBhjcPDgQRw6dAg7duzAtm3b\nsGTJEvzxxx/o3r07evToAeWlp4CeOtbMIJD345yM6noBTM1a0E9Mgxo2CtKuk92RiIi8SxtA/b38\nx4rmPoFcI4O1Pkp0LNSY+2AWPAL96D1QN98FifLfxhxERJbQf5+D6U6NrPQS2a1bt+KNN95A165d\nkZaWhtDQqk1+Hj16FOvWrcOGDRtw2WWXoUmTJlV6v685YQmQ1dPr5scC6CemQS64FKrnQMuu6+/L\nAMpzSlbmtJZTcgLOyervOc3hg9DjrkPDZR9acr1gqpFOqI+AtZ9BozXM0oUwmz+Fyr4HUusMS64L\n+P+/lWOY03pOycqc1vP3rPqVxYAoNBgx1q33V6oCvvvuuygpKcH48eMhbh66GR4eju7duyM9PR2v\nvvoqdu3ahfPOO8+ta5F3yJlJpcuB5kyF3v07ZPB1PCuTiAJTuRlMT7FGBj5RCjL4euj42tAzJkDd\nnANp0tzuWERE3qENEOb+apoK33nw4EE0bNgQvXv3drtwlhcaGopLLrkEEX+dreJPgqlL3qlIrTOg\n7pwB88sP0E9Mhzl4wO5IRETWs6iLbLDUSNbHUqrnhVBXjYB+bCr0Jx/YHYeIyDtMiUddZCucwYyM\njETz5tY/pWvTpo3l1/QUD9guJVExUNlTYJY8Df3gHVCjJ0Fq17U7FhGRdbQum8H0ZF9hsNRI1se/\nSduOUAlnQP/fdOhff4JcdKXXOswSEdlCa+CvVYy2nYP5888/W3EZ8iMSGgp11UhIxgDoGRNgvtps\ndyQiIut4uYtseayRgUcaNoG6axbM11ugn5wBc+Sw3ZGIiKxj/m7y444qH1OyfPlylJSUlLu/wbZt\n2zBhwgS3Q5D/Uhn9YerUh37qIcjFV0F1u8DuSEREntPao+J5KqyRwUNia0CNmQazeC70gxNKV/sk\nJNodi4jIc1qXLZF1R5UHmPXq1UPTpk2Pe61atWpuB/AnwdqGvSJyTpvS5j9PTIPe/gMk6x+QKnZI\nJCLyK9octwfTqiWggVojWR9PTsLCgGtuBVa/Av3AeKgRd0CSWtgdi4jIM8azY0qqPEpo1aoVIiMj\nj3utd+/eVb2MX+Iek1OTM+pBTZwJ/czDMLMnQd10B6RGvN2xiIjcY6zZg3miQK2RrI+nJiKQPpkw\ndRtCz30AMuBySI8BljR9IiKyhS4pewjrkz2Y5Qtnfn4+gMB4OksVk8jo0iVALdpCTx8LU5BvdyQi\nIvdo7+zBZI0MXtLqXKg7Z8J8tBrm2UdgjhyxOxIRkXs83IPpUXX15wNCyTtEKaiBQ6CuHgU99wHo\nNa/BGGN3LCKiqjHe2YNZHmtk8JHadaHunAkA0A+Oh9n5m82JiIjc4OEeTPbVLofnfFWetEqDmvgQ\nzEfvwDzzMDvoEZGznLAH091zMIMF62PlSUQE5PrbIV37Qj94B8zmDXZHIiKqmhP2YFp+DmYw4R6T\nqpHEOqXLgRbNLW1ucPNESO16dsciIqqYl/ZgBirWx6oREUjGAJiGTaGfnAk5vzfkwiE8L5OInEGX\nlM1g+vwcTC6NpNIntbdBul0A/eAEmM/y7I5ERFSxcjOY3sIaSZJ0DtSkh0vPy5wzBWbff+2ORERU\nMW0AFeL22z2qri1asBU3lT6pVRkDoG6dDL30eejF82CKjtodi4jo1EyJ1/dgskYSAEhcTagx0yBN\nzoK+73aYL7+wOxIR0ekZG/dgxsbGevJ2CjDSuDnUpEeAwn3Q94+H2bHd7khERCfngxlM1kg6RkJC\noDKHQl2XDb3gEehX/gmjS+yORUR0ctrGLrJEJ5LIKMhNd5QumZ0xAXr9WrsjERH9r3J7MIl8RVqk\nli6ZLciHfngyzB977I5ERPS/7JzBDDTskmcNEYHq3g9qzH0wK5fg4PyZPA+MiPwLu8hWCeujdaRG\nPNTtUyFnt4KeNgZFmz6xOxIR0fG0LtuD6dMustu2bcPmzZvx/fffY8+ePThw4AD2798PAIiOjkZ0\ndDQSExPRtGlTtG7dGo0aNXL3Vj7DLnnWkoZNoCbNhlmyAHra7VDDx0IaNbM7FhFRaYe8v5b/eKOL\nbKDVSNZHa4kKgQwcAnNWSxx86iHg3C6QQcMgYeF2RyMigjEayoNO61UaYB49ehRvvvkm1q9fj9q1\na6NJkybo1asXEhISEBUVhaioKGitcejQIRw6dAg7duzAtm3bkJubi//+97/o1q0bevXqBcU23UFD\nqkUiavRd+HP1a9CP3APpOwjSJxPiQWcqIiKPGc/2l5wMayRVlZzdEtEPPoN98x6EmT4Wavg4SP0z\n7Y5FRMHOwz4FlR5gbty4Ee+99x66deuG6dOnn7YAhoeHIy4uDnXq1EHbtm0BlBbe9evXY9asWcjM\nzMRZZ53ldmhyHtWxO0zSOdDPPgKzdSPUdbdBEhLtjkVEwUprS5v8sEaSu1RsHNTIiTB570DPyoEM\nyIL0uJBnZhKRfTzsU1CpAeaqVatQvXp1jB071u0bhYeHo2vXrujSpQtee+017Ny5E+np6W5fj5xH\nap0BNW6CjEtpAAAgAElEQVQ6zKpl0NPHQIYMhzqvq92xiCgYWdjkhzWSPCUikPTeMGelQD/zMMyW\nz6CuuxVSI8HuaEQUjLQGxIvnYB44cADJycno1q2b2zcpLyQkBJmZmYiLi7PkeuQsokKg+g+GunUy\nzKsvQS94GObgAbtjEVGwsWgGkzWSrCS160Hd8SCk2dmlZ2Zu/NjuSEQUjLw9g3ls38ipHD58GF9/\n/TV+++03HDp0CBEREahRowaSk5MRHx9/yve1atXKvcQUEKRxc6i7H4H517PQ92ZDXZcNOZufCSLy\nEaMt2YPJGklWk9BQyEVXwqS0g17wMGTzBkjWDZDIU3/OiIgs5eFDWLe7yP7yyy9YtWoViouLceaZ\nZ6JmzZqoX78+jh49iv3792PlypU4ePAgWrdujc6dO7sd0JdcLhdcLpdXOgrS/5KIapChN8Ns3gD9\nzGxIu86QS66GRFSzOxoRBboTjimxuktqoNVI1kffk2bJUJMfhfnXc9BTb4W6ZjSkRardsYgoGJRr\nhOdOjXRrgLlu3TocOXIE11xzDcLCwk77swUFBVixYgX69++P8HD/br/NNuz2kNbtoaY8DvPy06VF\n9NpsyFn834GIvKjc8h+rB02BWCNZH+0h1SIhw0bBbN0I/fzjkJbnQgZfB6kWaXc0IgpkWgPio2NK\njjnrrLNQq1atSv1sUlISmjZtin379vl18SR7SVQM5IYxMJvWQz/1EKR9OiRzGCQiwu5oRBSIPGzB\nfjqskWQ1adkO6p7HYXIXQE+5FeqaWyDntLE7FhEFKqMBD44UdKu6VlQ4d+7cefxNlEKNGjXcuRUF\nGWnbEWrKY8C+P6HvzYYpyLc7EhEFIlNi+TmYx7BGkjdIZBTUtbdCXTUS+rk50IvnwRw+ZHcsIgpE\n5WYw3WFZdd2xYwc+++wzFBYWlp3nReQOiY6FGj4W6tKroefPgM5dAHP0iN2xiCiQeHEG82RYI8kq\n0urc0gexR49CT7kF5qvNdkciokBTbg+mO9xu8nOif/3rX6hWrRoWLVqEunXrIj4+Hh07drTq8hSE\npF1nqOYpMP98EnpqNtTVoyFnt7Q7FhEFAgvPwawM1kiykkRGQ67LLm2St+ARSOv2kEuvYadZIrKG\nhzOYlg0w09PTkZpa2t3sm2++gXgQiugYiYmD3HRH6d7MBQ9DWraDXHotJCra7mhE5GBGG4gPZzBZ\nI8kbpHV7qKmPwyx9Hvqe0VBX3Ahp18nuWETkdNqGPZgnExoait9//x1AaYOD5s2bW3Vpor/2Zj4O\nhIRATxkN89k6GGPsjkVETqW9twfzZFgjyVskMhpq2Cio4WOhl72Akrn3w/yxx+5YRORkHq7ysWwG\nc+3atfj2229xxhlnoE2bNmjXrh3q1atn1eWJIJFRkKtGwpzXDfqFJ4D1a6GuvAlSM8HuaETkNMa3\nezBZI8nb5KyWUPfMgXk9F3pqNuTiqyBd+0J8+CCFiAKEh30KLPt/nfbt2+Oxxx7DDTfcgPDwcKxc\nudKqSxMdR5q3gJo8B9KgMfS92dBr34TR2u5YROQk2rd7MFkjyRckLBwqcyjU2Gkw69ZAP3QXzG+/\n2B2LiJzGX2Yw1V9PyOrUqYM6depYdVmfcrlccLlclh+6TdaTsDDIxVfCpKVDv/A4zH/ehxp6M6R+\nI7ujEZETlJvBzM3NRUpKClJSUrx2O6fXSNZHZ5EGjaHunAHz3pvQMydAegyE9LsUEhpmdzQicgKt\nASndg+lOjbRsgPnll19i1apVSE9PR5s2bZCQ4Lxli97+A4OsJ/UbQU14EOb9t6Bn3QXp0gsycAgk\noprd0YjIn+mSsqezvhg0Ob1Gsj46j6gQSM8LYdp2gP7nfJip2aXbSs5pY3c0IvJ35WYw3amRli2R\nbdiwIYYOHYp9+/Zh7ty5mDJlilWXJjotUSFQGf2hpj4O/LkXevIomM/XswkQEZ2aMR51yKsq1kiy\niyQkQo2eVHq29MLHoJ+eBfPHXrtjEZE/83APpmUzmM2bN8eePXuQmZmJzMxM/nFPPiexNSE3jIH5\najP04vnAh2+XtmxPdN5yNCLyMg/P+Koq1kiyk4gAbTtCndMW5vUl0FNvgVw4BNK9PyTEdw9aiMgh\nPNyDaekMZtu2bcu+5hlfZBdJbg11zxxI0jnQ94+Ffj0XpqjI7lhE5E+M9ukxJayR5A8kohrUJddA\njX8A5vP10NPHwHz3ld2xiMjfaO0fXWSJ/ImEhkH1Hwx112yYbV9D33srzJdf2B2LiPyFNj6dwSTy\nJ1KvEdTYaZC+l0DPexD6hSdg9u+zOxYR+QsPH8J6ZYBpjMHmzZuxdetWHD161Bu3IKoUSazz196T\na0r3njw5E2bPLrtjEZHddIlP92CWxxpJ/kBEoDp0g7r3CSAsHHryKOgP3+axX0QElGiPaqRXBpi/\n//47nnvuObRo0QKffvopDh8+7I3bEFWKiEDadoS6dy5QpwH0fbdBr3wZpoh/2BEFrZJiIMSyNgRV\nwhpJ/kQio6GuuBEqewrMh29DP3gHzLav7Y5FRHYqKQY82J/tlQFmnTp18Mgjj0Aphc6dO6NaNR4Z\nQfaTiAioi6+EypkN8/P37DZLFMxKSjwqnp5gjSR/JGc2g7pzJqR7P+i5D0A/+yi7zRIFq5ISjx7C\nWvb4du3atVi+fDliYmIwcOBAdOjQwapLE1lKEusgZOREmPxN0C8/Dax9A2rIcEjdhnZHIyJf8fEM\nJmskOYEoBencEya1E8zruaXdZi+4FNJzoN3RiMiX/GUGs7i4GDNnzsTQoUPxxRdfYM2aNVZdmsgr\npEVbqMlzIK3SoGdOhF6yAObgfrtjEZEv+HgGkzWSnESqR0Jddi3UhJkwX2+FvucWFH2+3u5YROQr\nHs5gWjbAjIuLQ0REBJKTk3HjjTf6zbLDnTt3Yv78+Xj44YftjkJ+SEJDoXpdBDX1CeDwQewbcw30\nR6vZ5IAo0Pl4BtMfayTrI1VE6tRHyK2ToYYMx6EX5qLksXthdmy3OxYReZExxn9mML/88kvMnj0b\na9euxc6dOxEWFgYAOHDggFW3cEvt2rUxYsQIWzOQ/5PYGlDX3IKo8dNLmxzcPw7m23y7YxGRt/h4\nBtMfayTrI1WWtDoXMQ8tgCS3gp5xB/S/n4M5dNDuWETkDX+dgSkeHFNi2ePbBg0aID09HV988QXm\nzZuHP/74Az/++CP27duH0aNHW3UbzJs3Dxs3bkRcXBxmzZpV9vqmTZuwcOFCGGOQkZGBzMxMy+5J\nwSO0WTLUhBkwn7wP/cwsSJOzIZdeA0msY3c0IrJSSTEQ6rsZTF/USNZH8iYJDYPqMwimQ3eY5S9A\n3z0SctGVkPReEJuO/CEiL7CgPlpWXZs3b449e/Zg0KBBGDRoEA4fPoytW7di5cqVVt0CAJCRkYF+\n/frhiSeeKHtNa40FCxZg8uTJqFmzJiZOnIj27dujfv36lt6bgoMoBemYAZPaGWb1CujpY0sLaP8s\nSGSU3fGIyArFvl0i64sayfpIviBxNSHXZsP8WACd+yzMuyuhBl8PSUm1OxoRWaHYs+WxgIVLZBs2\nbIi2bduWfV2tWjWkpaXhpptusuoWAIDk5GRERR3/R35BQQHq1q2LxMREhIaGokuXLtiwYQMAYP/+\n/Xj66afxww8/YMWKFZZmocAmERFQF14ONeUxYH8h9N0jode+AVNSYnc0IvKUj5fI+qJGsj6SL8mZ\nSVDjpkNdfBX0P+ejZM5UmO0/2R2LiDzlYYMfwI0ZzJ07d+Lbb79Fly5dKvXz0dHRWL16NXr37l3l\ncJW1d+9eJCQklH0dHx+PgoKCsvsPHz78tO93uVxwuVxlX2dlZSEmJsY7YS0UHh7OnBY7adaYGOCW\nHBT/UIDDL86Ffv9NVBt6M8LanmdPSDjnd8qc1nNKVn/PeUAJwqKiAQC5ubllr6ekpCAlJcXt6/pb\njQzW+gj4/2fwGMfn7NobpnN3HH37FRx+eBJCz+uKaoOvhYqr6fuQcM7vE3BOVua0nj9n1cVHURga\nWpbPnRpZ5QFm7dq1AQCLFi1CrVq1kJKSggYNGkBEyn7m8OHDKCgowJYtWxATE4P+/ftX9TY+dbJf\nVmFhoU1pKi8mJoY5LXbarAlnwGRPAb74BAeenQMknlG6LKj+mT7NCDjnd8qc1nNKVn/PqQ8fRnFR\nEYDSQZNVAq1GOrU+Av7/GTwmYHKe3xfSrjOKVi7B0bHXQPoMgvS6CBIW7ruQcM7vE3BOVua0nj9n\nNX/+CaNCUFhYiJiYGLdqpFvzn7Vr18bQoUPx448/YsOGDXj55Zdx9OhRaK2hlEJcXBxatGiBgQMH\nIjo62p1bVEl8fDx2795d9vXevXsRHx/v9ftS8BERoG0HqJbtYNa+CT17EqRdp9JGB7E17I5HRJVk\nSoqhvLQH059qJOsj+ZJExUAu/wdM9/7QSxfC3H0zZNAwyHldj3vIQkR+zMMjSgAPm/xMmzYNAwYM\nQIcOHdC1a1ePglSFMea4M8SSkpKwY8cO7Nq1CzVr1kReXh6ys7OrfN1jS4GsfJpNgUlCwyC9LoLp\nlAGzcgn0PaMgvS4u/U9EhN3xiKgi5fZg5ubmerw09mTsqJGsj+QP5Ix6CLn5LphvtpY2AlrzGlTW\nDZCkc+yORkQVOWEPpjs10qMBZs+ePX3e7nzOnDnIz89HYWEhRo4ciaysLGRkZOD666/HtGnTYIxB\njx490KBBgypf2xt/YFBgK3tam9EfZtmL0JNGQC66AtK5J8SHDUSIqIpK/u4i661Bk69rJOsj+Rs5\nqyXUXbNg/vM+9NMPAWcmQV1yNaRO1T+DROQjJ8xg+myJ7DGJiYkASjvRffzxx4iOjkZKSgpiY2M9\nuexpnerJa2pqKlJT2SKb7CG160FGTID5/hvofy+EWf0K1KXXAK3bc1kQkT8q9v45mL6ukayP5I9E\nKUinDJhzO8OsWQk9405Iu86QgUMgNbhcm8jvWHAOpkfHlBz7wzk6Ohq9e/fGZ599hp07d3oUyE4u\nl+u4TklEVSVNzipt237ZtdDLXoB+aCLMtq/tjkVEJzphiWz5TqlWCaQayfpInpLwCKh+l0LdNxeo\nVg16yi3QryyGOXzQ7mhEVN4J50S7UyM9Gp6uWbMGxcXFOOuss9CoUSO0bNkSSUlJAIA//vgDNWo4\nq+kJlwCRFUQEaN2+tBHQuneh5z0INDsbatDVkDPq2R2PiACfLJENpBrJ+khWkehYyODrYXpcCLNi\nMXTOCMiALEjXvpDQMLvjEdEJ50T7fIlsdHQ0tmzZgiVLlqCoqAiJiYnYv38/2rVrB5fL5dWzL4n8\nnagQSHpvmPZdYda8Cv3geEj78yEXDmHHWSK7nVBAvYE1kujUJKE25IbbYX7+vrTj7DuvlnacTUvn\n1hIiO5UcP4PpDo/enZWVhWbNmgEAfvjhB+Tn52Pr1q1YunQpioqKHFc82SWPvEEiIiD9B8Oc3xfm\n9b86zvYcCOmdCYmoZnc8ouBUroB6q4tsINVI1kfyFmnYBCG3TYX58ovSHgZvLYe67FpIcmu7oxEF\npxMewPq8i+yxwgkAjRs3RuPGjdG/f38YY7B48WJPLm0LLgEib5KYWMiQ4TA9B8KsWFS6LGjgEEh6\nb3acJfK1cgXUW4OmQKqRrI/kbXJOG6ic2TAbPoR+/nGgTgOoS6+GNGhidzSi4HLCMSXu1MhKN/nJ\nz8/Hr7/+WqmfzcvLQ1paWpXDEAUDSawDNXwc1OgcmE8/gp4yGubTj2C0tjsaUfCwYAlQeayRRJ4T\npaA6dIO6dy4kJRX64cnQzz4Ks2eX3dGIgscJx5S4o8IB5po1a3DHHXdg3759qFevcg1K2rVrh6VL\nl+K+++7DBx984FFAokAljZtDjbkP6ooboVctg75/HIzr8+MOSSciL7FoDyZrJJH1JCwMqtdFUNPm\nAzVrQd93G/SSZ2AK/7Q7GlHAMxbUxwoHmE899RSSk5PRsWPHSl80MjISo0aNQkFBAebOnetRQF9i\nG3byNRGBtEiFypkN1e9S6Jeegp49iUebEHlbuXO+PDmmJFhqJOsj2UEio6AGDYWa+gRQUgx9983Q\nr/4T5hCPNiHympIiSMjfHZ29dkxJv379qhYMQI0aNdChQwe8//77VX6vXbjHhOwiIsC5XaDadoRZ\ntwZ6/gzgzCSozKGQ+o3sjkcUeCw8piQYaiTrI9lJ4mpCrhwB0zsT5tV/QufcBLngUkhGf0hYuN3x\niAKLBceUVDiDGR8fj7p161b5wgDQuvXJO4AdPXrUresRBToJCYE6vw/U9PmQ5i2gZ+eU7j/Z/bvd\n0YgCi0VLZFkjiXxHEutA3TAGasx9MN9shZ40AvrDt0uX9BGRNcqt8HFXhe9u0KCB2xc/cT/Kxx9/\njMLCQiQlJeHgwYNo2bKl29cmCmQSFg7pkwlzfh+Yt1dATxsD6dgd0n8wz9AksoJFTX5YI4l8Txo0\nRsjoSTDffQW97AWYt5dDXXwV0K6z3dGInM8XezBjYmLcvnh0dPRxX5eUlKBBgwZYv349tm7d6vZ1\nvYV7TMjfSPVIqIuvhLr3/wAAevIo6FcWwxzcb3MyIocrV0A92YMZLDWS9ZH8kTRLhho3Hery4dBv\nLoW+fxyKvtjAZnlEnjjhAaxX9mBGRERUPdhfwsOPXxcfHR2NRo0aoUWLFm5f05u4x4T8lcTWKD1D\ns9dFMK++hH23Xw30yYR07w8Jd//fKFEwMsYcd86XJ3swg6VGsj6SvxIRoGU7qBZtgY3rcGjh49Cx\nNaAGDYM0S7Y7HpHzWLAHs8IBZlhY2Cm/V1xcjI8++gjdunUr/Qd+gpATple/+uorbNmyBYcOHUJs\nbCyGDBlS5cBEwUxqnQG5/jZE/rEbhYufgnnnNcjAIZDOPSEW7CcjCgolJYBSJ61bVcUaSeQfRCkg\nLR3RXftg31sroJ+cCZzZ7K9meWfaHY/IOYo9PwfT7Q0ov/zyC+bMmYOffvoJH3/8MW6//XZUq1bt\ntO9p06YNzjnnHAClS4GIyD0hDZsgZNRdpftPlr943P4TURWufCcKbuVmL72FNZLIHsea5ZkO3WDW\nvgE9exIkpR3koisgiXXsjkfk/yyokW69e/Xq1XjhhRdw9OhRDBw4EOnp6ZgzZw5uuOEG1KpV65Tv\n27BhA95//33UqVMHLVu2RFJSktvBieiv/SdjpwFfboJe9iLw5lKoQcOAlFRLZmeIAlJJkccd8k6H\nNZLIfhIeAekzCOb8vqXN8qaPhZzXFXJhFiS2pt3xiPxXSTHg4fE/Vaqw+/fvx/z587FhwwbUqFED\n48ePL2uzPmLECDz55JMYNGgQmjdvftL3p6WloW7duiguLsZ3333H4klkAREBWqRCnVO6/0QveRqI\nrQl1ydXcf0J0MhYdUXIi1kgi/yPVIyEXXwmT0R/mjX9BTx4N6XYBpO8gSGR0xRcgCjYlxUC1SI8u\nUekBZn5+Ph5//HHs3bsXqampuPnmmxEbG1v2/bi4ONx+++14+umnsWvXLnTu/L+toss3LkhMTPQo\nuDe4XC64XC6PD90msoOIAOd2gWrbEebjd6Gfegho2KR0/0mDxnbHI/IfJ+mQ52kTm0CvkayP5HRl\nzfJ6Xwzz2kvQOSMgfQZBelwI8aBZF1HAOeEhrDs1ssIBptYaL7/8Ml555RUopXDNNdegf//+J/3Z\nsLAw3HzzzcjNzcXSpUvRp0+fSgfxB+ySR4FAQkIg6b1L95+8/yb0I5Mh57SBXHQlpLZ7B8ITBRQL\nOuQdEyw1kvWRAoUk1IZcmw3z28/QKxbD5NwEGZAFOb83JPTUTbuIgsYJD2G90kV2zZo10FqjXr16\nyM7ORuPGjSu8aFZWFj744AM8+eSTVQ5ERNaQsHBIr4th0nvDvPMq9APjIGnpkAGXQ2rE2x2PyD4n\nFE9PsEYSOZPUbYiQkXfC/PAt9PJFMKtXlDYCOq8rRLErOwUxC7aRVGoGs0ePHrjuuuv+58yu0+na\ntStq166Nr776yqOAROQZqRYJuXAITPf+MG8uhZ5yC+T8PpALLoFEuX9IPJFjWbgHkzWSyNmkcXOE\n3D4V5ust0MtegFm1DCpzKNDmPDbLo+BU4oNjSqZOnYrkZPcahSQnJ2PGjBluvZeIrCXRsZDB18H0\nHAjz+hLoSSMhvS6C9BwIqVbd7nhEvmPhDCZrJFFgkLNbQd05E9i8AXr5i8Cb/4YaNAyS3NruaES+\nZcExJRUemNegQQPk5+e7fYOEhISTvv7pp5+6fU0icp/E14IaNqq0kG7/EXrSCOg1K2GKiuyORuQb\nxdYNMFkjiQKHiEDanAc1+VFIxgDoF55AySOTYX741u5oRL5TXOzxUV4VDjCjo6OxY8cOrFy5Elpr\nj24GAEePHsXLL79c4YHTRORdckY9qBvHQ916D4xrI/TdI6HXrYHRPOCdAlxRERBmTTMP1kiiwCMq\nBKpjd6h7/w+S2gn6/6ajZN4DML/9bHc0Iq8zRUchHp6DWeEAEwB69OiBZs2aYebMmXj33Xdx+PDh\nKt9o//79eO211/DEE0+ge/fuaNmyZZWv4W0ulwu5ubl2xyDyKWnUFCG3Toa6YQzMh6uhp9wKs/Fj\nGGPsjkbkHcVFQLlukbm5uXC5XG5fLhhqJOsjBSMJDYPq3g9q2pOQJmdBz5wI/dwcmD077Y5G5D0W\n1EgxVfgrsqSkBGvXrsXatWsRERGBpk2bomnTpqhVqxYiIyMRGRkJYwwOHjyIgwcP4vfff8e2bdvw\n/fffIzw8HL1790ZaWlqVAtrl119/tTtChWJiYlBYWGh3jAo5JSfgnKzeymmMAbZ+Br3sRSAsrHT/\nyTlt3L5esP8+vcEpWf05p9m8Afq9NxCSfQ/q1atn2XWDpUY6oT4C/v0ZLI85reXNnObgfpi3VsC8\n/yakY3dI/8GQ2BpuX4+/U2s5JSfg31lLHroLauAQSHJrt2tklRbYhoSEoGfPnujZsyd27tyJLVu2\nYMuWLdizZw/279+PAwcOQEQQHR2N6Oho1K5dG0lJScjMzDzuwGki8k8iArRKg0ppB/NZHvSieUBC\nYulAs8lZdscjskaxdUtky2ONJApsEhkNGTQUpucAmDf+DT15FKRbP0jfQZDIKLvjEVmjuAjwcIms\n2zs4a9euXVZIiSiwiFKQ9ufDpHaCWbcGet6DwJlJUJlDIfUb2R2PyCOmqMjrB6qzRhIFLomtCRky\nHKb3xTCvvgQ9aQSkTyYk40JIRITd8Yg8U3T0uCWy7qjUHkwiCk4SGgrVtS/UtHmQ5udAz86BfvYR\nmN2/2x2NyH0WPJ0lIpKE2lDXZUONvx/m+2+hJ90EvfZNmOJiu6MRua+42ONVPhUOMA8cOIBNmzZ5\ndJOT+fjjjy2/JhF5h4RHQPUZBDX9SSDhDOhpY6D/+STMvv/aHY2o6ix4OnsMayQRSd2GCBl5J9So\nHJjPP4aefDP0+rUwFnSWJvI5X8xgRkVFYf/+/Vi6dClKSjw/vuDw4cNYtGgRatRwf1M0EdlDqkdC\nXXwl1H1zgZAQ6MmjoZe/CHNwv93RiCrPwmNKWCOJ6Bhp3Bwht98LdfVomPdeh743G+aLT9iVnZyl\nyEd7MNPT01FQUICZM2ciLS0NXbp0QWRkZJVutG/fPrz77rv47rvvMGzYMNSuXdutwERkP4mJg1z+\nD5heF8OsfBl60khI70xID+4/IQeweIksayQRlSfJraHunAl88Qn08heBN/9d2izv7FZ2RyOqmAWN\n8Crd5CcpKQkTJkzABx98gAceeAChoaFlLdgTExMRGRmJ6tWrl7VgP3ToUFkL9m3btqF69ero06cP\nMjMzPQpMRP5DEhIh19wCs+MXmBWLoSfdBOmfBTm/t9ebqBC5zcIlssewRhJReSICtO0A1ToN5pMP\noBc+BtSuB3XJMMiZSXbHIzo1C2pklbrIKqXQvXt3dO/eHbt27cKWLVvw5Zdf4sMPPzxtC/ZLLrnE\nES3YXS4XXC4XsrKy7I5C5ChSpwFkxASYHwtKl8yuXgG56AqYngPsjkb0v4qLgGp/zzDm5uYiJSUF\nKSkpHl02kGsk6yORe0SFQDpmwKSlw3y0GvqJaUCzZKiLhwIx59gdj+h/nbDKx50aKYYLw0/KCQdJ\n+/MhreU5JSfgnKz+ntN8vRV6+QtQR48AF10JtDmv9Gmun/L332d5Tsnqzzn1kmeAmrWg+mS6fYh0\nMHNCfQT8+zNYHnNay99zmiNHYN5dCfP2coSndUHxBZdBEhLtjnVa/v47PcYpOQH/zWpKSqBHXoqQ\np1YAgNs1kseUEJHl5OyWUBNmoPqQf0CvWAQ9YwLM11vsjkVUiseUEJFNJCICqt+lUNPnQ2omQN93\nG/TLT8Ps+8PuaESW7L8ELBpgfvvtt9i3b58VlyKiACEiCDu3M9TkOZDu/aGffxwlj9wD82OB3dEo\n2BUdBUKrtEPEI6yRRHQiiYxG9ctvgLr3CcAY6MmjoFcsgjl4wO5oFMws6lFgyQDzlVdewY4dO457\nbefOnVZcmogcTpSC6tgd6t65kNSO0E9MQ8n8B2F++8XuaBSsLGjBXhWskUR0KhJbE+qKG6EmPQz8\ndw/0pBHQby2DOXrE7mgUjCyqj5YMMNPT07F3715s374du3fvxu7du7F06VIrLk1EAUJCQ6G694Oa\n9iSkcXPohyZCL3wMZs8uu6NRkDHFRRAfDjBZI4moIlLrDKjrsqHGTYfZ9jV0zgjoD9+GseB8XaJK\ns2iJrCVrhJ555hnUr18fSv09Xt2+fbsVlyaiACMREZALLoXp2hfmrRXQ990G6ZRRerxJjH930qQA\nUVRk+TElp8MaSUSVJfUaIWTkRJjvv4X+93Mwq1+BuvQaoHV7v26WRwHCoiWyHg0wtdbYt28fhg0b\nhm7duh33vXXr1nkUjIgCm0RGQwYNhek5AGblEujJN0P6ZEJ6DoSER9gdjwKZRU9oK8IaSUTukibN\noZdYJc8AACAASURBVMZNB7Z8Cr30eeDt5VCXXgtperbd0SiQ2d3kZ+TIkVi1ahV+/vnn/ymcANC5\nc2ePghFRcJDYmlBXjoC6cybMDwXQd4+EXrcGRnNZEHmJRU9oT4c1kog8JSKQ1u2h7pkD6dQDet6D\n0PNnwOx0xlFB5EAWrfBxe4DZvn179O/fH61atTrp910ul9uhiCj4yBn1EDLyTqgb74D54C3o+8bA\nuD63OxYFIh80+WGNJCKriAqBSu8NNW0+0Kgp9APjof/5JI82IetZdIyX2wPMGjVqnPb7GzdudPfS\nRBTEpFky1IQZUAOHQL/0FEoemQzz0za7Y1EgKS4Cwrx7TAlrJBFZTSIioPoPhrp3LqAU9D2joFcu\ngTly2O5oFCjs3oP5+uuvY8OGDaf8/m+//YZhw4a5e3kiCmIiArTrBNW6PcxHb0M/NhXSoi3k4qGQ\nhES745HTFR0FQr07g8kaSUTeIjFxkCHDYXpcCLNiEfSkEZCBV0C69IKEhNgdj5ysyOYusq1atUKf\nPn1O+j1jDF599VW3Q1npyJEjeOaZZxAWFoYWLVogPT3d7khEVEkSGgrp3h+mY3eYt5aXdpxN7w3p\nfxkkMtrueORUFi0BOh0n1EjWRyJnk9p1ITeO/7vj7DuvQl1yNdDmPHacJbdYdYyX2wPMevXqoUWL\nFqf8/pdffunupS31n//8B506dUK7du3w6KOPsoASOZBUi4RcfBVMtwtgXn0JetLI0kFmt/4QH3QD\npQBT5P0lsk6okayPRIHhpB1nB18PaXKW3dHIaexeIvvFF19g0KBBCDvFH3eXXXaZ26FOZ968edi4\ncSPi4uIwa9asstc3bdqEhQsXwhiDjIwMZGZmAgD27t2LM888EwCOO4OMiJxHaiRArh4N0/Mi6GXP\nw6xZCRk0DNL+fD6tpcor9v45mHbUSNZHouAlIkDr9lAt28HkrYGeez/k7FaQQVdzawlVnt3HlPTo\n0QNr1qzBRx995HGIqsjIyEBOTs5xr2mtsWDBAuTk5GD27NnIy8srO8Q6ISEBe/bsAVC6LImInE/q\nN0LILXdDXXsrzNsroB8YD/PdV3bHIqc4ehTw8lmrdtRI1kciEhUCdX4fqPvmAYl1oO+7DXrF/7d3\n5+FRlWcbwO/nTUhCQgCHfVGRTSRWjWy2KBCkCBUr1TatFURRagMoIFsBEUFkR0BZVSq4G+uC1spX\na4uUWD9AjEpAMYhUNkkIS3YI7/v9MZqPSGAmyZl5z0nu33V5XZnJzDn3NYY8ec55l+dhigptRyMv\ncKg+VvoO5vXXX1/lk1dGhw4dkJWVVea5zMxMNGvWDI0a+a/QdO/eHVu2bEGLFi3QtWtXrF69Gtu2\nbUOnTp1sRCaiEJFLfwI1eQHMRxugV86FtE+A3DKEV2vpnExJCWA0EBHaIbI2aiTrIxH9QGJq+6eW\nXNcX5o3noKemQG6+HfKz3hDFhYDoHE4WA1EW52C6SU5ODho0aFD62OfzITMzEwAQHR2N4cOHn/f9\nGRkZZfYkS05ORnx8fGjCOigqKoo5HeaVrMz5IzfcDNOzL4refhknZ45BrZ/fjJibb4PE1A7q7V75\nPAHvZHVrTlOQj+PRMahbt27pc6mpqaVfJyQkICEhwUa0kKip9RFw78/gjzGns7ySEwhT1vh4YPQ0\nlGTuROFzy4EP3kXM4OGolZAY9CG88pl6JSfg3qyFAKROXcScka0yNbJaNJhVVd6HlZubaylN8OLj\n45nTYV7Jypzn0O/XkK49cfL1Z1E8ehBk4CDIT3tDAswv88rnCXgnq1tzmuNHgVpRpdni4+ORnJxs\nOZV7ebU+Au79Gfwx5nSWV3ICYc7apCXM2EeBbR8if8VcoMXFUL++C9K0RcC3euUz9UpOwL1Zdd4J\noL4Pp6pYI6vFrH6fz4fs7OzSxzk5OfD5fBYTEZEt4msEdc9YqJRJMP/+O/SjY2F2bbcdi9ziZHHI\n51+6CesjEf1ARCCdukPNWAZp1xF67gTol5+CyXdfo0OWOFQjPdlgGmPKLEjQtm1bHDp0CFlZWSgp\nKUFaWho6d+5c4eNmZGSUuQ1MRN4lrS+FmjgXcsOvoFcvwukVs2EOH7Qdi2wrp3impqaWGQbqZayP\nRBSI1IqCuuEWqOnLgNMl0FOHQ//jLf8cdarZHKqRnhsiu2TJEuzYsQO5ublISUlBcnIykpKSMHTo\nUMycORPGGPTu3RstW7as8LGr29wboppORCBde8Bc1Q3mvXXQs8ZBru0D+UUyJDbOdjyyoZziWV2G\nyLI+ElFFSN36kNtTYHrdCP3qapgP1kP9bhikAvMzqXoxJ4uhHKiRnmswR40aVe7ziYmJSEzkPwgi\nOptERUNuTIbp3gfmze9X0xs4CNK9T8D5mVTNVOMhsqyPRFQZ0uIiqFEPA59thX5hhX9+ZvLdkEZN\nbUejcHNomxL+ZXUGDgEiqt6kvg/qzlFQ9z8Ek/YP//6Ze3bZjkXhVM2HyIYK6yNR9SYikCu7QE1f\nBrmkPfSssdBvPA9TXGQ7GoVTTR0iG0ocAkRUM8jFbaEmzPHvn7lsFgoSu8HcdBukbn3b0SjUqvEQ\n2VBifSSqGaRWLcgvfgNzTRLMa2uhHxqOk4NSYC7vDBGxHY9CzaEayTuYRFQjiVJQP+vtX00vrg70\ntJFc5KAGMMXFEAc2kSYiqs7E1xBq2Fioe8ah+K2XoBdMhvl2j+1YFGoniwEHaiQbzDNwCBBRzSOx\ncag9eDjUhNkwn22BfmQ0zBef2Y5FocIhspXC+khUM0m7jqgzayWka0/oRQ9Bv7ASJu+E7VgUKhwi\n6zwOASKquaTZhVBjZgDb/gO95nHIJe0hv7kL4mtkOxo5iUNkK4X1kajmEhUB1bMfTOfuMOtehH5o\nBOSXt0F63ABREbbjkZM4RJaIyFn+Tah/5t8brGlL6EdGQ7+TCnPqpO1o5JRqvIosEVEoSVw81O/v\nhXrgEZjNG6FnjYf55ivbschJDtVINphERD8i0dFQN/8eavJCmG8yoR++D2bHJ7ZjkRPYYBIRVYm0\nbAU1fjak9wDopTP9w2YL8mzHIiewwXQe55gQ0ZmkUVNEjJgM9dt7oJ9dBv3kfJhjObZjUVVwDmal\nsD4S0ZlExL9Q3vSlgNHQD42E/uhfMMbYjkaVZEpKAGOAiLIzKDkHs4o4x4SIyiNXdIG69AqYd16B\nnn6/f+5Jz36ce+JFnINZKayPRFQeiYuHDBoO87ProV9YAbPpH1C3/xHS7ELb0aiivq+PP96OhnMw\niYhCRKKjoW65A2rcLJitm/xzT/Zm2o5FFcUhskREjpPWl0JNXghJ/Cn0vEnQr6+FKS62HYsqwsH6\nyAaTiKgCpMVFUONmQXrfCP34DOgXV8EU5NuORUEyJ4shbDCJiBwnERFQ1w+AmvY4cCQLetoImPT/\ntR2LgsUGMzQ4x4SIguGfe3K9f+5JySnoaSOgt/ybc0+8gHMwK4X1kYiCJfV9UMPGQQ25D/ova3B6\n2SyYo0dsx6JAztFgcg5mFXGOCRFVhNSpC7ljJEzmDujnV8B8+D7UoOGQBo1tR6NzKSoEYmqXeYpz\nMANjfSSiipLLroSa9jjM316FnjEKcvPvIT36QRTvb7lScdFZ9RHgHEwiIiukbUeoBxdB2naEnjkG\n+v23YfRp27GoPMVFQHSM7RRERDWC1Krl3/Zr3CyYjzZAz58Ec/Bb27GoPEWFjtVHNphERA6QyEio\nG5OhJs6F+TgNes5EmP3/tR2Lfqy4CIhhg0lEFE7S4iKoCXMgXXv4FwF66yWYU6dsx6IzFRcC0Wff\nwawMNphERA6Spi39iwB17wO9YDL0uhdZRN2knCGyREQUeqIUVNKNUFMXw/x3N/Qjo2Eyd9qORd8z\nRYUQh+ojG8wzcBEDInKCKAXVsx/UQ0tgvv2aRdRNyrlCy0V+AmN9JCKniK8h1IgpUL+8DXrlXOgX\nVsIUFtiORecY4cNFfqqIixgQkZPkggZQI6YAH6dBr5wLufqnkFsGQ2JibUerkczp08CpkrNWyeMi\nP4GxPhKRk0QE6Hwt1GVXwby2BvrhkVCDR0Iuv9p2tJqrqKjcIbJc5IeIyGVEBNL5WqjpTwAni6Af\nvh9m56e2Y9VMxUVAdLT/DxsiIrJO4upA3THSv6XJc8ugn13KvaVtKSp0bI0CNphERGEgcfFQd46C\nuj0F+pkl/iFBRYW2Y9UsnH9JRORK0jER6uEnABHo6ffBbN9mO1LNU+xcjWSDSUQURvKTTlAPPw6c\nLIaeMQrmy+22I9UcxeUP/yEiIvukdizU4BG8m2mLgzWSDSYRUZhJbB2ou0ZB/XYY9NMLoF9+Cqa4\n2Has6s/BPb6IiCg0Su9mArybGU7cB5OIyPvkyi7+IpqfCz3jfpjMHbYjVW8ODv8hIqLQkdqx/rmZ\nd3x/N3PtEzBFXGk2lExxEbcpCQUuw05E4SZx8VB3PwB1653QK+dBv/pnmFMnbceqnoqLyr06y21K\nAmN9JCIbJOH7u5nGQE8fxQuxoVRUUO4QWW5TUkVchp2IbJGrfwrVLgH6+eUwj46FumcspGUr27Gq\nlXNtIs1tSgJjfSQiW6R2LOTO+2E++Qh6xRzItX0hN/0OEsk2xlHFReWO8uE2JUREHibxdaH+OBHy\n84HQCx+Efm8djNa2Y1UfHCJLRORZkngN1ENLYL79GnrOBJiD39qOVL0UlT/KpzLYYBIRuYiIQHW/\nHmryApiP06AXT4M5dsR2rOqhqJCryBIReZjUuwDqvqmQa38OPe9P0P/8K4wxtmNVD9ymhIioepNG\nTaHGz4a0T4CeMRrm4w9tR/I+7oNJROR5IgLVqz/UxHkw//kX9JKHeSHWCQ7WSDaYREQuJRERUAN+\nBzXyQejX10KvWQJTyFX0Kq2wAKgdazsFERE5QJq2gJo4F3LJpdAzRuPklk22I3mWOX0aOHWSQ2SJ\niGoKaX0p1NTFgCjkThwG8/WXtiN5U2E+G0wiompEIiOhbv491IgpKHp2GfSLq7gSe2UUFQAxsRAR\nRw7HBpOIyAMkpjbUkPsQM+he6KUzof/+BhcAqiBTkA+JjbMdg4iIHCZtOqDOnKdgjh+Fnj0e5tA+\n25G8pcDZC7BsMImIPCSqaw//AkBb06CXzoTJPWE7kndwiCwRUbWl4ur4V2Lv2R967p+gP/yn7Uje\nUZgP1HbuAiwbTCIij5GGTaAmzIE0vxD6kdEwuyq2AXKNVVjgaAElIiJ3ERGonv2gxs6EWf8a9OpF\nMEVcuyCgwgIglncwQyIjIwOpqam2YxARBSSRkVC/vgtqUAr0qrnQ76RyyGwg55iDmZqaiowMNunn\nw/pIRF4iLVtBTVkIREZCP/IAzH93247kbue5g1mZGhnpRKbqIiEhAQkJCbZjEBEFTa7oAjXlMein\nF8Ds2g519xhI3Qtsx3Knc9zBTE5OthDGW1gfichrJDoGMuQ+6M0boRc/DBl4O+S6GxxbyKY6MQUF\nkHNMIalMjeQdTCIijxNfQ6ixj0IuaQ/9yBiYXdttR3InriJLRFTjqK49oCbMgfnnOzDPLIEpLrYd\nyX0cXqOADSYRUTUgERFQAwdBDbkPeuVc6PfWwRhjO5ZrmFMnAQOgVpTtKEREFGbStAXUpPmAPg09\nZzzM4QO2I7kLF/khIqJzkcs7QU2aD/PRv2CeWgBTVGg7kjt8f/eSQ6OIiGomiY6B3P0ApGc/6DkT\nYdI/sh3JPXgHk4iIzkcaNYWaOBeoFfX9fmD7bUeyr4BblBAR1XQiAtXrF1AjH4R+6Sno19fCnD5t\nO5Z9vINJRESBSFQ05M77Ib0HQM+dCPNJDb9Syy1KiIjoe9L6UqgHF8Hs3Q29eBr3lOYdTCIiCkbp\nfmD3TYV++Uno15+F0TX0Sm1hHhDLBpOIiPwkvi7UqGmQVu2gZ42F2feN7UjWmIJ8iIM1kg0mEVE1\nV3qldvdO6OWza+Sm0yY/D4irYzsGERG5iKgIqFuHQAYOgl74IMy2/9iOZEd+LhAX79jh2GASEdUA\nEl8PaswMSN36/sUNsg7ZjhRe+XkQB4snERFVH6pbT6hR06BfeQr6rZdgtLYdKbwKnL0IywaTiKiG\nkMhakMEjINfd4J+X+WUN2i8zP5d3MImI6JykVTuoyQthdnwCvWpuzVqFPT+PdzCJiKhyRATq+gFQ\nQ0dDr5oLvXG97UjhUeBs8SQioupH6l0ANfZRSO04/4XYGjDax+jTQBEX+Qna4cOHsXLlSjz22GO2\noxARuYp0TISaMAfm7+ugX3qy+i/Tnp8LxPIO5plYI4mIzia1akGG3Ae5tq+/ydz9he1IoVVYAMTE\nQlSEY4es1g1m48aN8cc//tF2DCIiV5KmLaAmz4f5bj/04zNgCqvv4j+GczDPwhpJRFS+0tE+d4yE\nXvYozMdptiOFTgimkEQ6erQQWbFiBbZt24Z69ephwYIFpc+np6djzZo1MMYgKSkJAwcOtJiSiMh7\nJLYO1H0Pwby4CnreJKj7H4Jc0MB2LOc5vEKem7BGEhGFhlzRBWr0dOilMyHZ30H6/goiYjuWs/Lz\nHB/h44k7mElJSZgyZUqZ57TWWL16NaZMmYKFCxciLS0N+/fvBwBs3LgRa9euxdGjR23EJSLyFImI\ngAxKgXTtAT1nQvXcC6wab1PCGklEFDpyUWuoP82D+WgDzPMrqt+UkhBcgPVEg9mhQwfExZXd/DMz\nMxPNmjVDo0aNEBkZie7du2PLli0AgB49emDIkCGoVasWnnrqKXzzzTd48803bUQnIvIEEYHqfyvk\nljugH5sKs/NT25Gc5fAKeW7CGklEFFriawg1cQ5MzmHopY9Uqykl/ikkNXCIbHlycnLQoMH/D+Py\n+XzIzMws85o6depg2LBhAY+VkZGBjIyM0sfJycmIj3f/HyJRUVHM6TCvZGVOZ3klJxCGrH0GoKR5\nS+QvmYGY2+9FVI8bKnUYt32mxwpyEd+0GSQq+qzvpaamln6dkJCAhISEcEYLCadqpFfrI+C+n8Fz\nYU5neSUn4J2szPm9+HiYSfNQ+MzjKFk4BXGT5kHV91XqUG76TItLTuH0BT7EniNPZWqkZxtMJ5X3\nYeXm5lpKE7z4+HjmdJhXsjKns7ySEwhT1gvbQMbORMGS6Sjc/y3kxuQKzzlx02dqiosAEeQVnwSK\nT5b5Xnx8PJKTky0lcz+v1kfAXT+D58OczvJKTsA7WZmzLJN8N8w7qTgxdQTU6OmQxs0qfAw3fab6\nSBYQFVNunsrWSE8MkS2Pz+dDdnZ26eOcnBz4fJW7ikBERGVJswv9c04+/hDmladhtLYdqfJOHAPq\n1LOdIqxYI4mIQkNEoAb8FnLDLdDzJ8H892vbkaom97jjNdIzDaYxBsaY0sdt27bFoUOHkJWVhZKS\nEqSlpaFz585VOkdGRkaZ28BERDWZ1PdBjX8UZm8mzDOLYUpKbEeqnLwTQPy5i2dqamqZYaBeFOoa\nyfpIRFSW6tkP6nfDoBdPg/lyu+04lZd7HIive85vV6ZGemKI7JIlS7Bjxw7k5uYiJSUFycnJSEpK\nwtChQzFz5kwYY9C7d2+0bNmySuepLnNviIicIrF1oEbPgF41F2bFbKh7J5Q7j9HVAhRPrw+RDUeN\nZH0kIjqbdOoOFVsHetVcqMEjIInX2I5UYSbvBNR5LsJWpkZ6osEcNWpUuc8nJiYiMTHRsfP8sJiB\n1//YICJykkRHQw2fDPPMEujF06BGToXExgV+o0uY3BOQAHcwvdxAhaNGsj4SEZVPLrsSatQ06Cdm\nQvJOQF3X13akisk9HnCUT0VrpCcazHDx8h8YREShJJGRwN1jgJefgl4wGWrMjPM2ba6Sd/7iyaYp\nMNZHIqJzk4vbQo2fBb3oIeiTJ6GuH2A7UvACNJg1apEfIiIKL1EKctsfID/pAr1gCsyJo7YjBeeE\n8wsYEBERnUmaNIcaPwvm/beg/+cN23GCYrQG8nMd3yeaDSYREQVNRKB+NQjS+Vro+VNgjuXYjhRY\ngDmYRERETpAGjaHGzYL599+h//qK7TiB5ecBMbH+UUoOYoN5Bq6SR0QUHHXT7yDX9IKePxkmJzvw\nGywyeYHnYHp9FdlQY30kIgqO+BpCjXsUZvNG6DefL7PCt+sEmEICVONVZMOFc0yIiIKnbkyGjoz0\nz8kcOxPSoLHtSOU7cYxzMKuI9ZGIKHj+bb5mQT82FSg5Bdx6J0TEdqyznQjcYHIOJhERhZW64RZI\n7xv9dzKzDtmOU77jR4F6PtspiIioBpH4ev47mTs/g3ltjSvvZJrjOZB6Fzh+XDaYZ+AQICKiilN9\nbob0HQj92FTXDZc1WvvnYNarf87XcIhsYKyPREQVJ3HxUGOmw2R8AvPWS7bjnO1YDlD//BdgOUS2\nijgEiIioclTvAdCnTkE/NhVqwiwg3tkV6Sot7wRQuzYkstY5X8IhsoGxPhIRVY7UqQs1Zgb0ginQ\nkZFQN7qo5hw/CgS4g8khskREZI264VeQbj2hH3sIOve47Th+HB5LRESWSd36UA88AvPhP6H//qbt\nOP/veE5IaiQbTCIicowM+C3kJ52RP2s8TEGe7TghK55EREQVIfV9UGNnwvzrHRS7ZJ9Mc/wopD7n\nYBIRkYuJCOSWOxBx6eXQj8+AKSq0msccPxqSBQyIiIgqSnwNocbORNHbr0Bves92nJCN8mGDeQYu\nYkBEVHUigtp3jIQ0vwh62aMwp07ZC3MsBwhwdZaL/ATG+khE5Axp2AR1Js+HefN5mPSP7IY5nhNw\nDiYX+akiLmJAROQMUQoyKAXmyfkwf14EDBsHURauaR47AjS78Lwv4SI/gbE+EhE5J6L5hVAjH4R+\nfAZUbDykffh/v5qiQuDUKSDu/IvycZEfIiJyDVERUHc/AHPiGMwrT1vZA8wcyYI0aBz28xIREZ2P\ntGoHdc9Y6JVzYPbtCX+Ao9mArxFExPFDs8EkIqKQkVpRUCOmwOzaDrP+tfAHyMkCfI3Cf14iIqIA\npONVkNv+AL1kBkzWofCe/EgW4GsYkkOzwSQiopCS2DioUdNgPlgPnfZ+eE+ek80Gk4iIXEt1uQ7S\n/1boxQ/DnDgWtvOanCxIiOojG8wzcBEDIqLQkPoNoEY/DPPGszCfbw3LOU1BPqA1EBt33tdxkZ/A\nWB+JiEJH9R4A6XKtf/X14qLwnDTIET5c5KeKuIgBEVHoSNOWUMMnQy+dCXX/Q5BW7UJ7whz/8J9A\n80u4yE9grI9ERKElN98O5GRDP/0YVMqfQr8wXk4WcOlPAr6Mi/wQEZGrSetLoQaPgF42CyYnO7Qn\nO5IFNODwWCIicj8RgdwxAijIg3ltbcjPZ45wiCwREVUTkngNpM9N0Esf8S+THiIm6wCkUbOQHZ+I\niMhJElkLavgkmE83Q29cH9qTHT4INA5NjWSDSUREYSd9fwW5uC300wth9OnQnOS7g0CT5qE5NhER\nUQhIXDzU/VNh1r0Ik/FJSM5hiouB/FzgAq4iS0RE1YSIQG7/I1BUCPOXNSE5hzl8EBKiq7NERESh\nIo2bQ907EXr1YzD7/+v8CbIOAA2bhGyeJxtMIiKyQiJrQaVMgvlsK/QHIRgKdPgA0Jh3MImIyHuk\nfQIkeSj0EzNgThx19uDfhW54LMAGswwuw05EFF4SV8c/FOitF2F2ODcUyJScAo7lAA0aB3wttykJ\njPWRiCj81DVJkJ/2hl4+G+bUSceOaw4fCHqED7cpqSIuw05EFH7SuDnUHyZAr5oLNWEOpGmLqh/0\n0D6gUVNIZOAyx21KAmN9JCKyQ276HXBoH8yzS4GhYwJuvRWUA/8NaosSgNuUEBGRR8mll0N+NRh6\n6UyY/LwqH8/s2wtpcbEDyYiIiOwRpSB3joI5uA/mb686ckx/jWzlyLHKwwaTiIhcQV3XF3L51dCr\n5sKcruLKsvu+AdhgEhFRNSDR0VAjp8BsXA/z8YdVOpYpKQEO7weaX+RQurOxwSQiIteQ3wwFlIJ5\n5ekqHcfs/wbSkg0mERFVD1K/AdTwKdDPL4fZu7vyB/ruAFC/ISQ62rlwP8IGk4iIXEMiIqD+MAFm\n56fQG96t1DGMMcDe3cCFbRxOR0REZI9c3AZq8HDoZY/CHDtSqWOYvZmQi1o7nKwsNphEROQqEhsH\nNfJB/8qyOz+t+AGyDgKRtSANGjkfjoiIyCK5+meQnv2glz4KU1xc8QPs/gJo08H5YGdgg0lERK4j\nTZpDDRsH/dQCmMMHKvRek/kFJMTFk4iIyBb5xW8gzVrCPLMYRusKvdfs3glpc1mIkvmxwSQiIleS\ny66E/PL30E/MhCmowMqyX2UAbUNbPImIiGwREcgdI2GOHYF5+6Wg32fyc4Hsw8CFrUIXDmwwiYjI\nxVSv/pDLroB+cn5QK8sarWE+/xhyeafQhyMiIrJEakVBDZ8M859/Qf/vB0G9x3z+MdDhJ5DIWiHN\nxgbzDBkZGUhNTbUdg4iIziC/HQZoDfPqnwO/+NuvgZjakCbNgz5+amoqMjIyqpCw+mN9JCJyH6lb\n379mwStPw+z+IvAbPtsCubJrhc5RmRrJBvMMCQkJSE5Oth2DiIjOIBERUPdOhPn8Y+gNfzvva81/\n/gXp9LMKHT85ORkJCQlViVjtsT4SEbmTtGwFNeR+6JVzYI5knfN1piAPJmNbhRvMytRINphEROR6\nElcHatRDMO+8Cp32frmvMUWFMB9tgFzXN8zpiIiI7JEru0B+fjP049NhThwt9zXmw39CEq6G1K0f\n8jxsMImIyBOkcXOosY/AvPk89P+87t/v8gzmb69CLr8a0rCJpYRERER2yM8HQjp1h547CebwwTLf\nM/l5MO/+BXLDLWHJwgaTiIg8Q5q2hJo0D+Z/P4BZNQ8m6xCMMdBp7/uvzv76TtsRiYiIwk5EoH55\nG+Tnv4SePR56w7swp07B5B6HXjEb0rUH5OI2YckSGZazEBEROUR8jaAmzoN591XoWWOBk8VAus0i\nfAAADG1JREFUkxZQo6dB6jewHY+IiMga1esXMO0SoF95GiZ1NSCA9PoF5NYhYcvABpOIiDxHoqMh\nAwfB3Hw7UFwIiYm1HYmIiMgVpMXFiHjgEZiTxYCKgESGt+Vjg0lERJ4lIgCbSyIiorNIVLSV83IO\nJhERERERETmCDSYRERERERE5gg0mEREREREROaJaz8HcsmULtm3bhsLCQvTu3RtXXHGF7UhERESu\nwBpJREShUK0bzC5duqBLly7Iz8/Hc889x+JJRET0PdZIIiIKBU80mCtWrMC2bdtQr149LFiwoPT5\n9PR0rFmzBsYYJCUlYeDAgeW+/7XXXkO/fv3CFZeIiChsWCOJiMhNPDEHMykpCVOmTCnznNYaq1ev\nxpQpU7Bw4UKkpaVh//79AICNGzdi7dq1yMnJwQsvvICrr74arVq1spCciIgotFgjiYjITTxxB7ND\nhw7Iysoq81xmZiaaNWuGRo0aAQC6d++OLVu2oEWLFujRowd69OiBd999F9u3b0dhYSEOHTqEPn36\n2IhPREQUMqyRRETkJp5oMMuTk5ODBg0alD72+XzIzMws85r+/fujf//+4Y5GRERkFWskERHZ4tkG\n00kZGRnIyMgofZycnIzmzZtbTBS8+Ph42xGC4pWcgHeyMqezvJIT8E5Wr+RMTU0t/TohIQEJCQkW\n07iLl+sj4J2fQeZ0lldyAt7JypzO80rWytRIT8zBLI/P50N2dnbp45ycHPh8vkodKyEhAcnJyaX/\nnflBuhlzOs8rWZnTWV7JCXgnq5dynvn7v7o0l07VSK/WR8BbP4NewJzO80pW5nSeV7JWtkZ6psE0\nxsAYU/q4bdu2OHToELKyslBSUoK0tDR07tzZYkIiIiI7WCOJiMgtPDFEdsmSJdixYwdyc3ORkpKC\n5ORkJCUlYejQoZg5cyaMMejduzdatmxpOyoREVFYsUYSEZGbeKLBHDVqVLnPJyYmIjEx0fHzeWWI\nFHM6zytZmdNZXskJeCcrc4ZPOGuklz4vr2RlTmd5JSfgnazM6TyvZK1sTjFnjqkhIiIiIiIiqiTP\nzMEkIiIiIiIid2ODSURERERERI7wxBzMUEhPT8eaNWtgjEFSUhIGDhxY5vsHDhzA8uXLsWfPHtx2\n220YMGCApaSBs27atAnr1q0DAMTExGDYsGG46KKLXJdz69ateOWVVyAiiIiIwJAhQ9ChQwfX5fxB\nZmYmpk6ditGjR6Nbt25hTukXKOuOHTswb948NGnSBADQtWtX3Hrrra7LCfj301u7di1Onz6NunXr\nYtq0aa7L+dZbb2HTpk0QEZSUlGD//v1YvXo14uLiXJWzoKAATzzxBLKzs6G1xk033YRevXqFNWOw\nWfPz87FixQp89913iIqKQkpKStgXm1mxYgW2bduGevXqYcGCBeW+5s9//jPS09MRHR2NESNGoFWr\nVmHN6DZeqZFeqY8Aa6TTWB+dxxoZ3pxuqI9AiGqkqYFOnz5tRo4caQ4fPmxOnTplxo0bZ/bt21fm\nNcePHze7d+82L730knn77bctJQ0u65dffmny8/ONMcZ88sknZvLkya7MWVRUVPr13r17zejRo8Md\nM6icP7xu+vTpZvbs2eajjz4Ke84fMgTKmpGRYebMmWMl3w+CyZmfn2/GjBljjhw5Yozx//tyY84z\nbd261cyYMSOMCf2Cyfn666+bF154wRjj/yzvuusuU1JS4sqszz33nHn11VeNMcbs37/fyme6c+dO\ns2fPHjN27Nhyv79t2zYza9YsY4wxu3btsvI71E28UiO9Uh+NYY20kZP1sWJYI8Of0w310ZjQ1Mga\nOUQ2MzMTzZo1Q6NGjRAZGYnu3btjy5YtZV5Tt25dtG7dGhEREZZS+gWTtX379oiNjQUAtGvXDjk5\nOa7MGR0dXfp1UVERRCTcMYPKCQDr16/HNddcg7p164Y94w+CzWosr9MVTM5NmzahW7dupRu92/hc\ng/08f5CWlobu3buHMaFfMDlFBIWFhQD8/5bi4+Ot/K4KJuu+fftw+eWXAwCaN2+Ow4cP48SJE2HN\n2aFDh/NeYd+yZQt69uwJwP87tKCgAMeOHQtXPNfxSo30Sn0EWCOdxvroPNbI8Od0Q30EQlMja2SD\nmZOTgwYNGpQ+9vl81opOIBXN+v777+Oqq64KR7Qygs25efNmjBkzBnPnzkVKSko4IwIILmdOTg62\nbNmCvn37hjveWTmC+Uy/+uorjB8/HrNnz8a+ffvCGRFAcDkPHDiAvLw8TJ8+HZMmTcLGjRvDHbNC\n/5ZOnjyJ9PR0K8O+gsnZr18/7Nu3D/feey/Gjx+PO++8M8wp/YLJevHFF2Pz5s0A/AU3OzsbR44c\nCWvOQLxUE8LBK5+HV+ojwBrpNNZH57FGOqu61EegcjWhRjaY1dX27duxYcMG3H777bajnFPXrl2x\naNEijB8/Hi+//LLtOOVas2ZNmc/Q9hXQ82ndujWWL1+O+fPno1+/fpg/f77tSOXSWmPPnj2YNGkS\nJk+ejNdeew2HDh2yHeuctm7dGvCKnk3p6em45JJLsGrVKsydOxerV69GUVGR7VjlGjhwIPLy8jBx\n4kSsX78el1xyCZRi6aHw8kJ9BFgjncT6GDqskc6ozvWxRi7y4/P5kJ2dXfo4JyendGiC2wSbde/e\nvXjyyScxefJk1KlTJ5wRAVT8M+3QoQMOHz6MvLy8sOYNJufXX3+NxYsXwxiD3NxcfPLJJ4iMjETn\nzp3DljPYrDExMaVfJyYm4umnn3blZ+rz+RAfH4+oqChERUXhsssuwzfffIOmTZu6KucPPvzwQytD\nf4Dgcm7YsKF0sYCmTZuicePG2L9/P9q0aeO6rLVr18bw4cNLH48YMaJ00Q238Pl8Za4aHzlyxLU1\nIRy8UiO9Uh8B1kgbOVkfnc/6A9ZIZ3J6oT4ClauR1aNNrqC2bdvi0KFDyMrKQklJCdLS0s77i9Hm\n1blgsmZnZ2PhwoUYOXJk2H8hVSTnmVfkvv76a5SUlIS92AeTc+nSpVi6dCmWLVuGa665Bvfcc0/Y\nm8tgs545Bj4zMxMAXPmZdunSBV988QW01iguLsZXX30V9pXSgv13X1BQgB07dqBLly5hzfeDYHI2\nbNgQn3/+OQD/z8DBgwetFKVgshYUFKCkpAQA8I9//AMdO3Ys84dfuBhjzvm7vHPnzvjggw8AALt2\n7UJcXBzq168fzniu4pUa6ZX6CLBG2sjJ+uh8VoA10smcbqmPgPM1UoxbxzaEWHp6Op555hkYY9C7\nd28MHDgQ7733HkQEffr0wbFjxzBp0iQUFhZCRBATE4NFixZZ+R8fKOvKlSuxefNmNGrUCMYYRERE\nYPbs2a7LuW7dOmzcuBGRkZGIiorC4MGD0b59e9flPNPy5cvRqVMnq9uUnC/r+vXr8d577yEiIgJR\nUVEYMmQI2rVr57qcgH958w0bNkApheuvvx79+/d3Zc4NGzbg008/xahRo8KeL9icR48exfLly3H0\n6FEA/mE21157rSuz7tq1C8uWLYNSCi1btkRKSkrpoivhsmTJEuzYsQO5ubmoV68ekpOTUVJSUub/\n++rVq5Geno6YmBikpKSgdevWYc3oNl6pkV6pj8FkZY10NifrY2iyskY6l9MN9REITY2ssQ0mERER\nEREROatGDpElIiIiIiIi57HBJCIiIiIiIkewwSQiIiIiIiJHsMEkIiIiIiIiR7DBJCIiIiIiIkew\nwSQiIiIiIiJHsMEkIiIiIiIiR7DBJCIiIiIiIkewwSQiIiIiIiJHsMEkIiIiIiIiR0TaDkBE9nz7\n7bf46quvsG/fPlx22WU4fvw4IiMj0atXL9vRiIiIrGF9JKo83sEkqsGOHDmCVq1aISsrC126dMF1\n112HN954w3YsIiIiq1gfiSqPDSZRDXbVVVfh008/RadOnQAAe/bsQXx8vOVUREREdrE+ElUeG0yi\nGu6zzz5Dx44dAQAbN27ETTfdZDkRERGRfayPRJXDOZhENVhRURGOHTuGnTt34rPPPkObNm3QrVs3\n27GIiIisYn0kqjw2mEQ12Pbt25GYmIiePXvajkJEROQarI9ElcchskQ11MGDB/HXv/4VJ06cQH5+\nvu04RERErsD6SFQ1YowxtkMQERERERGR9/EOJhERERERETmCDSYRERERERE5gg0mEREREREROYIN\nJhERERERETmCDSYRERERERE5gg0mEREREREROYINJhERERERETmCDSYRERERERE54v8A7NOjg8b2\nBjsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc9bf457668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15, 6))\n",
    "for idx, fi in enumerate([afi, nfi]):\n",
    "    plt.subplot(1,2, 1 + idx)\n",
    "    plt.semilogy(modelparams[:, 0], crb(fi))\n",
    "    plt.ylabel(r'$\\operatorname{Tr}\\left(\\left(\\sum_m F(p, m)\\right)^{-1}\\right)$')\n",
    "    plt.xlabel('$p$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we note that the numerical FI calculations are not much slower than the analytic calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 18.4 ms per loop\n",
      "10 loops, best of 3: 36.3 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit analytic_model.fisher_information(modelparams, expparams)\n",
    "%timeit numerical_model.fisher_information(modelparams, expparams)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  },
  "widgets": {
   "state": {},
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}