non_python_operations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial is generated from a [Jupyter](http://jupyter.org/) notebook that can be found [here](https://github.com/elfi-dev/notebooks). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using non-Python operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your simulator or other operations are implemented in a programming language other than Python, you can still use ELFI. This notebook briefly demonstrates how to do this in three common scenarios:\n",
    "\n",
    "* External executable (written e.g. in C++ or a shell script)\n",
    "* R function\n",
    "* MATLAB function\n",
    "\n",
    "Let's begin by importing some libraries that we will be using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.io as sio\n",
    "import scipy.stats as ss\n",
    "\n",
    "import elfi\n",
    "import elfi.examples.bdm\n",
    "import elfi.examples.ma2\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    ".. note:: To run some parts of this notebook you need to either compile the simulator, have R or MATLAB installed and install their respective wrapper libraries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## External executables\n",
    "\n",
    "ELFI supports using external simulators and other operations that can be called from the command-line. ELFI provides some tools to easily incorporate such operations to ELFI models. This functionality is introduced in this tutorial.\n",
    "\n",
    "We demonstrate here how to wrap executables as ELFI nodes. We will first use `elfi.tools.external_operation` tool to wrap executables as a Python callables (function). Let's first investigate how it works with a simple shell `echo` command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   3.,    1.,  123.])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Make an external command. {0} {1} are positional arguments and {seed} a keyword argument `seed`.\n",
    "command = 'echo {0} {1} {seed}'\n",
    "echo_sim = elfi.tools.external_operation(command)\n",
    "\n",
    "# Test that `echo_sim` can now be called as a regular python function\n",
    "echo_sim(3, 1, seed=123)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The placeholders for arguments in the command string are just Python's [`format strings`](https://docs.python.org/3/library/string.html#formatstrings).\n",
    "\n",
    "Currently `echo_sim` only accepts scalar arguments. In order to work in ELFI, `echo_sim` needs to be vectorized so that we can pass to it a vector of arguments. ELFI provides a handy tool for this as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.93678222e+00,   0.00000000e+00,   7.43529055e+08],\n",
       "       [  9.43846120e-01,   0.00000000e+00,   7.43529055e+08],\n",
       "       [  2.67626618e-01,   0.00000000e+00,   7.43529055e+08]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Vectorize it with elfi tools\n",
    "echo_sim_vec = elfi.tools.vectorize(echo_sim)\n",
    "\n",
    "# Make a simple model\n",
    "m = elfi.ElfiModel(name='echo')\n",
    "elfi.Prior('uniform', .005, 2, model=m, name='alpha')\n",
    "elfi.Simulator(echo_sim_vec, m['alpha'], 0, name='echo')\n",
    "\n",
    "# Test to generate 3 simulations from it\n",
    "m['echo'].generate(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So above, the first column draws from our uniform prior for $\\alpha$, the second column has constant zeros, and the last one lists the seeds provided to the command by ELFI.\n",
    "\n",
    "## Complex external operations $-$ case BDM\n",
    "\n",
    "To provide a more realistic example of external operations, we will consider the Birth-Death-Mutation (BDM) model used in [*Lintusaari at al 2016*](https://doi.org/10.1093/sysbio/syw077) *[1]*.\n",
    "\n",
    "### Birth-Death-Mutation process\n",
    "\n",
    "We will consider here the Birth-Death-Mutation process simulator introduced in *Tanaka et al 2006 [2]* for the spread of Tuberculosis. The simulator outputs a count vector where each of its elements represents a \"mutation\" of the disease and the count describes how many are currently infected by that mutation. There are three rates and the population size:\n",
    "\n",
    "- $\\alpha$ - (birth rate) the rate at which any infectious host transmits the disease.\n",
    "- $\\delta$ - (death rate) the rate at which any existing infectious hosts either recovers or dies.\n",
    "- $\\tau$ - (mutation rate) the rate at which any infectious host develops a new unseen mutation of the disease within themselves.\n",
    "- $N$ - (population size) the size of the simulated infectious population\n",
    "\n",
    "It is assumed that the susceptible population is infinite, the hosts carry only one mutation of the disease and transmit that mutation onward. A more accurate description of the model can be found from the original paper or e.g. [*Lintusaari at al 2016*](https://doi.org/10.1093/sysbio/syw077) *[1]*.\n",
    "\n",
    "<img src=\"resources/bdm.png\" alt=\"BDM model illustration from Lintusaari et al. 2016\" style=\"width: 400px;\"/>"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    ".. For documentation\n",
    "\n",
    ".. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.5/images/bdm.png\n",
    "   :width: 400 px\n",
    "   :alt: BDM model illustration from Lintusaari et al. 2016\n",
    "   :align: center"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This simulator cannot be implemented effectively with vectorized operations so we have implemented it with C++ that handles loops efficiently. We will now reproduce Figure 6(a) in [*Lintusaari at al 2016*](https://doi.org/10.1093/sysbio/syw077) *[2]* with ELFI. Let's start by defining some constants:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Fixed model parameters\n",
    "delta = 0\n",
    "tau = 0.198\n",
    "N = 20\n",
    "\n",
    "# The zeros are to make the observed population vector have length N\n",
    "y_obs = np.array([6, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype='int16')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's build the beginning of a new model for the birth rate $\\alpha$ as the only unknown"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Prior(name='alpha', 'uniform')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = elfi.ElfiModel(name='bdm')\n",
    "elfi.Prior('uniform', .005, 2, model=m, name='alpha')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g++ bdm.cpp --std=c++0x -O -Wall -o bdm\r\n"
     ]
    }
   ],
   "source": [
    "# Get the BDM source directory\n",
    "sources_path = elfi.examples.bdm.get_sources_path()\n",
    "\n",
    "# Compile (unix-like systems)\n",
    "!make -C $sources_path\n",
    "\n",
    "# Move the executable in to the working directory\n",
    "!mv $sources_path/bdm ."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    ".. note:: The source code for the BDM simulator comes with ELFI. You can get the directory with `elfi.examples.bdm.get_source_directory()`. Under unix-like systems it can be compiled with just typing `make` to console in the source directory. For windows systems, you need to have some C++ compiler available to compile it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 19.,   1.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,\n",
       "         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Test the executable (assuming we have the executable `bdm` in the working directory)\n",
    "sim = elfi.tools.external_operation('./bdm {0} {1} {2} {3} --seed {seed} --mode 1')\n",
    "sim(1, delta, tau, N, seed=123)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The BDM simulator is actually already internally vectorized if you provide it an input file with parameters on the rows. This is more efficient than looping in Python (`elfi.tools.vectorize`), because one simulation takes very little time and we wish to generate tens of thousands of simulations. We will also here redirect the output to a file and then read the file into a numpy array. \n",
    "\n",
    "This is just one possibility among the many to implement this. The most efficient would be to write a native Python module with C++ but it's beyond the scope of this article. So let's work through files which is a fairly common situation especially with existing software."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Assuming we have the executable `bdm` in the working directory\n",
    "command = './bdm {filename} --seed {seed} --mode 1 > {output_filename}'\n",
    "\n",
    "\n",
    "# Function to prepare the inputs for the simulator. We will create filenames and write an input file.\n",
    "def prepare_inputs(*inputs, **kwinputs):\n",
    "    alpha, delta, tau, N = inputs\n",
    "    meta = kwinputs['meta']\n",
    "\n",
    "    # Organize the parameters to an array. The broadcasting works nicely with constant arguments here.\n",
    "    param_array = np.row_stack(np.broadcast(alpha, delta, tau, N))\n",
    "    \n",
    "    # Prepare a unique filename for parallel settings\n",
    "    filename = '{model_name}_{batch_index}_{submission_index}.txt'.format(**meta)\n",
    "    np.savetxt(filename, param_array, fmt='%.4f %.4f %.4f %d')\n",
    "\n",
    "    # Add the filenames to kwinputs\n",
    "    kwinputs['filename'] = filename\n",
    "    kwinputs['output_filename'] = filename[:-4] + '_out.txt'\n",
    "    \n",
    "    # Return new inputs that the command will receive\n",
    "    return inputs, kwinputs\n",
    "\n",
    "\n",
    "# Function to process the result of the simulation\n",
    "def process_result(completed_process, *inputs, **kwinputs):\n",
    "    output_filename = kwinputs['output_filename']\n",
    "    \n",
    "    # Read the simulations from the file.\n",
    "    simulations = np.loadtxt(output_filename, dtype='int16')\n",
    "    \n",
    "    # Clean up the files after reading the data in\n",
    "    os.remove(kwinputs['filename'])\n",
    "    os.remove(output_filename)\n",
    "    \n",
    "    # This will be passed to ELFI as the result of the command\n",
    "    return simulations\n",
    "\n",
    "\n",
    "# Create the python function (do not read stdout since we will work through files)\n",
    "bdm = elfi.tools.external_operation(command, \n",
    "                                    prepare_inputs=prepare_inputs, \n",
    "                                    process_result=process_result, \n",
    "                                    stdout=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's replace the echo simulator with this. To create unique but informative filenames, we ask ELFI to provide the operation some meta information. That will be available under the `meta` keyword (see the `prepare_inputs` function above):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"67pt\" height=\"148pt\"\n",
       " viewBox=\"0.00 0.00 66.52 148.28\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 144.2833)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-144.2833 62.5152,-144.2833 62.5152,4 -4,4\"/>\n",
       "<!-- sim -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>sim</title>\n",
       "<ellipse fill=\"#cccccc\" stroke=\"#000000\" cx=\"29.2576\" cy=\"-22.8841\" rx=\"22.7687\" ry=\"22.7687\"/>\n",
       "<text text-anchor=\"middle\" x=\"29.2576\" y=\"-18.6841\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">sim</text>\n",
       "</g>\n",
       "<!-- alpha -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>alpha</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"29.2576\" cy=\"-111.0257\" rx=\"29.0172\" ry=\"29.0172\"/>\n",
       "<text text-anchor=\"middle\" x=\"29.2576\" y=\"-106.8257\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">alpha</text>\n",
       "</g>\n",
       "<!-- alpha&#45;&gt;sim -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>alpha&#45;&gt;sim</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M29.2576,-81.6758C29.2576,-73.5067 29.2576,-64.5751 29.2576,-56.1962\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"32.7577,-55.9774 29.2576,-45.9775 25.7577,-55.9775 32.7577,-55.9774\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x10baa8240>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create the simulator\n",
    "bdm_node = elfi.Simulator(bdm, m['alpha'], delta, tau, N, observed=y_obs, name='sim')\n",
    "\n",
    "# Ask ELFI to provide the meta dict\n",
    "bdm_node.uses_meta = True\n",
    "\n",
    "# Draw the model\n",
    "elfi.draw(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[13  1  4  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [19  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [14  3  2  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]]\n"
     ]
    }
   ],
   "source": [
    "# Test it\n",
    "data = bdm_node.generate(3)\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Completing the BDM model\n",
    "\n",
    "We are now ready to finish up the BDM model. To reproduce Figure 6(a) in [*Lintusaari at al 2016*](https://doi.org/10.1093/sysbio/syw077) *[2]*, let's add different summaries and discrepancies to the model and run the inference for each of them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Distance(name='d_sim')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def T1(clusters):\n",
    "    clusters = np.atleast_2d(clusters)\n",
    "    return np.sum(clusters > 0, 1)/np.sum(clusters, 1)\n",
    "\n",
    "def T2(clusters, n=20):\n",
    "    clusters = np.atleast_2d(clusters)\n",
    "    return 1 - np.sum((clusters/n)**2, axis=1)\n",
    "\n",
    "# Add the different distances to the model\n",
    "elfi.Summary(T1, bdm_node, name='T1')\n",
    "elfi.Distance('minkowski', m['T1'], p=1, name='d_T1')\n",
    "\n",
    "elfi.Summary(T2, bdm_node, name='T2')\n",
    "elfi.Distance('minkowski', m['T2'], p=1, name='d_T2')\n",
    "\n",
    "elfi.Distance('minkowski', m['sim'], p=1, name='d_sim')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"206pt\" height=\"341pt\"\n",
       " viewBox=\"0.00 0.00 205.55 341.29\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 337.2863)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-337.2863 201.5479,-337.2863 201.5479,4 -4,4\"/>\n",
       "<!-- d_T2 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>d_T2</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"168.774\" cy=\"-28.774\" rx=\"28.5497\" ry=\"28.5497\"/>\n",
       "<text text-anchor=\"middle\" x=\"168.774\" y=\"-24.574\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">d_T2</text>\n",
       "</g>\n",
       "<!-- alpha -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>alpha</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"98.774\" cy=\"-304.0287\" rx=\"29.0172\" ry=\"29.0172\"/>\n",
       "<text text-anchor=\"middle\" x=\"98.774\" y=\"-299.8287\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">alpha</text>\n",
       "</g>\n",
       "<!-- sim -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>sim</title>\n",
       "<ellipse fill=\"#cccccc\" stroke=\"#000000\" cx=\"98.774\" cy=\"-215.8871\" rx=\"22.7687\" ry=\"22.7687\"/>\n",
       "<text text-anchor=\"middle\" x=\"98.774\" y=\"-211.6871\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">sim</text>\n",
       "</g>\n",
       "<!-- alpha&#45;&gt;sim -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>alpha&#45;&gt;sim</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M98.774,-274.6788C98.774,-266.5097 98.774,-257.5781 98.774,-249.1992\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"102.2741,-248.9804 98.774,-238.9805 95.2741,-248.9805 102.2741,-248.9804\"/>\n",
       "</g>\n",
       "<!-- T1 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>T1</title>\n",
       "<ellipse fill=\"#cccccc\" stroke=\"#000000\" cx=\"28.774\" cy=\"-125.2755\" rx=\"19.8612\" ry=\"19.8612\"/>\n",
       "<text text-anchor=\"middle\" x=\"28.774\" y=\"-121.0755\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">T1</text>\n",
       "</g>\n",
       "<!-- sim&#45;&gt;T1 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>sim&#45;&gt;T1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M84.6084,-197.5505C73.822,-183.588 58.8848,-164.2525 47.1742,-149.0937\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"49.9304,-146.9364 41.0471,-141.1625 44.3909,-151.2159 49.9304,-146.9364\"/>\n",
       "</g>\n",
       "<!-- d_sim -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>d_sim</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"98.774\" cy=\"-125.2755\" rx=\"31.9567\" ry=\"31.9567\"/>\n",
       "<text text-anchor=\"middle\" x=\"98.774\" y=\"-121.0755\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">d_sim</text>\n",
       "</g>\n",
       "<!-- sim&#45;&gt;d_sim -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>sim&#45;&gt;d_sim</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M98.774,-192.5523C98.774,-184.7884 98.774,-175.8887 98.774,-167.1594\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"102.2741,-167.0684 98.774,-157.0684 95.2741,-167.0685 102.2741,-167.0684\"/>\n",
       "</g>\n",
       "<!-- T2 -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>T2</title>\n",
       "<ellipse fill=\"#cccccc\" stroke=\"#000000\" cx=\"168.774\" cy=\"-125.2755\" rx=\"19.8612\" ry=\"19.8612\"/>\n",
       "<text text-anchor=\"middle\" x=\"168.774\" y=\"-121.0755\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">T2</text>\n",
       "</g>\n",
       "<!-- sim&#45;&gt;T2 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>sim&#45;&gt;T2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M112.9395,-197.5505C123.7259,-183.588 138.6631,-164.2525 150.3737,-149.0937\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"153.157,-151.2159 156.5008,-141.1625 147.6175,-146.9364 153.157,-151.2159\"/>\n",
       "</g>\n",
       "<!-- d_T1 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>d_T1</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"28.774\" cy=\"-28.774\" rx=\"28.5497\" ry=\"28.5497\"/>\n",
       "<text text-anchor=\"middle\" x=\"28.774\" y=\"-24.574\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">d_T1</text>\n",
       "</g>\n",
       "<!-- T1&#45;&gt;d_T1 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>T1&#45;&gt;d_T1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M28.774,-105.2803C28.774,-94.477 28.774,-80.7559 28.774,-67.9254\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"32.2741,-67.8993 28.774,-57.8994 25.2741,-67.8994 32.2741,-67.8993\"/>\n",
       "</g>\n",
       "<!-- T2&#45;&gt;d_T2 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>T2&#45;&gt;d_T2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M168.774,-105.2803C168.774,-94.477 168.774,-80.7559 168.774,-67.9254\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"172.2741,-67.8993 168.774,-57.8994 165.2741,-67.8994 172.2741,-67.8993\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x118d33a20>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "elfi.draw(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.11 s, sys: 143 ms, total: 3.26 s\n",
      "Wall time: 5.56 s\n",
      "CPU times: user 29.9 ms, sys: 1.45 ms, total: 31.3 ms\n",
      "Wall time: 31.2 ms\n",
      "CPU times: user 33.8 ms, sys: 500 µs, total: 34.3 ms\n",
      "Wall time: 34 ms\n"
     ]
    }
   ],
   "source": [
    "# Save parameter and simulation results in memory to speed up the later inference\n",
    "pool = elfi.OutputPool(['alpha', 'sim'])\n",
    "# Fix a seed\n",
    "seed = 20170511\n",
    "\n",
    "rej = elfi.Rejection(m, 'd_T1', batch_size=10000, pool=pool, seed=seed)\n",
    "%time T1_res = rej.sample(5000, n_sim=int(1e5))\n",
    "\n",
    "rej = elfi.Rejection(m, 'd_T2', batch_size=10000, pool=pool, seed=seed)\n",
    "%time T2_res = rej.sample(5000, n_sim=int(1e5))\n",
    "\n",
    "rej = elfi.Rejection(m, 'd_sim', batch_size=10000, pool=pool, seed=seed)\n",
    "%time sim_res = rej.sample(5000, n_sim=int(1e5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Results after 100000 simulations. Compare to figure 6(a) in Lintusaari et al. 2016.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4Tdf6wPHvyiQxRYipSYh5TESEJtQshk60cltKUVqi\ntPQayu+iVHVAW6VVRV2K0quoOWKsuYomMSTmuYZMSEhkWr8/VnKaSCInycm8Ps9zHrL3Omu/MbxZ\nZ+213yWklGiapmnFn1lBB6BpmqblD53wNU3TSgid8DVN00oInfA1TdNKCJ3wNU3TSgid8DVN00oI\nnfA1TdNKCJ3wNU3TSgid8DVN00oIi4K6sL29vXR2di6oy2uaphVJx48fD5NSVs7Jewss4Ts7O3Ps\n2LGCurymaVqRJIS4mtP3ZjmlI4SwFkIcFUIECiFOCyGmZdCmgxDivhAiIPk1JacBaZqmaXnDmBH+\nY6CTlDJaCGEJHBBCbJNSHnmi3X4p5YumD1HTNE0zhSwTvlTlNKOTv7RMfukSm5qmaUWMUXP4Qghz\n4DhQF/hOSvlHBs1aCyGCgJvAWCnl6Qz6GQoMBahRo0aOg9Y0zTjx8fHcuHGD2NjYgg5FyyZra2sc\nHR2xtLQ0WZ8iO/XwhRAVgPXAe1LKU6mOlweSkqd9nge+kVLWe1pfHh4eUt+01bS8dfnyZcqVK0el\nSpUQQhR0OJqRpJSEh4cTFRVFrVq10pwTQhyXUnrkpN9srcOXUt4D9gDdnzj+QEoZnfz7rYClEMI+\nJwFpmmY6sbGxOtkXQUIIKlWqZPJPZsas0qmcPLJHCGEDeAMhT7SpJpL/RQkhWiX3G27SSDVNyxGd\n7IumvPh7M2YOvzqwLHke3wz4n5RysxDCF0BKuQDwAYYLIRKAGKCP1HsnFllSSn755RciIiJo0KAB\nDRo0wMHBQScOTSvijFmlEwQ0z+D4glS//xb41rShaQUhLi6O4cOHs2TJkjTH+/bty08//YSFRYE9\nq6cVY/v378fX1xdLS0sOHz6MjY1NQYdULOlaOppBZGQk3bt3Z8mSJUyZMoUbN26wa9cuxo8fz6pV\nq+jfvz8JCQkFHaZWREkpSUpKyvDcypUrmThxIgEBAUYl+6f1pT2FlLJAXi1atJBa4REfHy+bNWsm\nrays5PLly9OdnzlzpgRknz59ZHx8fAFEqOXEmTNnCvT6ly9flvXr15dvvvmmbNy4sVy6dKn09PSU\nzZs3lz4+PjIqKkouWrRI2tnZSWdnZ/nGG29IKdW/Nw8PD+ni4iKnTJmSYV9XrlyR27dvT9eflFLW\nrFlTTpkyRTZv3lw2bdpUBgcHSymljIqKkoMGDZJNmzaVLi4u8tdff5VSykz7KWgZ/f0Bx2QO867+\nfK4B8PPPPxMYGMiqVavo06dPuvPjxo1DSsmHH36Ik5MTM2fOLIAotdwYPXo0AQEBJu3Tzc2NOXPm\nPLXN+fPnWbZsGXXr1uXVV19l586dlClThi+++IKvvvqKKVOmcODAAV588UV8fHzw9/fn/PnzHD16\nFCklL7/8Mvv27aNGjRqGvjw9PQkLC+OTTz7JsD8Ae3t7Tpw4wfz585k9ezaLFy9m+vTp2NracvLk\nSUB9qs2qn+JEJ3yNhIQEpk+fTrNmzXjttdfSnb8VdQshBOPHj+fcuXN8/fXXvPXWWzRq1KgAotWK\nmpo1a+Lp6cnmzZs5c+YMbdq0AdT9Ii8vr3Tt/f398ff3p3lzdeswOjqa8+fPU6NGDUNfAEeOHHlq\nf6+++ioALVq0YN26dQDs3LmT1atXG9rY2dkZHVdxoBO+xsqVK7lw4QLr16/HzCztbZ2VQSsZvHEw\ncYlxNKnchDYvtMHG34bRo0fj5+enV+4UIVmNxPNKmTJlADV97O3tzapVq57aXkrJxIkTGTZsWJrj\nV65cMfRlTH+lSpUCwNzc/Kn3noyNqzjQN21LuJTRffPmzenZs6fhuJSST/Z9Qv/1/Wnt1JrPO3/O\nM+WeYdnpZZR+pzT+e/zZsGFDAUauFTWenp4cPHiQCxcuAPDw4UPOnTuXrl23bt1YsmQJ0dGqhNfN\nmze5e/dujvtLzdvbm++++87wdWRkZI76Kap0wi/hli9fzsWLF5k6dWqa0fqIrSOYvGcyb7q+yfb+\n2/nwuQ/xf9Ofrf22cjfpLhX6VuCDDz4gJiamAKPXipLKlSuzdOlS+vbti6urK15eXoSEhKRr17Vr\nV9544w28vLxwcXHBx8eHqKioHPeX2qRJk4iMjKRp06Y0a9aMPXv25KifIiund3tz+9KrdApeUlKS\nrFevnmzRooVMSkoyHN91aZdkKvIDvw/SHE/xfzv/TzIVSRPk559/np8ha9lU0Kt0tNwx9SodPcIv\nwYKCgjh//jy+vr6G0X2STGLcjnHUsK3Bp50/zXCOfmqHqXg5emHxigVfLvmSuLi4/A5d07Qc0Am/\nBPvtt98QQvDSSy8Zjq0+tZoTt04wo9MMrC2sM3yfpbklP/f+GatSVoQ2C+WXX37Jr5A1TcsFnfBL\nsA0bNtC6dWuqVq0KQGxCLP+36/9oXq05b7i88dT3OldwZsSzI6ARfLbgM6QunaRphZ5O+CXU1atX\n+euvv+jVq5fh2HdHv+Pq/avM8p6Fmcj6n8b7z76PuZk5weWC2bt3bx5Gq2maKeiEX0KlLKlMWYoZ\nmxDLpwc+pXvd7nSu3dmoPhzLO/Kvxv9CtBB88c0XeRarpmmmoRN+CfXbb7/RuHFj6tVTG5OtPbOW\niJgIxnqNzVY/454bh7SSbA/dXmzXLmtacaETfgkUERHBvn370kznLDqxiDp2dehYq2O2+nKv7k7r\nZ1qDJ3zz7TemDlUrxpydnQkLCzNJXwsWLOCnn34CYOnSpfz99995cp01a9bQpEkTzMzMKIpbtOqE\nXwJt2bKFxMREQ8I/G3aW36/+ztvubxs1d/+kie0nQnlY9ucyHj9+bOpwNe2pEhIS8PX1ZcCAAUD6\nhG9KTZs2Zd26dbRr1y5P+s9rupZOCfTbb7/h4OBAixYtAFh8YjEWZhYMchuUo/6er/c81UpV43b9\n22zZssVQtErTAHr16sX169eJjY1l1KhRDB06NF2b6dOns2LFCipXroyTkxMtWrRg7NixBAQE4Ovr\ny6NHj6hTpw5LlizBzs6ODh064ObmxoEDB+jbty9RUVGULVsWZ2dnjh07Rr9+/bCxseHw4cMAzJs3\nj02bNhEfH8+aNWto2LAhU6dO5fLly1y6dIlr167x9ddfc+TIEbZt24aDgwObNm3C0tIyTZxFvWCg\nTvglTGJiIjt27KBv376YmZnxOOExSwOX8nKDl6lWtlqO+jQTZrzT6h2mx07n+5Xf64RfSI0eDSau\njoybG2RVk23JkiVUrFiRmJgYWrZsSe/evalUqZLh/J9//snatWsJDAwkPj4ed3d3w2BkwIABzJs3\nj/bt2zNlyhSmTZtmKAIXFxdnmFaZOnUqAD4+Pnz77bfMnj0bDw8PwzUyKpUMcPHiRfbs2cOZM2fw\n8vJi7dq1zJw5k1deeYUtW7akmfYsDvSUTglz+vRpoqKiaNu2LQAbzm4g7FEY77i/k6t+B7gNAAG7\nQndx+/ZtU4SqFRNz586lWbNmeHp6cv36dc6fP5/m/MGDB+nZsyfW1taUK1fO8CDg/fv3uXfvHu3b\ntwdg4MCB7Nu3z/C+119/3egYUpdKvnLliuF4jx49sLS0xMXFhcTERLp37w6Ai4tLmnbFhR7hlzAp\nH3Fbt24NqJu1NW1r4l3bO1f91q1YF7dKbgS4BLBixQrGjs3eah8t7xVEdeS9e/eyc+dODh8+TOnS\npenQoQOxsbEm6Tt1qeSsZFYqOeW4mZkZlpaWhlIiZmZmxXI7zyxH+EIIayHEUSFEoBDitBBiWgZt\nhBBirhDighAiSAjhnjfharl16NAhqlSpQq1atbgdfZtdl3YxsNlAzM3Mc9330GeHQhVYsGGBfvJW\nA9Qo3c7OjtKlSxMSEsKRI0fStWnTpg2bNm0iNjaW6OhoNm/eDICtrS12dnbs378fUJVdU0b7T1Ou\nXLkMq2tqxk3pPAY6SSmbAW5AdyGE5xNtegD1kl9Dge9NGqVmMocPH6Z169YIIfgt5Dckkn81+ZdJ\n+n6tyWuYY87FMhc5ceKESfrUirbu3buTkJBAo0aNmDBhgmG3qtRatmzJyy+/jKurKz169MDFxQVb\nW1sAli1bxrhx43B1dSUgIMCobQcHDRqEr68vbm5uJi/fvX79ehwdHTl8+DAvvPAC3bp1M2n/eS47\npTWB0sAJ4Nknjv8A9E319Vmg+tP60uWR89/du3clIL/44gsppZSdl3WWDeY1yLAEck69sPwFyVjk\nyPdHmqxPLeeKSnnklE3DHz58KFu0aCGPHz9ewBEVDgVSHlkIYS6ECADuAjuklH880cQBuJ7q6xvJ\nx57sZ6gQ4pgQ4lhoaGh2fi5pJpDycdrLy4uwR2HsvbIXn8Y+Jt2mcHCLwVAWVh5ZSWJiosn61Yq3\noUOH4ubmhru7O71798bdXc8K5wWjbtpKKRMBNyFEBWC9EKKplPJUdi8mpVwILATw8PDQk7z57NCh\nQ1hYWODh4cHPwT+TKBPp3ai3Sa/xQr0XKGNWhkiHSA4cOGDUnKum/fzzzwUdQomQrWWZUsp7wB6g\n+xOnbgJOqb52TD6mFSKHDx+mefPm2NjY8Gvwr9S2q41bNTeTXqOURSleafQKNISff9H/iTWtMDFm\nlU7l5JE9QggbwBt4csPHjcCA5NU6nsB9KeUtk0er5Vh8fDxHjx6ldevWRMZEsvPSTnwamXY6J0Xf\nZn2hFKz+c3WxXNqmaUWVMSP86sAeIUQQ8CdqDn+zEMJXCOGb3GYrcAm4ACwC3s2TaLUcCwoKIiYm\nBi8vLzad20RCUgK9G5t2OidFl9pdKGNWhgeOD3SdfE0rRLKcw5dSBgHNMzi+INXvJTDCtKFpppT6\ngasRB0bgVN6Jls+0zJNrWZlb8UrjV1jxaAU///IzXbp0yZPraJqWPbq0Qglx6NAhHBwcsK9mz45L\nO3il4St5Mp2Toq9rX7CGNSfWEB8fn2fX0Yquolgeedy4cTRs2BBXV1deeeUV7t27l2nb8PBw3Nzc\ncHNzo1q1ajg4OBi+jouLY/DgwVSpUoWmTZuaJDZj6IRfQhw5cgQvLy9+v/o7sQmxPF/v+Ty9Xpfa\nXShjXoboGtHs3LkzT6+llWz5WR7Z29ubU6dOERQURP369fnss88ybVupUiUCAgIMFT8/+OADw9dW\nVlYMGjQIPz+/PIkzMzrhlwAPHjzg8uXLuLu7s/3CdqwtrGlXM+t63klJEBICy5fD/v3qa2NZmVvx\nauNXoSGsXrM6F9FrRV2vXr1o0aIFTZo0YeHChRm2mT59Og0aNOC5556jb9++zJ49G4CAgAA8PT0N\nI+rIyEgAOnTowOjRo/Hw8OCbb75h6tSpzJ49m19//dVQHjn1k7bz5s3D3d0dFxcXQkLUmpOpU6cy\ncOBA2rZtS82aNVm3bh3jx4/HxcWF7t27Z/jJtGvXrlhYqJlwT09Pbty4keM/l3bt2lGxYsUcvz8n\ndMIvAU6dUo9MuLi44HfRj/Y122NjaZNp+wcPYOBAsLeHRo1gwABo1w6cnFSJ3cuXjbtuH5c+YA3r\ng9br1Tol2JIlSzh+/DjHjh1j7ty5hIeHpzmfujzytm3b0uwkNWDAAL744guCgoJwcXFh2rR/Snml\nlEceM2aM4ZiPjw8eHh6sXLmSgIAAbGzUv/OU8sjDhw83/DABVR559+7dbNy4kf79+9OxY0dOnjyJ\njY0NW7ZsyfL76tGjR67+bPKbTvglQFBQEAB2tewICQuhW53M639cugReXvDzz9CrFyxeDIGB6utW\nrWDBAmjRAnbtyvq6Kat1opyiOHjwoKm+Ha2IKY7lkWfMmIGFhQX9+vUzOobCQCf8EuDkyZOUL1+e\nUzFqpN+97pPPzSn79qmkfusW+PvDkiUwZAi4ukLfvrB+PZw5Aw4O0K0bzJsHTyuKaWVuRc9GPaEh\nrFm/Ji++Na2QS10eOTAwkObNmxf58shLly5l8+bNrFy5Mk8XPuQFnfBLgJSPw9svbsepvBMN7Rum\na3PxIjz/vJrG+eMP6JjJXua1a8OhQ/DCC/D++zBp0tOv3a9ZP7CG//35P10yuQQqbuWR/fz8mDlz\nJhs3bqR06dKG4zdv3qRz5855ck1T0gm/mJNScvLkSZq4NmHX5V10r9s93agkMVHN2VtYwM6dUK/e\n0/ssV06N9t9+Gz79FFasyLxtl9pdKC1KE1ollL/++ssE35FWlBS38sgjR44kKioKb29v3Nzc8PVV\nz57eunXLcDPXWH379sXLy4uzZ8/i6OjIjz/+aNJYM5TTMpu5fenyyPnj2rVrEpD/nvNvyVTk2jNr\n07WZNUtKkPKnn7LXd1yclB06SFmqlJRHjmTe7vVVr0s+RE6cNDGb0Wu5pcsj54958+bJDRs2mLzf\nAimPrBVdKTdswyqEYS7M6Vwr7cfO06fhP/+BV16B/v2z17elJaxZA888o27wZrZC7U33N8EGfj6i\ni6lpGSvq5ZFHjhzJyy+/XNBhZEnvaVvMnTx5EoCgR0F4OXlha21rOCclvPUW2NrCDz9ATu4/2dvD\npk3g6QlvvqlW75g9MYzwruONNdZcLXuVCxcuULdu3dx8S1oOjfYbTcDtAJP26VbNjTndc79Zri6P\nnD/0CL+YCwoKwrGOI0GhQelG91u2wJ9/wsyZULlyzq/RpAl8/TXs3QvfZ7C5pZW5FS/UeQEawJp1\nerVOSRMTE0P79u3zdUOcZcuWUa9ePerVq8eyZcsMx/v06ZNuWWiJktO5oNy+9Bx+/mjatKls9UYr\nyVTkzos7DceTkqT09JTS2VnNxedWUpKU3bpJWbq0lBcvpj+/+exmyVRkw5ca5v5imtEKwxz+t99+\nK+fMmZNv1wsPD5e1atWS4eHhMiIiQtaqVUtGRERIKaXcu3evfPvtt/MtltzSc/ia0eLi4ggJCcGi\ntgXmwpxnHZ81nNuzB44cgQ8/VHPxuSWEekjLwkJNEz1ZhiFlWifELIRbt/RWCSXJypUr6dmzJwCz\nZs2iZcuWuLq68tFHHwFw5coVGjZsSL9+/WjUqBE+Pj48evQIgAkTJtC4cWNcXV0ZO3asUdfbvn07\n3t7eVKxYETs7O7y9vQ01a9q2bcvOnTtL7JPfOuEXYyEhISQkJBBRLgL36u6UtSprODdjBlSvDoMG\nme56jo5qamffvvRTO1bmVnjX8IaGsO63daa7qFaoxcXFcenSJZydnfH39+f8+fMcPXqUgIAAjh8/\nbnhy9uzZs7z77rsEBwdTvnx55s+fT3h4OOvXr+f06dMEBQUxKfmhj5UrVxqqTqZ++fj4AGpNvJPT\nPxvwOTo6cvOm2oDPzMyMunXrEhgYmM9/EoWDTvjF2MmTJ8EcLj2+RNsabQ3HDx+G3bth7Fiwtjbt\nNd96C7p0USt/7txJe26I1xCwgSV7lpj2olqhFRYWRoUKFQDw9/fH39+f5s2b4+7uTkhIiGE+3cnJ\niTZt2gDQv39/Dhw4gK2tLdbW1gwZMoR169YZHnTq16+foepk6tevv/5qVExVqlTJs2qahZ1O+MVY\nUFAQ5jXMiUuKo23NfxL+p59CpUowbJjprymEKrnw6BFMmJD2XLe63bCUlvz1+K+n1hHXig8bGxtD\nKQUpJRMnTjQk6AsXLjBkyBCAdA8DCiGwsLDg6NGj+Pj4sHnzZkOdm6xG+A4ODly/ft3Q140bN3Bw\ncDB8HRsbayiqVuLkdPI/ty990zbv9ejRQ1b7VzXJVOTd6LtSSimvXpVSCCmnTMnba48frx7mOnQo\n7fEuC7pIxiGXr1ietwFoUsrCcdPW0dFRxsTEyO3bt8tWrVoZHrK6ceOGvHPnjrx8+bIE5KHkfyxD\nhgyRs2fPllFRUfLOnTtSSinv3bsnK1asaNT1wsPDpbOzs4yIiJARERHS2dlZhoeHG843bdpU3rp1\ny8TfZd7QN201o506dQozZzMa2jekchm17nLZsn/W3+elyZNVkbURI1TphhRDnxsKZWDRjkV5G4BW\naHTt2pUDBw7QtWtX3njjDby8vHBxccHHx8dQ86ZBgwZ89913NGrUiMjISIYPH05UVBQvvvgirq6u\nPPfcc3z11VdGXa9ixYpMnjyZli1b0rJlS6ZMmWKoO3/nzh1sbGyoVq1ann2/hVpOf1Lk9qVH+Hkr\nOjpaIpDWU63lOxvfkVKqpZO1a0vZsWP+xLB6tRrlz5+fKq7H0dLiIwtp8bKFfPToUf4EUoIVhhH+\n8ePHZf/+/TM9f/nyZdmkSZN8ieWrr76SixcvzpdrmUK+j/CFEE5CiD1CiDNCiNNCiFEZtOkghLgv\nhAhIfmVd4UjLU+fPn4cqEEus4Ybt/v2q3n1ej+5TvPYadOigRvspU/ZlrMrgWcmThHoJ+O/wz59A\ntALl7u5Ox44d8/XBq8xUqFCBgQMHFnQYBcaYKZ0EYIyUsjHgCYwQQjTOoN1+KaVb8utjk0apZdu5\nc+eghvp9yg3b//5XVbpM3gsi56RUj+j+5z/w8ceqdOaFC+mK4wuhlmlGRMD06f8cH9Z2GJSDRdv0\ntE5JMXjwYMzNzTM85+zsbNiVLa+99dZb2a5qWZxk+Z1LKW8Bt5J/HyWECAYcgDN5HJuWC2fPnoWa\n4FDOgZq2NYmOVoXO+vSBbOwbkVZiInz+Ofz4o9rn0NxcPWGVkujbtoWFC6HhP/X23dzUJipz56pV\nQfXrQ89GPTFbb8bu27tJTEzMNBFommZa2bppK4RwBpoDf2RwurUQIkgIsU0I0SST9w8VQhwTQhwL\nDQ3NdrCa8c6dO4e5sznP1XgOIQS//goPH+ZiOicmBnx81I4ndeqo7bBCQyEqCo4eVUP5U6egWTOY\nNg3i4gxv/eQTsLFR6/4BypUqh3t5d2JqxXDo0KHcf7OaphnF6IQvhCgLrAVGSykfPHH6BFBDSukK\nzAN+y6gPKeVCKaWHlNKjcm6qdWlZOnX5FIllE3nWQZVTWLpUja5bt85BZ2Fh0LkzbNgAc+bAjh3q\nJ4ednfq40LKl2t08OBh694apU+GllyB5/XXVqmr2Z9Mm9VaAd9q8A7awYNMCk3y/mqZlzaiEL4Sw\nRCX7lVLKdM/FSykfSCmjk3+/FbAUQtibNFLNaFJKzkWfA6DFMy24dUuVO+jXLwclkB8+hPbt4cQJ\nNSc0Kt09+39Urap2O1+8WG2K6+NjGOmPGgW1asG//w0JCfCa22uIJMHWK1uRUm99qGn5wZhVOgL4\nEQiWUma4EFYIUS25HUKIVsn9hpsyUM14d+/e5VGFRwgEzas1Z8MGNc2eo5u148apkfvGjWr0bowh\nQ1QxnS1b1E2D+HisrWHWLDXrs3gxVLCuQBObJtx75l6+3bDTChdnZ2fCwsJM0teCBQv46aefALXJ\neOrSCaa8zrhx42jYsCGurq688sorT31iPDw83PAUcLVq1XBwcDB8ffXqVTp27Ejjxo1p0qQJ33zz\njUniy1JW6zaB5wAJBAEBya/nAV/AN7nNSOA0EAgcAVpn1a9eh5939u3bJ+mDdPrcSUoppbe3lPXq\nqXX42bJ1q1pIP2ZMzgKZM0e9/913pZTq+u3aSWlvL2VkpJSzd8+WTEUO/3h4zvrXslQY1uFnpmbN\nmjI0NDTX/cTHx6f5un379vLPP/80+XWklHL79u2G640fP16OHz/eqPd99NFHctasWYav//77b8M2\njg8ePJD16tWTp0+fTvc+U6/DN2aVzgHgqRMBUspvgW9z+kNHM62zZ8/CM2o6JzJSlUIeMyab0znh\n4TB4MDRtqu665sSoUXD9Onz5JXTqhOjdmzlzoEUL1eWEjwcydu9Y1p9fz3zm5+waWqHXq1cvrl+/\nTmxsLKNGjWLo0KHp2kyfPp0VK1ZQuXJlnJycaNGiBWPHjiUgIABfX18ePXpEnTp1WLJkCXZ2dnTo\n0AE3NzcOHDhA3759iYqKomzZsjg7O3Ps2DH69euHjY0Nhw8fBmDevHls2rSJ+Ph41qxZQ8OGDZk6\ndSqXL1/m0qVLXLt2ja+//pojR46wbds2HBwc2LRpE5ZP1A7v2rWr4feenp5GF2x7UvXq1alevToA\n5cqVo1GjRty8eZPGjTNa8W46urRCMXTi/AkoD23rtGXzZjVn/sor2ezk3XdV0l+xInclNT/9FFq1\nUtM8V67QvLn6OTJ3LkTetKeuRV1u293m2rVrOb+GVqgtWbKE48ePc+zYMebOnUt4eNrZ3j///JO1\na9cSGBjItm3bOHbsmOHcgAED+OKLLwgKCsLFxYVp06YZzsXFxXHs2DHGjBljOObj44OHhwcrV64k\nICDAUCTN3t6eEydOMHz4cGbPnm1of/HiRXbv3s3GjRvp378/HTt25OTJk9jY2LBly5Ysv68ePXrk\n6s8G1H4Af/31F88++2zWjXNJJ/xi6PjfxwFo5diKdetUTZuWLbPRweHD8L//qUdkmzXLXTBWVrBq\nlbqJ0LcvxMfzySdQqpS6PdDfvT/Yww/rf8jddbRCa+7cuTRr1gxPT0+uX7+ebovBgwcP0rNnT6yt\nrSlXrhwvvfQSAPfv3+fevXu0b98egIEDBxrq5wO8/vrrRsfwavINrBYtWnDlyhXD8R49emBpaYmL\niwuJiYmGipwuLi5p2j1pxowZWFhY0K9fP6NjyEh0dDS9e/dmzpw5lC9fPld9GUMn/GLoUswlkFCv\nnBvbt6vR/ZMbiz/V1Klqd/IPPjBNQLVrqzu1R47A9OlUq6aWaW7YAE3MfEHCqqBVprmWVqjs3buX\nnTt3cvjwYQIDA2nevLmhXHJulcnGE4SlSpUCwNzcPM1uVynHzczMsLS0NJRpNjMzy3RXrKVLl7J5\n82ZWrlyZrqxzdsTHx9O7d2/69etn+IGU13TCL2YSEhIItw7HXtpzcE9ZYmKyOZ1z8KBaUvnhh1C2\nbNbtjfWvf0H//upJ3eBgRo8GZ2f4+MOqOCXV4LL1ZSIjI013Pa1QuH//PnZ2dpQuXZqQkBCOHDmS\nrk2bNm1BFIv/AAAgAElEQVTYtGkTsbGxREdHs3nzZgBsbW2xs7Nj//79ACxfvtww2n+acuXKGapw\nmpqfnx8zZ85k48aNhg1ZQO2y1blzZ6P7kVIyZMgQGjVqxL///e+8CDVDOuEXM1euXEFWk9QvV591\n66BiRWjXLhsdfPQRVKkCw4ebPrgvv1Q/RHx9sS4lmTULTp6E+omDoAos/m2x6a+pFaju3buTkJBA\no0aNmDBhAp6enunatGzZkpdffhlXV1d69OiBi4sLtra2ACxbtoxx48bh6upKQEAAU6ZkXZdx0KBB\n+Pr64ubmRkxMjEm/n5EjRxIVFYW3tzdubm74+voCcOvWrWzV6Dl48CDLly9n9+7dhqWaW7duNWms\nGcrp8p7cvvSyzLyxbN0yyVTkeytHywoVpBw4MBtv/v13tYzyq6/yKjwpFy1S11iyRCYlSdm2rZR2\nNa9LpiKbDm+ad9ctoQrzsszUUjZFefjwoWzRooVhyWJRMW/ePLlhwwaT95vvyzK1omXfeXVTq4Z5\nd+7dgxdeyMabp06FatUgedSSJwYPVruwjBuHeOkl5syxx8PDkdLhjTnDGWJjY7E29Ua7WqE3dOhQ\nzpxRf/8DBw7E3d29oEPKlpEjRxZ0CEbRCb+YCbgbAGUg9NRzmJmpDcWNcuqUWrA/c6aqdJZXzMzg\nhx/U6p+JE3FftIihQ2Hh8cHIrmNZuXUlQ14dknfXL8lGj4aAANP26eam6ivl0s8//2yCYLSs6Dn8\nYuZK3BVsHtqwa1sZvLxUfTOj/PCDWkKZH7ujNG4MI0eqiptBQXz2GdjdVmUbfti/MO+vr2klVU7n\ngnL70nP4ecNsnJmsObqBBCk/+cTIN0VHS1m+vJT9+uVpbGmEh0tpZydlly5SJiXJ5cul5B0PaTas\nmkxISMi/OIq5wjCH/+jRI9muXbt8/XtdunSprFu3rqxbt65cunSp4fjrr78uz507l29x5JbexFzL\n1M3ImySVSaLcI/XEXvIzJFn75Rd48EDtUJJfKlZUK4J27oRt2+jXD6qEdyOp+m2W/fb0Jxy1omXJ\nkiW8+uqr+bbRTUREBNOmTeOPP/7g6NGjTJs2zbDkd/jw4cycOTNf4iiMdMIvRnad3AVA4t+dqFIF\nmjc38o0//ACNGsFzz+VdcBkZPhzq1YMxYxAJ8Xz/ng8AE5fvzd84tDy1cuVKevbsCcCsWbNo2bIl\nrq6ufPTRR4BaStywYUP69etHo0aN8PHx4dGjRwBMmDCBxo0b4+rqytiUHXSysH37dry9valYsSJ2\ndnZ4e3vj5+cHQNu2bdm5c2emD1UVdzrhFyOHL6pCUdePd6B7dyOfrv3rL7Vjla9vDorl55KVlaqZ\nHBICixfzagc3rCJqcrfSH+zerWvkFwdxcXFcunQJZ2dn/P39OX/+PEePHiUgIIDjx48bSiWcPXuW\nd999l+DgYMqXL8/8+fMJDw9n/fr1nD59mqCgICZNmgSoHyApa9dTv3x81IDh5s2bODk5GWJwdHTk\n5s2bgHqCtm7dugQGBubzn0ThoBN+MRJ0JwgemRN9q4bx0zk//KCKo735Zp7GlqmXX1Z74X78MTx8\nSDenFlDjEANGXuLx44IJSTOdsLAwKlSoAIC/vz/+/v40b94cd3d3QkJCDHV1nJycaNOmDQD9+/fn\nwIED2NraYm1tzZAhQ1i3bp3hydZ+/foREBCQ7mVs5coqVaqkqZdfkuiEX4xcfngZEeqEmRmkquKa\nuZgYtUPVa69lYzmPiQkBn30Gt2/DvHlM7Kk+tt8st41URQ21IsrGxsZQO0dKycSJEw0J+sKFCwwZ\nopbgPlmTRgiBhYUFR48excfHh82bNxsKm2U1wndwcOD69euGvm7cuIGDg4Ph69jYWEMVzRInp3d7\nc/vSq3RMKykpSVpMtpDmL70uPT2NfNPateqp1x078jQ2o7zwgpQVKkgZESHLjC8jzYa4SGvrJHnh\nQkEHVrQVhlU6jo6OMiYmRm7fvl22atXK8FTtjRs35J07d+Tly5clIA8dOiSllHLIkCFy9uzZMioq\nSt65c0dKKeW9e/dkxYoVjbpeeHi4dHZ2lhERETIiIkI6OzvL8PBww/mmTZvKW7dumfi7zBv6SVst\nQ9fuXyPBPAH+bo/380a+afVqVTenQ4e8DM04M2aoh3hmzaJDlQ5ssd6CWblbDBv2DDt25P/tBc10\nunbtyoEDB+jatSvBwcF4eXkBULZsWVasWIG5uTkNGjTgu+++Y/DgwTRu3Jjhw4dz//59evbsSWxs\nLFJKvvoqwx1W06lYsSKTJ0+mZXJN8ClTplCxYkUA7ty5g42NDdWqVcubb7awy+lPity+9AjftDaG\nbJRMReJ0QO7ZY8QbHjyQ0tpayhEj8jo04/XtK2Xp0vLA779KpiKbDBwuQcoffyzowIquwjDCP378\nuOzfv3+m5y9fviybNGmSL7F89dVXcvHixflyLVPQ6/C1DB2+pFbomEc0JIOChOlt3AixsWqT8cLi\n448hLo7Wv+zG5pENl62X0a6d5N//hhJ6j61YcHd3p2PHjiQmJhZ0KFSoUIGBAwcWdBgFRif8YuLY\ntWNwrxpN60rjdiRcvRocHaF16zyPzWh168KQIYhFi/ARLXhU9RHDx/zJ48cwYoTaNEsrmgYPHpzp\ng1fOzs6cOnUqX+J46623slXGuLjJMuELIZyEEHuEEGeEEKeFEKMyaCOEEHOFEBeEEEFCiKJV6q4Y\nOB16Fu62oEMHI7JiRARs365G99naCisfTJ4M5uZ8droUmMGqv2bw8cfw22+Qw/2iNU1LZsz/9gRg\njJSyMeAJjBBCPLm1eg+gXvJrKPC9SaPUnio+MZ7bCTfhblN69jRiX8z16yE+vnBN56RwcID33+eZ\nLbtxv2bJ7lu7+eADaNFC1Vt7Yv9rzQhSfzQqkvLi7y3LhC+lvCWlPJH8+yggGHB4ollP4KfkewpH\ngApCiOomj1bL0NnwsySJRLhbn+eeK5X1G1avVtMnhbXm+IcfIsqXZ+5hW6KrRPNX8J/8+KP6YJKP\nu8EVC9bW1oSHh+ukX8RIKQkPDzf53hDZmswSQjgDzYE/njjlAFxP9fWN5GO3nnj/UNQnAGrUqJG9\nSLVMnbqr5j/LxVhhaZlF44gIVfd+/PjCu9axYkUYN442kybh+TfM3DCTNZPWMGECfPIJ9O2bjcJw\nJZyjoyM3btwgNDS0oEPRssna2hpHR0eT9ml0whdClAXWAqOllA9ycjEp5UJgIYCHh4cecpjI4Usn\nIcmcenZGJPAtWyAxEXr1yvvAcmPUKOTcuXy2I5Re3f0BmDRJzeMPG6b2aylXroBjLAIsLS2pVatW\nQYehFRJG3bETQliikv1KKeW6DJrcBJxSfe2YfEzLB/tCAiGsAa3cY7NuvGEDVK8OHh55H1hulC2L\nmDSJDtckraIeEBgSSKlSas+U69fh//6voAPUtKLHmFU6AvgRCJZSZvao20ZgQPJqHU/gvpTyViZt\nNRM7f/8U3G1Amzaln94wNhb8/FTBssK2OicjQ4fysHplPt0LX6z9HAAvL3j/ffj2WzhwoECj07Qi\nx5j/9W2AN4FOQoiA5NfzQghfIUTKbtdbgUvABWAR8G7ehKs9KSY+hoeW1yC0FPXq1X5649274eFD\nSK5NXuiVKoXNjC/wuAVWv28yHP7kE3B2hiFD1M8wTdOMY8wqnQNSSiGldJVSuiW/tkopF0gpFyS3\nkVLKEVLKOlJKFynlsbwPXQPYd/o8CAlh96hdO4uEv3EjlC0LnTrlT3AmYDZgAJcq2zDh+EOCQ04C\n6ltYuBDOnVMP52qaZpwi8Llee5oNh4IBKBV1E3t7+8wbJiWphN+9O5QyYulmYWFuzvXRQ2kYAQcn\n/bMFo7e32m995ky1h4umaVnTCb+IO3Q2BKSgrl36muJpHDsGt24VnemcVNp8OIuj1aHb9iNp5nC+\n/BLs7dVmXYWgTIumFXo64RdhUsL5yBDE/arUr1336Y03bABzc3je2NrJhYeFuSWre7jiFC25MmWC\n4bidnUr6R4/CokUFGKCmFRE64RdhFy/CozLBEFrKuPn7du3UQ01FUI/Rs9lZC+y+mw9RUYbjb7wB\nHTvCxIlw924BBqhpRYBO+EXYrt1JUOksMvTB0xP+9evqSaUXXsi/4EzM28Wb6a1tsH0Uj/zyS8Nx\nIWD+fLX4aNy4AgxQ04oAnfCLsM0HroFlLIRFPj3h+/mpX3v0yJ/A8kg1j5dY2wgSZs2EsDDD8YYN\nVbL/6SfYt68AA9S0Qk4n/CJKSjgQolboEAZ16tTJvLGfHzg5QaNG+RNcHvn4Xx8zqQOYxcTA9Olp\nzv3nP1CjBowapW/galpmdMIvooKD4Z55iPoiDGrWrJlxw/h42LlTLccsrMXSjNTAoQGR1lVZ7Apy\n/nwICTGcK11aLdEMCID//rcAg9S0Qkwn/CJqzx7APgSrBGtq2NfAysoq44ZHjsCDB8WmvGTfxn2Z\n7A1xVuYwdmyac6+9Bm3aqNH+gxyV99O04k0n/CJq924o5RhMqSjLrKdzLCygc+f8Cy4PTX5tMqHm\n8FXL0qry544dhnNCwJw5arXOjBkFGKSmFVI64RdBSUmwdy+IyiEk3E54+g3bbdvUvrW2tvkWX16q\nWL4iTeKbMLV1JPE1a6gdURISDOc9PGDgQJX4L14swEA1rRDSCb8ICgqCiJhwYs1DibkWk3nCv31b\n1R0oJtM5KcZ3GU9cKZjXxUktN/3hhzTnP/1UfaiZNKmAAtS0Qkon/CIoZf4egDAyT/j+auOQ4pbw\n+3frT6k7pZhc6QSyc2c1aX/njuH8M8+o1TqrV6ubuJqmKTrhF0G7d0OVJv8k/Ezn8P38oGpVaNYs\n/4LLB2ZmZvSo3INHpWPwG/YiPHqU7qmrceOgQgX1s0DTNEUn/CImIUE9XFStaQgWWMC9TEb4iYlq\nhN+tW9HY7CSbPnnjE4iBiSFL1f68y5fD778bztvZwYQJsHWr3ihF01IUv0xQzJ04oZYcmlcNxjbB\nlvLlylMxo/o4x45BeHixm85J0aRBE6qHVicoPojQUe+oHVGGD4e4OEOb995TuzlOmKAeVNO0kk4n\n/CJmzx71a6R5CFb3rahdu3bGZZH9/NQ6RW/v/A0wHw3zGIY0l3y8azbMm6eeRps503C+dGmYPBkO\nHlSLlTStpNMJv4jZvRsaucRyLeoycbfinj5/36qVKhhfTI1+YzRm18xYdnoZic/3gNdfV1tgnTxp\naDNkCNSsqQ7rUb5W0umEX4TExan5aLdOF0iSSdy/eD/j+fvwcFUkvphO56SwtbWlbam2RFlEsSZo\njRrlV6gAgwapkhKAlZUqnfzHH6rChKaVZFkmfCHEEiHEXSHEqUzOdxBC3E+1wfkU04epgcrhjx6B\nQzNVNC3hViYPXe3cqZ7OKuYJH2Bqn6lwH6b5TYPKleH779WNjlRTO4MGgaNjunprmlbiGDPCXwpk\nlTn2p9rgXG8rnUf27FHT8pbVk5dkhmeyQsfPTy1TadkyfwMsAO3btqfylcqEPA7h9N3T0Lu3Kqoz\nbZp6Qg21he/48bB/f5qFPJpW4mSZ8KWU+4CIfIhFy8Lu3eDmBlcfhVDJvBLEZ7AGPylJJfyuXdWW\nhsWcEIKRrUdCQvIoH+Dbb9XOXn37qo9EwNtvq0cS9ChfK8lMNYffWggRJITYJoRoYqI+tVRiYuDw\nYbWdX3BoMBUSKmBmZkaNGjXSNgwKUiUVSsB0Top3B72L2Skz1l9aT2RMpJraWb5crdoZPRoAGxv1\nMNauXerPUdNKIlMk/BNADSmlKzAP+C2zhkKIoUKIY0KIY6GhoSa4dMlx+DA8fgwdOiZxNvwsFvcs\nqFGjBpaWlmkbpuxu1a1b/gdZQOzt7elSvgsJIoH5R+erg97eagH+okXwyy8A+Pqqgf8XXxRgsJpW\ngHKd8KWUD6SU0cm/3wpYCiEyXAsopVwopfSQUnpUrlw5t5cuUfbsUTM0td1u8Cj+EXF/x2U+f+/m\npp44KkHGDxgPl+DL/V8Sl5j88NW0aeDlBe+8AxcvUqYMjBih9nM/e7Zg49W0gpDrhC+EqCaSn/wR\nQrRK7jM8t/1qae3erUr/3ohVK3Qiz0Wmn79/8EA9ZVSCRvcpOnbsyDNXnyEyMZJVJ1epg5aWsGqV\nKp3ZqxdERzNypFqqmWofdE0rMYxZlrkKOAw0EELcEEIMEUL4CiF8k5v4AKeEEIHAXKCPlPoRF1OK\njlZLMjt2hJAwtULn3oV76Uf4e/aoYjslaP4+hZmZGeN6jYM7MH33dAz/BGvWVFM6Z87AgAFUsU9i\n4EC14XmqApuaViIYs0qnr5SyupTSUkrpKKX8UUq5QEq5IPn8t1LKJlLKZlJKTynlobwPu2Q5cEDl\n8U6dVMIvb1keHmawJNPPD8qWVRuelECDBg3C6rgVF6Musv3i9n9OeHurIf369fDJJ4wZox5i+/bb\ngotV0wqCftK2CNizR81OtGkDwWHBPGP5DPBEwpdSJfxOndScRQlUoUIFBjQfAFHw6e+fpj05apTa\nCuujj6gfuIaePWH+fHj4sGBi1bSCoBN+EbB7N3h6qmJgIWEhlI8vDzyxBv/8ebhypUTO36f2/oj3\n4Qjsv7GfE7dO/HNCCFiwQP3U7N+f6V1+JyICli4tsFA1Ld/phF/I3bunKgV07AiRMZHceXgH80hz\nKlSogJ2d3T8NtydPYZTA+fvUXFxcaFOqDSJO8On+J0b51tZqiU6dOjT5T09eb3ySb7/VRdW0kkMn\n/EJu3z718GzK/D3A4xuPM56/r1sXnraheQkxevho5BHJuuB1BIcGpz1ZsSL4+SHKlOHHWz14FHKV\nXbsKJk5Ny2864Rdye/aogamn5z8JP+JcRNrpnMePYe/eEj+6T9GzZ0+euf4MIlHw+cHP0zeoUQP8\n/CidFM3vZh1ZNfN6/gepaQVAJ/xCbvduNe1cqpRK+FbmVtw8fTPtCP/AAVUzpoTP36ewtLRkjO8Y\nko4msTJoJZciL6Vv5OKC2LGDapbhTNzRkWuHbuR/oJqWz3TCL8Tu3FGlcbp0UV8HhwVTq1wt4h/H\np034fn5qGU+HDgUSZ2H09ttvU/ZkWWSiZObBmRk3atmSB2v8qcJdbF7sBDdv5m+QmpbPdMIvxFI2\n7EhJ+CFhIVSzqAY8sSRz+3Zo21atwdcAKF++PMP7DSfpRBL//eu/3HyQcTKv8tKzfNFxO9aRt0jq\n2Alu3crnSDUt/+iEX4jt3KnuMTZvDrEJsVyMvEj5uCeWZN68qbb009M56bz//vuYHzEnITGBzw9k\nMJefrOtHXnTHj8RrN9VyqNu38zFKTcs/OuEXUlLCjh3QubMqmnY+/DxJMgnzCHPMzc1xcnJSDf39\n1a/6hm06jo6OvNHjDcyCzFh4fCHX7l/LsF27dhDl0oahjtuQN26oJVG67oJWDOmEX0iFhKjBu7e3\n+vpM6BkAYq/HUrNmTSwsLNQJPz9VGdPFpYAiLdzGjBlDwu4EEpMSmbFvRoZthID33oOlF9sS9NlW\nuHpV3Q/R0ztaMaMTfiG1Y4f6NXXCNxNmhJ0N+2f+PjFRNezWTWUtLZ1mzZrxwnMvYBFkwZKAJRmv\n2AHeeEPtfz5jfzvYtg2uX1dJX9/I1YoRnfALqR07oE4dcHZWXweHBVPbrjZXLlz5Z/7+zz8hMlLP\n32dh8uTJPN7xGJJg+r6M9zgsUwaGDIF16+BG7XbqRvjff6ukf0Mv2dSKB53wC6H4ePUcVcroHtQI\nv36F+oSFpRrhb9+uRvapG2rpPPvss3h7emMZYMlPgT9xLvxchu3efVc91fzDD6iHH/z91Vx++/Zw\nLeP5f00rSnTCL4T++EPVwE/J4wlJCZwLP0cVsyoAaRN+y5ZQqVIBRVp0TJ48mZgdMVhIC6bunZph\nm9q14YUXYOFC9fAyXl7qo1Z4uBrpX72anyFrmsnphF8I7dgBZmZqhSDAxYiLxCfFUy62HJCc8CMj\n1U8GPZ1jlLZt29Leoz1WAVasPrWaU3dPZdjuvffg7l343/+SDzz7rFofGxmpRvpXruRbzJpmajrh\nF0I7dqjtDFOKYaas0BFh6sZsnTp1VBJKStLLMbNh8uTJRG+PppQolekov0sXqFcPvv8+1UEPD9i1\nC+7fV+tk//47X+LVNFPTCb+QCQ9XA/fUeTwl4T+89pCKFStia2urlmPa2kKrVgUUadHTqVMn2rdq\nj/mf5qwNXkvA7YB0bczMYPhwOHwYAlKfdndXf+Z376qfCqGh+Re4ppmITviFjL+/Grg///w/x4LD\ngqlhW4PrF66r6ZyU3a26dFEbdGtGEULw2Wef8XDnQ6yxZsqeKRm2GzQIbGyeGOWDmt7ZskVN63Tt\nqjYr0LQiRCf8QmbrVrC3V7MIKc6EnqGRfSMuXryoEn5AgJpWePHFggu0iPLy8uIl75eQByWbzm3i\njxt/pGtjZwd9+sDKlWoWJ4127dTeuKdPQ8+eEBubP4FrmglkmfCFEEuEEHeFEBne5RLKXCHEBSFE\nkBDC3fRhlgyJiWrg3r27KqcAkCSTCAkLoUHFBly+fJkGDRrApk1qOWbqjwGa0WbMmMHjfY8pLUvz\n4c4PkRlsefXuu2q/2+XLM+igWzf46Se1O03//uovTtOKAGNG+EuBp90Z7AHUS34NBZ78IKwZ6dgx\nCAuDHj3+OXb13lViEmKolFSJpKQk6tevD5s3q+mFKlUKLtgizMXFhX4+/YjfEc/vV39n87nN6dp4\neKjX/PmZbIHYpw98/TWsXas2SNf7JGpFQJYJX0q5D4h4SpOewE9SOQJUEEJUN1WAJcm2bWrgnnql\nZcoNW6v7VgA0qVRJPWGrp3Ny5eOPP0Yel5SPL8+HOz8kISkhXZt334XgYDWQz9Do0TBuHHz3Hcye\nnbcBa5oJmGIO3wFIvUfcjeRjWjZt3aq2Mkz9HFVKwn988zEADS5eVCd0ws+V2rVrM/q90TxY94Dg\nsGCW/LUkXZvXX1fz+fPnP6Wjzz9XDT/8EH77Le8C1jQTyNebtkKIoUKIY0KIY6F6WVsad+6ogfuT\n0/LBYcFULVOVG+dvUKVKFUrv3g2OjuDqWjCBFiOTJk3CPsyecvfK8dHej4iOi05zvnRptWJn3bqn\nlMg3M4P//lctj+3XD06cyPO4NS2nTJHwbwJOqb52TD6WjpRyoZTSQ0rpUblyZRNcuvjYvl39+mTC\nPxN6hkaVG3H27Fma1q2r1m2++KKujmkCtra2fDrjU6LWRnE7+naGWyH6+kJCAixe/JSObGxgwwa1\nvOqll3SFTa3QMkXC3wgMSF6t4wncl1LqQuLZtHUrVK0Kbm7/HEuSSZy6ewqXKi6cO3eOF8uVU0tH\nXnqp4AItZgYPHkyzis0ofak0sw7N4sq9K2nO16+vahr98INK/JmqWlXdTL9/H3x8kovxaFrhYsyy\nzFXAYaCBEOKGEGKIEMJXCOGb3GQrcAm4ACwC3s2zaIupuDi1HPP559UMQYqLERd5GP+QeuXrcefO\nHdpFRanRZEqRHS3XzM3N+eabb3j02yMS4hMY6z82XZvhw1WF5C1bsujMxQWWLYMjR9TKHU0rZIxZ\npdNXSlldSmkppXSUUv4opVwgpVyQfF5KKUdIKetIKV2klMfyPuziZc8eNTB85ZW0x1Me/beNsQWg\n0fnz6ulaG5v8DrFYa9++Pf1e6ofcJ1kbvJY9l/ekOf/SS+DgkMXN2xS9e8OECeojwVPngTQt/+kn\nbQuB9evVBhxPlrUPvBOIuTAn8VYi7kDp0ND0PxU0k/jqq68od7Ic1jHWjPIblWaZpoUFDBumbp+c\nP29EZ598okovjBihCiNpWiGhE34BS0xUq/l69ABr67TnAu8E0sC+AVcuXKE3IM3N9fx9HqlSpQoz\nP51J7MZYTt49yfw/0w7n335bJX6jRvnm5rBqlfpY0Lu33hBdKzR0wi9gR46ofPDqq+nPBd4OpFnV\nZpw7d47XLC0R7durlSBanhgyZAitK7bG8qol/9n1H24++Ge1TfXq8K9/wZIlanOaLFWsqD66RUSo\nN8bH513gmmYknfAL2Pr1YGmZfjlmREwE1x9cp1nVZsQFBlI3Pj7jnwqayZiZmbHwh4XIzZKYxzF8\nsP2DNOffew8ePMikvk5GmjWDH3+E/fthzBjTB6xp2aQTfgGSUiX8zp1VafvUgu4EAeBa1RXXCxfU\nwV698jnCkqdJkyZ8/MHHJO5JZM2ZNWw7v81wztMTWrSAefOyUTqnb1/44AP1JqN/Umha3tAJvwCd\nPAmXLmV8HzbwdiAA1UQ1XoyP53atWmpOWMtz48aNo2V8S8wizBi2aRiP4h8B6lm3995T9XV27cpG\nhzNnqj1xhw2DwMA8iVnTjKETfgFav14lkZ49058LvBNI5dKViTkaTAvgwZNLeLQ8Y2FhwfKly7HY\nZsH1qOtM2j3JcO7119VtlHnzstUhrF6tCvO8+qraH1fTCoBO+AVo7Vpo3Vo9pPmkwDuBNKvWTP1U\nAMr275/P0ZVsDRo0YOaImfAnzDkyh4PXDgJqJdXQoWpLgsuXs9Fh1aqwZg1cuwYDBqhtzTQtn+mE\nX0BOnlSvvn3Tn0tISuD03dM0q9qM6ocOESQE1dq0yf8gS7j33nuPTomd4B688b83DFM7w4erlZdz\n52azw9atVQ39zZvh009NH7CmZUEn/AKycqVKGq+9lv7c2bCzPE58TOtEB2r9/Td7q1bFzEz/VeU3\nMzMzVi1bhd1+O649vMZ4v/GAKlb6+uvqQdpsb2s7YoSqqjllinqSS9Pykc4iBSApCX7+WW10klHR\n0MA76sae14ErAJx117tGFpQqVaqwdvZaOArfHf+OfVfUbihjxqj1+AsXZrNDIVTZhaZN1ce7K1dM\nHrOmZUYn/AJw4ABcv64GehkJuB2ApbCg8jp/DgDVvbzyNT4trQ4dOjCx1US4Bz4rfHgY95DmzaFT\nJ6d+YaUAABw5SURBVDWtExeXzQ7LlFFF9hMTVWVNvRG6lk90wi8AK1ao//MZrc4BNcJ/Jb4OFiEh\nrARc9WYnBW76pOm0+rsVoYmhDFw+EFCj/Js34X//y0GHdeuqjdCPH1d3gfWeuFo+0Ak/nz1+rBZr\n9Oqlkv6TpJQc+/sYg05ZkWhmxhqgWbNm+R6nlpa5uTnbF27HNsSWtTfW8r+j/6N7d2jUCL78Mof5\n+uWXYdo09UDWV1+ZPGZNe5JO+Pls2zZ1oy+z6Zxz4eeIfBhBu0M3OOPkRIKtLTVq1MjfILUMVahQ\ngd3/2Y24J3jz1zcJe3CXf/8bAgKy+SBWapMmqQJr48erTRE0LQ/phJ/PVqxQN2oze47q0PVDtL8K\nZe5G8mupUri6uiL0doaFhntTdz5r9RlxpeNoNbkVffokUL16LlZZmpnB0qXqJm6fPhASYspwNS0N\nnfDz0Z07auvT/v3Vw5cZOXT9EG+dtkKWLcuCv//W8/eF0Id9P8Tbxpur9lfp/Z8XGDtWsmcPHDyY\nww7LllX/MEqVUnWydTllLY/ohJ+PfvxR7Ys6bFjmbQIuHMDnVCLR3bpxNzpaz98XUpvHbqZaUjX8\nrf25fn8S9vYwY0YuOnR2Vg9k3bmj9jx4+NBUoWqagU74+SQxUS2/7tQJGjTIuM292Hu47Q2hdGwi\nf3l4AHqFTmFlZW7F3vf2Ym5tzpwrn9K6zWG2bVOLbnKsZUtVc+f4cbVGPzHRZPFqGuiEn2/8/FQZ\nleHDM2/zx40/GHYMohrUYu/jxwghaNq0af4FqWVLA/sGfPfid1AbNoZ3oHTpx7kb5YNaufPNN6pY\nz9tv65o7mkkZlfCFEN2FEGeFEBeEEBMyON9BCHFfCBGQ/Jpi+lCLtu+/h2rVMl97D3Bl56943ALL\n4SMIDAqibt26lMlo7aZWaAz1GIpPQx9E53ge2Q9n/XoICsplpyNHwkcfqZu5o0bpNfqayWSZ8IUQ\n5sB3QA+gMdBXCNE4g6b7pZRuya+PTRxnkXb1KmzdqgZslpaZt3NavZUYKzOsB71NUFCQnr8vAoQQ\n/NjrR2rb1cbqjVVQ+iIDB17PfccffaSe7Pr2W5g4USd9zSSMGeG3Ai5IKS9JKeOA1cBTxqnakxYu\nVCVU3nkn8zaJkRG0P/Q3JzrUJ9rcnIsXL+r5+yKifKnyrHltDZRJonS/XgQEOjBq1JrcdSoEzJoF\nvr7wxRfw4Yc66Wu5ZkzCdwBSD1luJB97UmshRJAQYpsQoolJoisGoqNhwQJ48UV42vNTt3/4kjLx\ncH9QX06dOoWUUo/wi5Dm1ZvzdbeveeRwCvOOk5g7tzpTpnyEzE2SFgK++07d+Jk1Sy3v0jdytVww\n1U3bE0ANKaUrMA/4LaNGQoihQohjQohjoaGhJrp04bZwIUREqE/lmUpKwmbRfzleHRp060dg8jZ4\nOuEXLcM9htPftT9J7T6H+pFMn36cd955h7hsV1dLxcxMJf2JE2HRIvWI9uPHpgtaK1GMSfg34f/b\nO/foqKp7j3/2zOSdkCcmAUECAQICIoSXhStYKY/aK9e2VKxKW7vQalF6L71otV2s4uXhpVxL1QoF\nFR/LB4oIFy5qFXkEAiKEN+GVCDGBwBCSkDBJZuZ3/9hDHpDHkMdMQvZnrb3OPrP3Oed7zpz9O3v/\nzj5706Xa+s2e3yoRkSIRueSJrwcClFJxV+9IRJaKSKqIpHasbVzgGwyHAxYuhDFj9ATYdbJuHTEn\n81g+ugPdo7uzb98+Is2QCm0OpRRL71nKwITbsfzkQaJ7PsPy5a8zfvx4CpoyraFS+lPeBQvg/fd1\n317zcZahEXhj8L8GeiqlkpRSgcD9wJrqGZRSCcrz/b9Saqhnv/bmFtvWWLEC8vLg2WfrySSCzJvH\nqRgr+feMQSlFRkYG/fv3N0MqtEFCAkJYff/HhIcEUTDuER58ZANpaWkMHz6cI00dNuE//1Mb/D17\nIDUVdu9uHtGGdkODBl9EnMBvgU+Bw8AHInJQKfWYUuoxT7afAAeUUnuBxcD90iTnZdvH6dQVsqFD\ndYWsTrZsQW3fzoLhLsan3IPD4eCbb75heL1NAkNrpmtkVz75+UpUzAneYx4ffPQ5BQUFDBkyhA8a\nNZZyNSZP1mM4KAUjR+r+vu27qBmuBxHxSxg8eLDcyLz1lgiIfPJJAxnHj5dLUWES/CySW5QrGzdu\nFEDWrl3rE52GluP5tW8Ks5HkWVPk21PfyogRIwSQJ598UsrKypq28zNnRMaN0zfZxIkieXnNI9rQ\n6gF2SSPtrvnStgUoL4c//xn699e9c+pkzx7YsIG3xsTQ75ZUEiMS2bRpE0opRo4c6TO9hpbh2Xse\n4nuX53E85F2e++olvvrqK5566ikWL17MsGHDOHToUON3Hh+vx9r+29/gyy/1aJvLl5svcw31Ygx+\nC/DSS3DsmHbp1Dv3+Lx5uCPCeabXaX7U60cAbNq0iYEDBxIVFeUbsYYW5X+fnkXogSd46+R/M3/b\nAl588UVWr15NTk4OgwcPZvHixbgba6SV0l/l7tkDKSn6y75hw2D79uY9CcONQ2ObBk0NN6pLJz9f\nJDJSZPz4BjKmp4uAZPz6HmE28k3uN+JwOCQ4OFhmzJjhE60G3/DmW05h0lRhNvLcF8+J2+2WvLw8\nmTBhggAyatQoOXLkSNMO4naLvPOOSKdO2s3zwx+KbN/ePCdgaFVgXDqthz/9SX9sVe+MdSLw7/8O\n8fG8MFLRKaITtyfczs6dO3E4HIwePdpXcg0+4MGfW/mh8zWsGb/m+S3PM+ufs4iPj2fdunUsW7aM\n/fv3M2DAAObMmUNZY/vYKwUPPACZmfD885CeDiNGwN1367H2nc7mPSlD26SxT4qmhhuxhr93r4jF\nIvLkkw1k/OADEZDyV1+R8LnhMm3NNBERmTNnjiilxG63t7xYg085c0YkNs4lcVMfF2YjD656UErL\nS0VEJC8vTyZPniyAJCcny5o1a8TtdjftgMXFIgsXVtX4O3cW+dOfRDIzm+FsDP6EJtTwjcFvJpxO\nkVGjRGJiROq115cviyQlifTvL59mrhdmI2szdY+cu+++WwYMGOAbwQafs2qVCLhl9B+fF2YjQ5YO\nkZzCnMr0DRs2SJ8+fQSQH/zgB5KRkdH0g1ZUiKxerX2MSukiP2iQyIIFIgcPaleQoU1hDH4r4IUX\n9NV8/XUvM37+uUxfP11Cng+R0vJSKSsrk9DQUJk+fbov5Br8xC9+oVuBz3+4WsLnhkviwkT57Phn\nlenl5eXy4osvSlRUlAAyZcoUOXbsWPMc/PRpkUWLRIYO1fcgiNxyi8hvfiOydq3IpUvNcxxDi2IM\nvp/JyBAJCBC5774GKkxHj4qEhorcc4+UOcskYWGC3PvuvSIikpaWJoB8+OGHvhFt8AuFhSK9e4vE\nxYl8ume/pLyUIsxGpq2ZJoWOwsp8Fy5ckGeeeUZCQ0PFarXKww8/LIcOHWo+IadOiSxZInLvvSJh\nYdoUBAWJjB2ra/9ff62brYZWhzH4fuTyZZF+/UQSEkTOnasnY0WFyLBhItHRIjk58vbet4XZyIZj\nG0REZO7cuQJIfn6+b4Qb/MaRI7on18CBIucKSmXmpzNFzVbS9X+6ynv73xOX21WZNy8vT2bMmCEh\nISGilJL77rtPNm3a1HQff3UcDpHPPxf53e9E+vSpqv1HRYlMmiSyeLFx/7QijMH3I089pa/i+vUN\nZJwzR2d8911xu90yZOkQ6f233pWFe+zYsXLrrbe2vGBDq2D9eu1SnzxZ29G0U2nS/5X+wmwkdWmq\nfHHyixpGPT8/X5577rlKV0+/fv3klVdekaKiouYXl5uru3g+8ohIt25VD4CEBJEHHhBZtkzk5Mnm\nP67BK4zB9xOvv66vYIO9cnbtErHZRKZMERGR7ae3C7ORl3e+LCIiubm5YrVaZdasWS0r2NCqWLBA\n3z8zZmij73Q55Y09b0iXRV0qDf+KjBVyueJy5TYlJSWyfPlyGTRokAASEREhTzzxhBw4cKDlhJ48\nqY38lCki8fFVD4CkJP1QeOcdM7SDDzEG3w9s3SoSGChy993aW1Mn586JJCfrbnEXLoiIyJQPp0iH\neR2kuKxYRETmzZsngGSaLnPtCrdbG3sQ+f3vqzwmlysuy8s7X5Y+L/URZiNxL8TJrz/5taw7uk4c\nFQ7Ptm5JT0+Xhx9+WIKCggSQgQMHyvz58yUrK6tlRR88qN08kyZpt8+VB0CfPiKPPiry9tsi337b\nchraOU0x+Epv73tSU1Nl165dfjl2c/Bf/6WHP96xA6Kj68h0+TJ8//v60/cvv4QRI/iu6Du6/bUb\n04dOZ9G4RYgIvXv3JiEhgc2bN/v0HAz+R0SPjvDKK/CHP+hvpq6Mii0ifJn1Jf/Y/Q/WHVvHpfJL\nhAWEMaTzEIZ1HkZqp1SSY5Lp4OrAmpVreP/990lPTwdg+PDh/OxnP+OnP/0pnTvXNkFdM+FyVd3f\nGzfCtm1QVKTTunaFUaPge9/Tw8b27w+BgS2npZ2glPpGRFIbta0x+I2nuBgiIupIdLv1ULarVsHK\nlfDjHwPw3JfPMXfLXI4/eZzu0d3ZvHkzd955J2+88QZTp071nXhDq8Ht1lPX/uMferl48bWT3Zc5\ny/gi6wvWH1vPju92sPfMXircFZXpHYI6EBMSQ6gKxVHowJ5jpzC/EByQEJ1Avx79GDZgGHcNuovk\nuGQ6R3TGarE2/8m4XLBvH2zZAlu36uWZMzotMBBuu02P5Z+aCoMHQ+/eEBzc/DpuYIzBb2243fDk\nk3pqur/8RQ+jAJy4cIIBrw5gfPJ4Ppr8EQBTp07l448/Ji8vj7CwMH+qNvgRt1vPYvjCC3r+hJUr\nISam7vwOp4OD+QfJvphN1sUsTheepsBRQGFZIRcdF7nouMj54vPYS+2Uce1wDTZs9I3tyx1JdzC0\n01BG3TKKHtE9mn/SHRHIzoZdu2qGK60AiwV69IC+fatCSgokJekLYCYBugZj8FsTDgdMnQoffAD/\n8R968mmlcLldjF4xmv1n93Pg8QPc3OFmCgsLSUxM5KGHHmLJkiX+Vm5oBaxYAdOmaW/Ie+/pSnBT\ncbldnM4/zdov1rJ+23rSM9O5aL0IiaA6KyRQ24COQR25q8ddjEkaw+huo+kV26tlZl1zu+H4ce0K\nOnSoKhw9WnPMn4gIbfi7datadu1aFTp2bJcPBGPwWwsFBTBpEmzerKtqM2dW3pB/2fYXZn4+kzcn\nvclDtz0EwJIlS3jsscfYsWMHQ4cO9adyQysiLU17A8+ehaefhj/+EYKCmm//IsLhw4dJS0tjy9Yt\nbNy/kRxrDnRDh3CdL8ISQf/Y/oxMHklq51R6xfYiOSaZsMAWaolWVOgHQWYmZGXpkJ1dFS8pqZk/\nKAi6dKn5EKi+3qUL3ICtZmPwWwMbN8Ijj0BOjq6mTZlSmXT43GFuX3I7E3pOYNXkVSilKC0tZdCg\nQQQEBLBv3z4zf62hBgUF8Lvf6VupXz/tGRw7tuUqtLm5uaSlpbFt+za2HdnGvqJ9ODo6IBHoSI2Z\nM2JsMfSI7kH/Tv1J6ZhCr9he9IztSffo7gTbWsgfLwJ2O5w+DadO6VA9fuqUnkD66rkFYmKufSB0\n6wbJydCzZz0v4VovxuD7k+JimDVLzy2anAxvvKF7JXjIKshiwjsTOF96noOPHyQ+PB6ARx99lKVL\nl/LZZ58xduxYP4k3tHbWrYMnnoBvv4U779S9w6rdXi2G2+0mKyuLPXv2sHP3TrYd28bR80c55z4H\nMUAselmtAq1QxNpi6R7ZnX6J/RjebTiDEgfR76Z+BNmasYlSFxUVkJtb8yFw9UOhsLDmNvHx2vBf\neQBcCcnJEB7e8pobgTH4/sDhgFdfhfnzIT9fV8fmzIHQ0MosO7/byY/e/RHlrnLW3L+GUbeMAmDl\nypVMnjyZWbNmMX/+fH+dgaGNUFYGy5bp2+vsWd3BZdo03Yj0tU0qLS0lMzOTw4cPc/jwYfYf20+m\nPZOc0hwuBV6qehjEAR4br9yKKGcUXWxd6BHRg55xPenbqS/9uvajV9deRIT6sJZdVKTdQ8eP62np\njh2riufl1cybkKCNf1ISJCZeG+LiICRE9z7ytukloudALS3V3bZLS2vGr17W8pt69dWWNfhKqfHA\nXwErsExE5l+VrjzpE4FS4Bcisru+fbZ5g//447pWf9ddMHeunlrOg9Pt5O19b/P4usdJCE9g/c/X\nkxKXAkB2djYDBw4kJSWFLVu2EHB1/zuDoQ5KSuC112DJEjh4ULunJ0yAf/1XmDgRYmP9q6+4uJjs\n7Gyys7M5cfIEe0/t5aD9IKcqTlEQXEB5dDnUZtudYHVaCXAFEEwwIdYQwmxhRARGEBkSSVRoFLHh\nscRFxBEbEUvHyI7cFHkTsaGxRAVHERUcRWRwZNPdSZcuwYkT1z4IsrN119KKitq3s1i04Q8N1cFm\n0y+fXa6ayyuGvjGVbItF7zskBHXuXMsZfKWUFTgKjAVygK+BKSJyqFqeicB0tMEfBvxVRIbVsrtK\n2rzBP3FCNxHHjKn8qbSilBUZK1i4fSEnC04y4uYRrL5/NTeF3YSI8PHHHzNz5kzsdjsZGRkkJSX5\n8QQMbRURPW3tihWwdq2umCqlezSOGAHDh2u/f58+0KGDv9VW4XK5OHLqCBlZGRzKPcTJ8yc5V3QO\n+yU7hY5CLjkv4XA7KFNlVKgK3XsoEN1SsDW8f4vbQoArgCAJIkSFEGYNI9QaSrgtnIiACCIDI4kM\njiQqOIrY0Fhiw2KJDY+lY4eOJEQmEB0eTVhYGMHBwViunoxaBC5c0Bf7SrDba62Ru8vLcVoVTuXG\nbbWALQCLzYYKDEKFhaFCw7CFRWANj6j5oLgSr+23gIDKVkRTXDpeXEaGAsdF5KTnYO8B9wKHquW5\nF3jT89lvulIqSimVKCJ51+7uxuBCp2iOBZ0n+8D77M7bzdbTW9mVu4tyVzlDOw9l4diFTOwxkeys\nbLYe3MqiRYtIS0ujb9++rFu3zhh7Q6NRCu64Q4e//x1274YNG/RHrh99pN0/V0hMhM6dr/VGxMfr\nh0F4uG4phIfrcMW2WK3N/4LYarVya9Kt3Jp0a4N5RYTS0lLsdjt2u52z586Sa8/lfNF5/X1BiZ2C\n0gIuXr5IYXkhxRXFlLpKuSyXcSgHdmXnrO0sBFMVGrJ2LsChgypTWMotWCosWN1WLDYLyqawWPVS\nWRTYwG1zIzbBFeTCFerCZdGhXtxAMdiKbQSpIIJVMMEqmBBLCMEWvQy1hRJq1SHMFkZYQFhlq6cp\neGPwOwOnq63noGvxDeXpDNywBv/HL/6Yr8q/0isuCLIHEXQmiMhTkZzOPc1vXL/Bbrfj9PQrjo+P\nZ8mSJfzqV7/CZvPmshsMDWOxVH24CrqTysmTNbu25+bql77p6XDunPf7ttl0CAioGa/+MKhrWV+a\nd9so9BvhMKBrw1qBDp5QHT2GjBu3W3CpUpy2AlwBF3HZLuIKLMAVUIAroBB3wEVcgReRwCLcgYW4\nAy8iQUW4Igpx2RzgtnmCFVw2cFqhLBAqQqE8BCpC9LI8FEt5COKJi8sCFhcoNygXWJxgcYOtDGdg\nCc6gEkqCSiCwBIIuQVApBJ7zxEsg8HKD5349+NTyKKWmAdMAunZt+E9szYxLGIdzq5MIZwQRrgiC\nLEFYIixY+luwDrRisViIiYkhJSWF3r17M2DAAEKrvdA1GFoCi0V3MElO1r79q6mo0C9+8/O1y/rq\nUFqq8zidOtQVhypX9NXL+tKuZ5vmQaFfPULVI+GW+jep8ISS+rPVwEJVS6Jeag5m5nYL4hKk1I2U\nXJvmFicVliKclmKclmLKLcV8w7jrEFYTb3z4I4DZIjLOs/4MgIjMq5ZnCfCViLzrWc8ERtfn0mnz\nPnyDwWDwA03x4VsazsLXQE+lVJJSKhC4H1hzVZ41wMNKMxwovJH99waDwdAWadClIyJOpdRvgU/R\nbaPXROSgUuoxT/qrwHp0D53j6G6Zv2w5yQaDwWBoDF758EVkPdqoV//t1WpxAZ5oXmkGg8FgaE68\ncekYDAaD4QbAGHyDwWBoJxiDbzAYDO0EY/ANBoOhnWAMvsFgMLQT/DY8slKqGMj0y8GbhzjgvL9F\nNAGj37+0Zf1tWTu0ff29RaRRg+r4c1CXzMZ+LdYaUErtMvr9h9HvP9qydrgx9Dd2W+PSMRgMhnaC\nMfgGg8HQTvCnwV/qx2M3B0a/fzH6/Udb1g7tWL/fXtoaDAaDwbcYl47BYDC0E3xm8JVSMUqpz5VS\nxzzL6FrydFFKbVRKHVJKHVRKPeUrfXWhlBqvlMpUSh1XSj1dS7pSSi32pO9TSg3yh87a8EL7zz2a\n9yultimlbvOHzrpoSH+1fEOUUk6l1E98qa8hvNGvlBqtlMrw3O+bfK2xPry4fyKVUmuVUns9+lvN\nKLlKqdeUUvlKqQN1pLfacgte6W9c2a0+w0pLBuAF4GlP/GlgQS15EoFBnngEevL0vr7SWIseK3AC\n6I6eTnnv1XrQw0L/H3pqneHADn/pbYT2O4BoT3xCa9Hurf5q+b5Ej+b6E3/rvs7rH4WeG7qrZ/0m\nf+u+Tv1/uFKOgY7ABSDQ39o9ev4FGAQcqCO9VZbb69DfqLLrS5fOvcAKT3wFMOnqDCKSJyK7PfFi\n4DB6blx/UTmBu4iUA1cmcK9O5QTuIpIORCmlEn0ttBYa1C4i20SkwLOaDtzsY4314c21B5gOfATk\n+1KcF3ij/wFglYicAhCR1nQO3ugXIELpCWjD0Qbf6VuZtSMim9F66qK1llugYf2NLbu+NPjxUjUL\n1hkgvr7MSqluwO3AjpaVVS91Tc5+vXn8wfXqegRd42ktNKhfKdUZ+Dfg7z7U5S3eXP9eQLRS6iul\n1DdKqYd9pq5hvNH/EtAHyAX2A0+JiNs38ppMay23jcHrstusX9oqpf4JJNSS9Gz1FRERpVSd3YOU\nUuHoWtsMESlqTo2Ga1FKjUHfNCP9reU6eRGYJSJuXclsc9iAwcD3gRBgu1IqXUSO+leW14wDMoC7\ngB7A50qpLabM+o7rLbvNavBF5O660pRSZ5VSiSKS52k61dp8VUoFoI39OyKyqjn1NYLvgC7V1m/2\n/Ha9efyBV7qUUgOAZcAEEbH7SJs3eKM/FXjPY+zjgIlKKaeIrPaNxHrxRn8OYBeREqBEKbUZuA39\n7srfeKP/l8B80Y7k40qpLCAF2OkbiU2itZZbr2lM2fWlS2cNMNUTnwp8cnUGjy9wOXBYRBb5UFtd\ntOUJ3BvUrpTqCqwCHmqFtcoG9YtIkoh0E5FuwIfA463E2IN3984nwEillE0pFQoMQ7+3ag14o/8U\nunWCUioe6A2c9KnKxtNay61XNLrs+vCtcyzwBXAM+CcQ4/m9E7DeEx+JfhG0D91UzAAm+voN+VW6\nJ6JrXCeAZz2/PQY85okr4GVP+n4g1Z96r1P7MqCg2rXe5W/N16P/qrxv0Ip66XirH/g9uqfOAbQL\n0++6r+P+6QR85rnvDwAP+ltzNe3vAnlABbol9UhbKbde6m9U2TVf2hoMBkM7wXxpazAYDO0EY/AN\nBoOhnWAMvsFgMLQTjME3GAyGdoIx+AaDwdBOMAbfYDAY2gnG4BsMBkM7wRh8g8FgaCf8P5vC9k+E\nVlG6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x118d24dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Load a precomputed posterior based on an analytic solution (see Lintusaari et al 2016)\n",
    "matdata = sio.loadmat('./resources/bdm.mat')\n",
    "x = matdata['likgrid'].reshape(-1)\n",
    "posterior_at_x = matdata['post'].reshape(-1)\n",
    "\n",
    "# Plot the reference\n",
    "plt.figure()\n",
    "plt.plot(x, posterior_at_x, c='k')\n",
    "\n",
    "# Plot the different curves\n",
    "for res, d_node, c in ([sim_res, 'd_sim', 'b'], [T1_res, 'd_T1', 'g'], [T2_res, 'd_T2', 'r']):\n",
    "    alphas = res.outputs['alpha']\n",
    "    dists = res.outputs[d_node]\n",
    "    # Use gaussian kde to make the curves look nice. Note that this tends to benefit the algorithm 1 \n",
    "    # a lot as it ususally has only a very few accepted samples with 100000 simulations\n",
    "    kde = ss.gaussian_kde(alphas[dists<=0])\n",
    "    plt.plot(x, kde(x), c=c)\n",
    "    \n",
    "plt.legend(['reference', 'algorithm 1', 'algorithm 2, T1\\n(eps=0)', 'algorithm 2, T2\\n(eps=0)'])\n",
    "plt.xlim([-.2, 1.2]);\n",
    "print('Results after 100000 simulations. Compare to figure 6(a) in Lintusaari et al. 2016.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interfacing with R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to run R scripts in command line for example with [Rscript](http://stat.ethz.ch/R-manual/R-devel/library/utils/html/Rscript.html). However, in Python it may be more convenient to use [rpy2](http://rpy2.readthedocs.io), which allows convenient access to the functionality of R from within Python. You can install it with `pip install rpy2`.\n",
    "\n",
    "Here we demonstrate how to calculate the summary statistics used in the ELFI tutorial (autocovariances) using R's `acf` function for the MA2 model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import rpy2.robjects as robj\n",
    "from rpy2.robjects import numpy2ri as np2ri\n",
    "\n",
    "# Converts numpy arrays automatically\n",
    "np2ri.activate()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    ".. Note:: See this issue_ if you get a `undefined symbol: PC` error in the import after installing rpy2 and you are using Anaconda.\n",
    "\n",
    ".. _issue: https://github.com/ContinuumIO/anaconda-issues/issues/152"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a Python function that wraps the R commands (please see the documentation of [rpy2](http://rpy2.readthedocs.io) for details):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "robj.r('''\n",
    "    # create a function `f`\n",
    "    f <- function(x, lag=1) {\n",
    "        ac = acf(x, plot=FALSE, type=\"covariance\", lag.max=lag, demean=FALSE)\n",
    "        ac[['acf']][lag+1]\n",
    "    }\n",
    "    ''')\n",
    "\n",
    "f = robj.globalenv['f']\n",
    "\n",
    "def autocovR(x, lag=1):\n",
    "    x = np.atleast_2d(x)\n",
    "    apply = robj.r['apply']\n",
    "    ans = apply(x, 1, f, lag=lag)\n",
    "    return np.atleast_1d(ans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  5.,  23.])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Test it\n",
    "autocovR(np.array([[1,2,3,4], [4,5,6,7]]), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load a ready made MA2 model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"104pt\" height=\"356pt\"\n",
       " viewBox=\"0.00 0.00 103.89 356.45\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 352.4503)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-352.4503 99.8937,-352.4503 99.8937,4 -4,4\"/>\n",
       "<!-- t2 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>t2</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"24.4469\" cy=\"-258.4503\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"24.4469\" y=\"-254.2503\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">t2</text>\n",
       "</g>\n",
       "<!-- MA2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>MA2</title>\n",
       "<ellipse fill=\"#cccccc\" stroke=\"#000000\" cx=\"47.4469\" cy=\"-175.672\" rx=\"28.5583\" ry=\"28.5583\"/>\n",
       "<text text-anchor=\"middle\" x=\"47.4469\" y=\"-171.472\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">MA2</text>\n",
       "</g>\n",
       "<!-- t2&#45;&gt;MA2 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>t2&#45;&gt;MA2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M29.3246,-240.8951C31.556,-232.8643 34.2919,-223.0173 36.9552,-213.432\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"40.3752,-214.1968 39.6802,-203.6248 33.6307,-212.3228 40.3752,-214.1968\"/>\n",
       "</g>\n",
       "<!-- t1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>t1</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"47.4469\" cy=\"-330.4503\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"47.4469\" y=\"-326.2503\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">t1</text>\n",
       "</g>\n",
       "<!-- t1&#45;&gt;t2 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>t1&#45;&gt;t2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M41.8792,-313.0211C39.2236,-304.7078 35.9868,-294.5751 33.0264,-285.308\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"36.3403,-284.1797 29.9633,-275.719 29.6723,-286.3098 36.3403,-284.1797\"/>\n",
       "</g>\n",
       "<!-- t1&#45;&gt;MA2 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>t1&#45;&gt;MA2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M49.1133,-312.4688C50.0021,-301.9749 51.0021,-288.4749 51.4469,-276.4503 52.1969,-256.1703 51.448,-233.5972 50.4036,-214.9639\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"53.8847,-214.5453 49.7822,-204.7769 46.8977,-214.9715 53.8847,-214.5453\"/>\n",
       "</g>\n",
       "<!-- S1 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>S1</title>\n",
       "<ellipse fill=\"#cccccc\" stroke=\"#000000\" cx=\"19.4469\" cy=\"-91.4469\" rx=\"19.3939\" ry=\"19.3939\"/>\n",
       "<text text-anchor=\"middle\" x=\"19.4469\" y=\"-87.2469\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">S1</text>\n",
       "</g>\n",
       "<!-- MA2&#45;&gt;S1 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>MA2&#45;&gt;S1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M38.2786,-148.0935C35.2657,-139.0306 31.9125,-128.9441 28.8803,-119.823\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"32.1044,-118.4264 25.6284,-110.0412 25.4619,-120.6348 32.1044,-118.4264\"/>\n",
       "</g>\n",
       "<!-- S2 -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>S2</title>\n",
       "<ellipse fill=\"#cccccc\" stroke=\"#000000\" cx=\"76.4469\" cy=\"-91.4469\" rx=\"19.3939\" ry=\"19.3939\"/>\n",
       "<text text-anchor=\"middle\" x=\"76.4469\" y=\"-87.2469\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">S2</text>\n",
       "</g>\n",
       "<!-- MA2&#45;&gt;S2 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>MA2&#45;&gt;S2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M56.9426,-148.0935C60.0966,-138.9332 63.6107,-128.727 66.7777,-119.5291\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"70.0982,-120.6359 70.0445,-110.0412 63.4796,-118.357 70.0982,-120.6359\"/>\n",
       "</g>\n",
       "<!-- d -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>d</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"47.4469\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"47.4469\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">d</text>\n",
       "</g>\n",
       "<!-- S1&#45;&gt;d -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>S1&#45;&gt;d</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M26.5124,-72.9131C29.7485,-64.4245 33.646,-54.2012 37.1961,-44.8887\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"40.5726,-45.8573 40.8644,-35.2664 34.0317,-43.3637 40.5726,-45.8573\"/>\n",
       "</g>\n",
       "<!-- S2&#45;&gt;d -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>S2&#45;&gt;d</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M69.2783,-73.2915C65.8425,-64.5897 61.6614,-54.0004 57.8813,-44.4269\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"61.0633,-42.9553 54.1353,-34.9395 54.5525,-45.5262 61.0633,-42.9553\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x118f72128>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ma2 = elfi.examples.ma2.get_model(seed_obs=4)\n",
    "elfi.draw(ma2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replace the summaries S1 and S2 with our R autocovariance function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Method: Rejection\n",
       "Number of samples: 100\n",
       "Number of simulations: 10000\n",
       "Threshold: 0.111\n",
       "Sample means: t1: 0.599, t2: 0.177"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Replace with R autocov\n",
    "S1 = elfi.Summary(autocovR, ma2['MA2'], 1)\n",
    "S2 = elfi.Summary(autocovR, ma2['MA2'], 2)\n",
    "ma2['S1'].become(S1)\n",
    "ma2['S2'].become(S2)\n",
    "\n",
    "# Run the inference\n",
    "rej = elfi.Rejection(ma2, 'd', batch_size=1000, seed=seed)\n",
    "rej.sample(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interfacing with MATLAB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a number of options for running MATLAB (or Octave) scripts from within Python. Here, evaluating the distance is demonstrated with a MATLAB function using the official [MATLAB Python cd API](http://www.mathworks.com/help/matlab/matlab-engine-for-python.html). (Tested with MATLAB 2016b.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matlab.engine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A MATLAB session needs to be started (and stopped) separately:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "eng = matlab.engine.start_matlab()  # takes a while..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly as with R, we have to write a piece of code to interface between MATLAB and Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def euclidean_M(x, y):\n",
    "    # MATLAB array initialized with Python's list\n",
    "    ddM = matlab.double((x-y).tolist())\n",
    "    \n",
    "    # euclidean distance\n",
    "    dM = eng.sqrt(eng.sum(eng.power(ddM, 2.0), 2))\n",
    "    \n",
    "    # Convert back to numpy array\n",
    "    d = np.atleast_1d(dM).reshape(-1)\n",
    "    return d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.41421356,  8.77496439,  1.        ])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Test it\n",
    "euclidean_M(np.array([[1,2,3], [6,7,8], [2,2,3]]), np.array([2,2,2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load a ready made MA2 model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"104pt\" height=\"356pt\"\n",
       " viewBox=\"0.00 0.00 103.89 356.45\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 352.4503)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-352.4503 99.8937,-352.4503 99.8937,4 -4,4\"/>\n",
       "<!-- t2 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>t2</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"24.4469\" cy=\"-258.4503\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"24.4469\" y=\"-254.2503\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">t2</text>\n",
       "</g>\n",
       "<!-- MA2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>MA2</title>\n",
       "<ellipse fill=\"#cccccc\" stroke=\"#000000\" cx=\"47.4469\" cy=\"-175.672\" rx=\"28.5583\" ry=\"28.5583\"/>\n",
       "<text text-anchor=\"middle\" x=\"47.4469\" y=\"-171.472\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">MA2</text>\n",
       "</g>\n",
       "<!-- t2&#45;&gt;MA2 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>t2&#45;&gt;MA2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M29.3246,-240.8951C31.556,-232.8643 34.2919,-223.0173 36.9552,-213.432\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"40.3752,-214.1968 39.6802,-203.6248 33.6307,-212.3228 40.3752,-214.1968\"/>\n",
       "</g>\n",
       "<!-- t1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>t1</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"47.4469\" cy=\"-330.4503\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"47.4469\" y=\"-326.2503\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">t1</text>\n",
       "</g>\n",
       "<!-- t1&#45;&gt;t2 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>t1&#45;&gt;t2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M41.8792,-313.0211C39.2236,-304.7078 35.9868,-294.5751 33.0264,-285.308\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"36.3403,-284.1797 29.9633,-275.719 29.6723,-286.3098 36.3403,-284.1797\"/>\n",
       "</g>\n",
       "<!-- t1&#45;&gt;MA2 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>t1&#45;&gt;MA2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M49.1133,-312.4688C50.0021,-301.9749 51.0021,-288.4749 51.4469,-276.4503 52.1969,-256.1703 51.448,-233.5972 50.4036,-214.9639\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"53.8847,-214.5453 49.7822,-204.7769 46.8977,-214.9715 53.8847,-214.5453\"/>\n",
       "</g>\n",
       "<!-- S1 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>S1</title>\n",
       "<ellipse fill=\"#cccccc\" stroke=\"#000000\" cx=\"19.4469\" cy=\"-91.4469\" rx=\"19.3939\" ry=\"19.3939\"/>\n",
       "<text text-anchor=\"middle\" x=\"19.4469\" y=\"-87.2469\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">S1</text>\n",
       "</g>\n",
       "<!-- MA2&#45;&gt;S1 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>MA2&#45;&gt;S1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M38.2786,-148.0935C35.2657,-139.0306 31.9125,-128.9441 28.8803,-119.823\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"32.1044,-118.4264 25.6284,-110.0412 25.4619,-120.6348 32.1044,-118.4264\"/>\n",
       "</g>\n",
       "<!-- S2 -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>S2</title>\n",
       "<ellipse fill=\"#cccccc\" stroke=\"#000000\" cx=\"76.4469\" cy=\"-91.4469\" rx=\"19.3939\" ry=\"19.3939\"/>\n",
       "<text text-anchor=\"middle\" x=\"76.4469\" y=\"-87.2469\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">S2</text>\n",
       "</g>\n",
       "<!-- MA2&#45;&gt;S2 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>MA2&#45;&gt;S2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M56.9426,-148.0935C60.0966,-138.9332 63.6107,-128.727 66.7777,-119.5291\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"70.0982,-120.6359 70.0445,-110.0412 63.4796,-118.357 70.0982,-120.6359\"/>\n",
       "</g>\n",
       "<!-- d -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>d</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"47.4469\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"47.4469\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">d</text>\n",
       "</g>\n",
       "<!-- S1&#45;&gt;d -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>S1&#45;&gt;d</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M26.5124,-72.9131C29.7485,-64.4245 33.646,-54.2012 37.1961,-44.8887\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"40.5726,-45.8573 40.8644,-35.2664 34.0317,-43.3637 40.5726,-45.8573\"/>\n",
       "</g>\n",
       "<!-- S2&#45;&gt;d -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>S2&#45;&gt;d</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M69.2783,-73.2915C65.8425,-64.5897 61.6614,-54.0004 57.8813,-44.4269\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"61.0633,-42.9553 54.1353,-34.9395 54.5525,-45.5262 61.0633,-42.9553\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x118f82908>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ma2M = elfi.examples.ma2.get_model(seed_obs=4)\n",
    "elfi.draw(ma2M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replace the summaries S1 and S2 with our R autocovariance function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Method: Rejection\n",
       "Number of samples: 100\n",
       "Number of simulations: 10000\n",
       "Threshold: 0.113\n",
       "Sample means: t1: 0.602, t2: 0.178"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Replace with Matlab distance implementation\n",
    "d = elfi.Distance(euclidean_M, ma2M['S1'], ma2M['S2'])\n",
    "ma2M['d'].become(d)\n",
    "\n",
    "# Run the inference\n",
    "rej = elfi.Rejection(ma2M, 'd', batch_size=1000, seed=seed)\n",
    "rej.sample(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, don't forget to quit the MATLAB session:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "eng.quit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verdict\n",
    "\n",
    "We showed here a few examples of how to incorporate non Python operations to ELFI models. There are multiple other ways to achieve the same results and even make the wrapping more efficient.\n",
    "\n",
    "Wrapping often introduces some overhead to the evaluation of the generative model. In many cases however this is not an issue since the operations are usually expensive by themselves making the added overhead insignificant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### References\n",
    "- [1] Jarno Lintusaari, Michael U. Gutmann, Ritabrata Dutta, Samuel Kaski, Jukka Corander; Fundamentals and Recent Developments in Approximate Bayesian Computation. Syst Biol 2017; 66 (1): e66-e82. doi: 10.1093/sysbio/syw077\n",
    "- [2] Tanaka, Mark M., et al. \"Using approximate Bayesian computation to estimate\n",
    "tuberculosis transmission parameters from genotype data.\"\n",
    "Genetics 173.3 (2006): 1511-1520.\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}