RunningMatlabInJupyter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running Matlab in Jupyter\n",
    "\n",
    "Every now and then I have to use some Matlab functions in a project. For this reason I still have my Matlab version 2012b around, from a time when I was working much more in Matlab and I still had a license therefore (the one good thing of Matlab: a once bought license is forever, you just cannot update). What I would like to do is running this particular Matlab code in Jupyter, doing my usual Python stuff before and after calling that function.\n",
    "\n",
    "There are a couple of packages out there to help using Matlab in Jupyter. There are even entire Matlab kernels, so you can use Matlab in your dedicated Jupyter notebook. Some repositories:\n",
    "\n",
    "- https://github.com/imatlab/imatlab;\n",
    "- https://github.com/Calysto/matlab_kernel;\n",
    "- https://github.com/arokem/python-matlab-bridge.\n",
    "\n",
    "Using the second package from this list, `matlab_kernel`, there was even a Matlab class created entirely in Jupyter Notebooks, available on ReadTheDocs and on Github, something we are familiar from Python: Harvard AM111 2017 Fall; *Introduction to Scientific Computing*, https://am111.readthedocs.io.\n",
    "\n",
    "**However, most of them are overshoot for my purposes, or require MATLAB 2016b or newer with the new [MATLAB API for Python](https://www.mathworks.com/help/matlab/matlab-engine-for-python.html).**\n",
    "\n",
    "This notebook shows two (of probably many) ways of running Matlab code in a Notebook without too much fiddling. The first one is simple, using the `%%bash`-magic of the notebook. Downside of the simple approach is that there is no realtime output of the Matlab shell. The second option uses `subprocess`, with which we can monitor the output from Matlab. And at the end an example to pass Python variables to Matlab and the result from Matlab back to Python.\n",
    "\n",
    "For this code to work you must have Matlab installed, and `matlab` must be in your bash-PATH, if not you have to replace the `matlab`-commands with the full path to your binary. They work on my Linux machine, and I assume it is the same on a Mac. It might work on Windows as well, but I don't know."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import numpy as np\n",
    "from scipy import io as sio\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create script-string for Matlab\n",
    "\n",
    "We basically just write the whole script in one string, on one line. In order to make it visually more attractive we split it over several lines, ending each line with `\\` to continue the line. This is rather fragile, and you really have to end each command with a `;` as the command is all in one line, otherwise it will fail. So it is only convenient for small scripts. (There might be better ways of which I am not aware of.)\n",
    "\n",
    "The example script doesn't make a lot of sense, and is just for demonstration purposes. I specifically entered the `pause(3)`-commands, so we can check if the output in the notebook is real-time or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "script = \"\"\"\\\n",
    "    x = 1:.5:pi; \\\n",
    "    disp(['x = ', num2str(x)]); \\\n",
    "    pause(3); \\\n",
    "    ; \\\n",
    "    y = sin(x); \\\n",
    "    disp(['y = sin(x) = ', num2str(y)]); \\\n",
    "    pause(3); \\\n",
    "    ; \\\n",
    "    disp('Finished demo in Matlab'); \\\n",
    "    quit; \\\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 1: `%%bash-magic`\n",
    "\n",
    "This method is straight forward, **but you'll have to wait till the end to see the output**.\n",
    "\n",
    "The `-nojvm` tells Matlab to start without Java, which implies that no GUI is started an the splash-screen is also suppressed. You can provide as many arguments as you want with \"\" in the first line, they are numbered sequentially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                            < M A T L A B (R) >\n",
      "                  Copyright 1984-2012 The MathWorks, Inc.\n",
      "                    R2012b (8.0.0.783) 64-bit (glnxa64)\n",
      "                              August 22, 2012\n",
      "\n",
      " \n",
      "To get started, type one of these: helpwin, helpdesk, or demo.\n",
      "For product information, visit www.mathworks.com.\n",
      " \n",
      "x = 1         1.5           2         2.5           3\n",
      "y = sin(x) = 0.84147     0.99749      0.9093     0.59847     0.14112\n",
      "Finished demo in Matlab\n"
     ]
    }
   ],
   "source": [
    "%%bash -s \"$script\"\n",
    "matlab -nojvm -r \"$1\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 2: Using `subprocess`\n",
    "\n",
    "By using `subprocess` and writing a small while-loop we gain real-time output from Matlab.\n",
    "\n",
    "At the same time we can suppress the initial few lines (`< M A T L A B (R) > ...`), as we don't want to see them every time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matlab(command, startup_info=False):\n",
    "    \"\"\"Run a command in Matlab, flushing stdout continuously.\n",
    "    \n",
    "    startup_info: bool\n",
    "        If True, it will print the initial Matlab copyright etc.\n",
    "    \"\"\"\n",
    "    \n",
    "    # We add `disp('outputstarting')` at the beginning, so we can check for\n",
    "    # it and suppress the MATLAB-startup message.\n",
    "    if not startup_info:\n",
    "        command = \"disp('outputstarting'); \"+command\n",
    "    \n",
    "    # Start the process\n",
    "    process = subprocess.Popen(['matlab', '-nojvm', '-r', command],\n",
    "                               stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "\n",
    "    # Monitor STDOUT and STDERR\n",
    "    while True:\n",
    "        line = process.stdout.readline()  # Get line\n",
    "\n",
    "        # If line is empty, Matlab finished, exit loop\n",
    "        if not line:\n",
    "            break\n",
    "            \n",
    "        # Print line depeding on startup_info\n",
    "        lline = line.rstrip().decode('latin')\n",
    "        if startup_info:\n",
    "            print(lline)\n",
    "        elif 'outputstarting' in lline:\n",
    "            startup_info = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can simply call `matlab(script)` to run our desired script in Matlab:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = 1         1.5           2         2.5           3\n",
      "y = sin(x) = 0.84147     0.99749      0.9093     0.59847     0.14112\n",
      "Finished demo in Matlab\n"
     ]
    }
   ],
   "source": [
    "matlab(script)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Passing arguments\n",
    "\n",
    "The reason why I want to run Matlab in the notebook is because I do some calculations before and after calling the particular Matlab script, so passing variables is crucial. Thanks to `scipy.io` this is rather easy. To make it easier we also write a little function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save4mat(variables, fname):\n",
    "    \"\"\"Store variables in a dict; save dict as `fname.mat` in current directory.\n",
    "    \n",
    "    variables : list of strings\n",
    "    fname : string\n",
    "    \"\"\"\n",
    "    pvars = {}\n",
    "    for name in variables:\n",
    "        pvars[name] = eval(name)\n",
    "    \n",
    "    sio.savemat(fname, pvars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do some calculation in Python, and store it in `pythonvars.mat`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 3*np.pi, 1001)\n",
    "mult = np.linspace(0, 1, x.size)\n",
    "\n",
    "save4mat(['x', 'mult'], 'pythonvars')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Call Matlab; load `pythonvars.mat`; do some calculation; store in `output.mat`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished calculation in Matlab\n"
     ]
    }
   ],
   "source": [
    "matlab(\"\"\"\\\n",
    "    load pythonvars; \\\n",
    "    y = sin(x.*mult); \\\n",
    "    save('output', 'y'); \\\n",
    "    disp('Finished calculation in Matlab'); \\\n",
    "    quit; \\\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load `output.mat` in Python, get `y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = sio.loadmat('output')['y'][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A small figure to check if this did what we wanted..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4FVX6xz9vekiFJHRCQgg91FAEBPWHig2sFAs2xN5ZV3RVll3rsmtb14aAlSIWcEFdEEEpgYReQyItIQFSSAglpNzz+2MGvMSUC7n3zk1yPs8zz52Zc86c78ydmXdOfUUphUaj0Wg0tcXLagEajUajqR9og6LRaDQap6ANikaj0WicgjYoGo1Go3EK2qBoNBqNxilog6LRaDQap6ANSj1CRDqJSJnVOjwJEblPRJY44TgHRWSwMzQ5kNdsEflLFWF3i8h37tChqRoRWSoio63WYc+53Osi8oqITHO2hgZnUETkmN1iE5GTdtu3WK3PmYhIkojcarWOhoSIBIiIEpHWrji+UuojpdQ155PWfIkoEbnGbl+wua+581SeOfZw8xk7JiJFIrKjvtyPSqlLlFJzziet+XFyUkTCKuzf4eh/4akfjw3OoCilgk8vwH7gGrt9n7tDg4j4uCMfjaYS8oG/iYi4Kb/d5rMWCrwAzBSRODfl7cnsB0ad3hCRftSD93GdPwFnIiIhIlIsIqHm9t9F5JSIBJrbU0XkFXO9iYh8ISI5IrJHRJ6q6iE1vwy/EJE5IlIEjBERbxF5TkR2i0iuiHwuIuFm/CCz2iNfRApEZI2INDbDzqp6qaroKiL/BPoC08wvxH+aef7b1FwoIptEpKOD12aoiKw30yWJSF+7sCQRecH8PSoii07rNcMvMsMKRWS/iNxs7g8UkTdEJMM8r7dFxL8aDQ+IyE7za3eLiCSY+583/4MiEdkqIldVc4weZnXFETPPJ839Z1UzmV/X6VUcY5D5nxSKSJaIvG73kfCL+ZtqXvdrzTTXichm8//8VUS62B2vn/lfFInIZ4BfNfrPVGvI76Whe0TkN/OcXq8qrckCwB+7l1mF459Vqq0iv/vM/I6KyF9EpKOIrDWvx+dSyQeTMpgLnAQ6i8hPInJPhbx3icjwGvRj3sd/Ne+lQyIyXURCzLBOIlImIneKSKZ5r//JLq2Pea/uNvUni1kiEJFudvfGjtP/XU3XybweP4nIW+b/+5uIDKvhND4FxtltjwM+qZDHdeZ9cVRE9onIM3bBvwDe8nvtSq9KNL5rXoOj5v8zoEKUIBH5yrzvkkWkaw2aa0Yp1WAXYC8wrMK+tcBV5vovwG/AxXZhV5jrc4EvgWCgPbAHuKWKfF4BTgFXYhjxQOBp4FegJRAAzARmmPEfBeaZ8XwwDEOQGXYQGFzh2NPM9U5AmV1YEnCr3fZIYDXG16IX0BVo6sB1agocxXgJ+QB3ADlAmF0+qUAcEASsAiabYe2BY8ANZtoooIcZ9p55nuFAGPAj8EIVGm4D9gG9AAE6Aq3NsNFAC/OcbgOKgEgz7D5gibne2NT9EMZLNRToa4bNBv5il99wIN1u+8x1B/qZ/4m3ec7pwH1mWACgTmsz9w0AsoE+ZpoJwC7zegQAWcADgC9wC1Bmr6XCdbA/n9N5fW2eSyxQAFxUzX04zfwfd5pags1jNK/inqksv9P3fS+gFPgf0BZoAqQBoyteQ/O/GQOUADEYL9Dldvn0N6+xtwP34wPADjPPUOC/wId2z4AC3jH19jXzbGeGPwdswLgvvcxzCDePk21ef28zXT7QvgoNZ66TeY1KzXPyBh4H9laj/yBwIcY7o535v2eb6/b/xf9hPKNeQG9Tz/DKnvWK/5W5PQ7jnvcFngUyAF+7e6EEGGGG/wXjGa7x+le36BLKH1kODDW/lOOBd83tEKA7sNIMuwH4s1LqmFIqHXgD42VW5XGVUouUUjal1EngXuBppVSWUqoY+CswWkQE4+aMAuKUUmVKqWSl1HEnnFspxoPTCeOjcZtS6rAD6UYCG5VSc009M4FM4Aq7OB8qpX4zdc4Depr7bwO+U0p9ZabNUUptMr9i7wIeVUoVKKUKMW7yMVVoGA+8pJTaoAxSlVKZGCcyRymVbV7bT4EDGC/vilyL8YL7t1LqlFLqqFIq2YHzPwul1FrzPylXSv2G8ZIeWk2Se4F/K6XWmWk+wDBofYAhwCml1H+UUqXKqHbdfI6SXjLPZQ/GR1DPGuJ/CRRz9hfyufCKed9vwDCMC5VS+5RS+RjGxf5rOVZECoBc4ClgrFJqL/AV0EtEos14twFfKKXKHcj/FuAfZp5HMV6Wt5jPzmleUEoVm//vToxnF4z76GmlVLp5v2xQShUA1wFblVKfm/9RMvAdxnPuCKlKqU9M/R8DbcWscagCBXxunvdVQApw1rOolPrJfEZtSqn1GB+x1d1nZ2dg6DmilCoFXgIiMIzWaVYppRaY4a8AkRiG67zRdfl/ZDnwPMYXUwqwFPgnxhfJFqXUURFpi/HVsN8u3T6gVTXHzTi9Yt74bYBFImI/O6cXxp/+EdAcmCciwRhF4eccfNiq43sMY/I+0EpE5gFPKaWO1ZCuJcb52VPxfA/arZ/A+IIF4zx/q+KYvsA2u/eAYHydV0ZVx0FE7sYo1Z1+OQVjPBwOH+NcMKur/onx8J0uRa6sJklbYJR91QtGtVYrDAOfWSF+xWtdE1Vd+0pRSikReQ54C5h/jnkBHLJbP1nJtn3+e5RS7SvRcFxEvsYwBFMxSpmXOZh/xftxH8b/0MTcLldK5dqFnwCCzeeuFZXfA22BIabxO40PcMRBTRX/AzCuQ0ElcU/zCbAQ6EGF6i4wqlYxDEEXjPvFH6OqzCFEZBJGbUJzDAMWgPFcpJpRzryTlFJlIpKFcW3PG11C+SO/YvzBV2EYl40YL+HLzG0wbh4bv7/AMNcPVHPcM4ZDGWXOA8AlSqlwuyVAKZVrfj0/r5TqhPEFexO/f7kfBxrZHbe6HiFnTSVtftn/SynVC+OLrQfGi7gmsjAeOHtqOt/TZGBUC1UkG8N4xNmdf5hSKuJcjiMiHYC3MaqRmiilwjGqoCprz6pKC5zbdf0QWG9qDwWm2OVX2fTdGcDzFf7rRkqprzGuQ8UeYdF/PIRzUUp9h3Efj68QdC7XobZ8DNyKUTV2yCzxOELF+zEaw5DlV5fI7rmr7B7IAP5X4T8KVko95qCmc0YptQvIAy7CaNuqyFxgDtBGKRWGUS1e3X12BhG5FHgYo+QVjmFsT3L2c9HGLr43hjHJOvcz+R1tUCpgVr1sA+7HqKayYZRUxmMaFKXUKeAb4CUxGtDjMF7Mn51DVu8Br4hIGwARaSpmd04RGSYiXUTEC6Ptogw4XTrZCIw1GxcHYFRHVcUh7Iq4IjJARBLN6qbjGHWo5WbYfSKys4rjLMConrjRzHccxkP8gwPn+QlwtdnA6C0iUSLS3SxmTwfeFJFIMWhjPgiVMQ14WoxGdRGRDmJ0zQ3GMO45gJeI3IdRP14Z3wLtReR+EfETkVD5vXPBRlNnuIi0wngYqyIEKFRKHTMbMs80Lpv3RiFnVy18ADxsXnsRo6vuCBFphFFFFWBefx8RGcvv1TOu5i/Anyvs2wjcKEYDfCeML1xXsQzj/3uRPzZIJ4nI01WkmwVMFJFosyr67xjVZY744piG8dy2M/+LXmbV1LcY9/hoEfE1748B5geLK7kN+D/zvjmDWZoKBvKUUsUiMhDjw/I0hzEa5av6+AjBqOLOwSjdTMEoodgzUESuFhFfjOrIPIwPpfNGG5TKWY5hydfbbQcBK+zi3Gv+7sOoFpuGUSfqKK8BS4ClYvT8WsXv9ZetMKoiioCtwCKMrxWAZ4AEjKL0JIzG5Kp4HRgnRq+V1zC+VGaaaXeb2t8y47ahimobpdQhjMa7ZzFuuoeAq82652ox2xhGmrqPYBjn071JHsP4IkrBeAn/QBXGwGwb+RdG+8xR8zfcrFt+zzxGNkbDdEoVxzgCXIpR2juMUfQ/3WNuOkbJZj9GI++sak7rcWC8iBzDaPytOB7heeBLMXr8jFBKrQQewahqLMBod7jZkKROYnxFPmBen6sw6u5djlLqJ2BThd2vYVT15GAYwnP5SDrX/BVGFU5X4IsKwa0xnonKeBejI8IqjOqrfOAJB7N9BaOaaSnGffQe4G/eG5cDd2LcR1kYhsrXweOeF0qptMpKZua1uQ+Yar4fnsJo+zodfgTjv1pn3mcV282+4/dORbsx2rByKsT5CqMd8whGW9ENp6vVxejt5ug1PYM4ZtQ19R0RWQbcbRoAjcYtiMgEYJRSapjdvvYYnTwutk6Z5nzQBkWj0ViCiJwu9b+sjDEqmjqOrvLSaDRuR0RGYFQ7pmNUX2rqAbqEotFoNBqnoEsoGo1Go3EKDWpgY2RkpIqJibFahkaj0dQp1q1bl6uUiqopXoMyKDExMaSkVNqjVKPRaDRVICIOzd6gq7w0Go1G4xS0QdFoNBqNU9AGRaPRaDROoUG1oVRGaWkpmZmZFBcXWy3FowgICKB169b4+rp05gmNRlOPaPAGJTMzk5CQEGJiYhC3eUX1bJRS5OXlkZmZSWxsrNVyNBpNHcHSKi8xXHceFpGtVYSLGG4108Vwn9rbLux2EUkzl9vPV0NxcTERERHamNghIkREROhSm0ajOSesbkOZieELoSquwPCaGI/h7+JdMPy5Ay9gOMHqB7wgdj7MzxVtTP6IviYajeZcsbTKSyn1i4jEVBNlJPCJOZVzkumrogWGQ5rFpstRRGQxhmGqbspxjaZKjuRkc2BXCsezdmE7WUCZDXa1v4uoEH+6hpUS07oV3j4NvoZY4wSUzUZq8hIKdq0AnwBaDLiJtrHxVstyCp7+hLTCzk0lhqvUVtXs/wPm9NgTAKKjXe4I77wQEZ544gn++c9/AjB16lSOHTvG5MmTa33syZMn8+GHHxIVFUVZWRkvvfQSI0aMqPVx6wPph4v4ev0Bumz4K1efWoR9ETfDFsVtOwcA8Invy0R57yY1bBD+PUeTMPR6xMvqwr2mLnIo8zdyP7mTriWGG5oy5cVFm1tycf9TPHtVZwJ8vS1WWDs83aBUVu+iqtn/x51KfYDhKIjExESPnAnT39+fr7/+mkmTJhEZWZkr9Nrx+OOPM3HiRHbs2MGFF17I4cOH8WqgL0Rls7Fp6Rx81vyHW4oe4rhXCA8168nqFm0JbtefqJhuhEU0p2VAIzaeKufg0WKOJB8i9bf/0aHgF8KW/4+9v04mp+9TJA6/TVcNahwm62AWfDScGNtR1nR+mk6X38OJ4lMMSy5g5qo9xO2fy9gJT+Mf0Kjmg3konv5WycTO7zGGF7esavbXSXx8fJgwYQKvv/76H8LuuOMO5s37fXbv4OBgAJYtW8bQoUMZNWoUHTp04Omnn+bzzz+nX79+JCQk8Ntvf/ST1blzZ3x8fMjIyCA2NpbS0lIAjh49SkxMzJnt+sq2lQv57cVEeq64jyZlh3j+wmBWT/o/Hnv0z1xwxyskDLmO5tHxBAaF4O3tTXgjPzo1D+WCa+6i72OzCZy0m+ReL6PwYuGKZEZ/kMRvOcesPi1NHeBkSTm3f76LBbYLOTByLv3HTCKscSQtWrRi8oiufHRRKXfkv8nGaQ9aLbVWeHoJZQHwkIjMxmiAL1RKZYvIjxh+oU/XUlyG4Q63Vvz1u21szzpa28OcRZeWobxwTdca4z344IN0796dp556yuFjb9q0iR07dtCkSRPatWvH+PHjWbt2LW+++SZvv/02b7zxxlnx16xZg5eXF9HR0Vx00UUsXLiQa6+9ltmzZ3PDDTfU2zEnh44cJWP6HSQW/USWNCW554v0vPIebvDzP6fj+PkH0HfkA5RfNYGOyfuZvzidt9/6Bzf3bEzf6x7R1WCaKnl1QQppOcfpcterdOjwxzkW/2/4dSTt+ZEBh2axaekV9LhklAUqa4/V3YZnAauBjiKSKSJ3i8h9InKfGWURhj/kdOBDDL/bmI3xfwOSzWXK6Qb6ukpoaCjjxo3jrbfeqjmySd++fWnRogX+/v7ExcVx2WWXAZCQkMDevXvPxHv99dfp2bMnEydOZM6cOYgI48ePZ8aMGQDMmDGDO++806nn4yks2pLNZW8lkVVYzOo299DkTxvoe+1D+J6jMbHH28eHsRe044fHhnBzcDL9trzA2nfupKy0xInKNfWFHWt+5MEtN/JsomJIJcbkNL3veoN9Xq2J/PUvFJ887kaFzsPqXl5jawhXQKVlQKXUdGC6M/U4UpJwJY899hi9e/c+6+Xu4+ODzWYDjAGHJSW/v7T8/X9/KXp5eZ3Z9vLyoqys7EzY6TYUewYNGsTevXtZvnw55eXldOvWzSXnZBXHiwrYPP0Rph68kJhWnek2ag7tmoY4NY9moQFETfyO1dMe4YLsz9j0ryw6PvwNAY2CnZqPpu5SXlaG7/8mYRNfbhl+YbVx/fwDOHrRSyQsHcfqL1/hgnF/c5NK56HL6B5EkyZNGDVqFB999NGZfTExMaxbtw6A+fPnO7WdY9y4cYwdO7belU6y9uzk8OtD6Z+/gD93zmXe/QOdbkxO4+XtzQX3vsOaLn8h4UQyaW9eQ3HxSZfkpal7bPxxJu3LfyOjz59pFBxWY/yEISP5b8goZu5tzImSshrjexraoHgYTz75JLm5uWe277nnHpYvX06/fv1Ys2YNQUFBTsvrlltu4ciRI4wdW21BsU6xdeV3BH48jAhbDtsumc7l4ybh6+3627z/qD+R0vNvJB1vwX1fbKa03ObyPDWejbLZaLzubfZ5tab3leMdTtf8hlf58URHvliz34XqXEOD8imfmJioKjrY2rFjB507d7ZIkbXMmzeP+fPn8+mnn1YaXteuza8/zmPAqglkebfE6+ZZtGmf4HYNX6zZzzPfbOHuPmH85YaBuqG+AbP+l+/ovfRWknu+SN9rHzqntPe+8x09j/zIhKff8IgBtSKyTimVWFM865VqLOHhhx/m+++/Z9GiRVZLcQrTft3NP372Zmrk9Qwd/yqh4RGW6Li5fzQF2bu5ccMtrDl1JwNunWyJDo31vLo9gki/v/HGlfecc9oJsYfpk/Mpm34ZUqd6fOnPpwbK22+/TXp6Oh06dLBaSq1QNhtLZvyVtxam8H8J0Vz22PuWGZPT3HfNhWSGdKdf2htsXbHAUi0aa0g/fIw1e4/QbeBV59WjMOH/biGXcFTyRzVH9iC0QdHUWZTNxpr37mPYvn/xcrvNvD22N/4+1k9d4eXtTcd7PyXDuxXNlzxM3qFMqyVp3EzWN8/zgu8n3Ni70hmhasTPP4C0liPpdmJtnbp/tEHR1EmUzcaa9x9gwOE5JDUdxZXjp+Dt5TnToASFhGO7YTrB6jhZ02/DVl5utSSNmyg+eZzu2XPoHnaSqNCA8z5O88G34SM20pd95kR1rkUbFE2dQ9lsrPnwEQYcmsWayOvpf9/7Htn4Hdu1P5u6TSL7hDB71U6r5WjcxLaf5xLOMfz61a47fmyXvuzybk/mvt1OUuZ6PO8p1Ghq4D/fJ9M663vWRFxLvwc+8khjcpp+NzzOZzEv8/fFGWTkn7BajsYNyLavyCWcLgOvqfWxlgz6gifzR9aZe8dzn8QGhIhw2223ndkuKysjKiqKq6++utp0GzduPKuX1syZM3nooeq7J06ePJmpU6fWTrCFfLRiD//4NZcZXabT94HpHm1MAMTLi5dv6E6MHGTT9EdQNj0+pT5TVJhP12NJpEdd6pTuvtf0NObA/WFT3RiT4tlPYwMhKCiIrVu3cvKkMcJ68eLFtGpVc2NeRYNS30n574f4/TCRK7pE8cxNQ/Dytr4B3hFahQfyQrdcrj72JSkL3rVajsaF/Lp9H3PLhxLW/1anHK9Nk0Z8GPoRA1c5PjDSSrRB8RCuuOIKFi5cCMCsWbPOGr2+du1aBg4cSK9evRg4cCCpqamUlJTw/PPPM2fOHHr27MmcOXPOOt53331H//796dWrF8OGDePQoUNnwjZt2sQll1xCfHw8H374oXtOsJZsXv4N3ZP/TJ9Gh3j9pq4e1QDvCH2vf5xUn47EbnyNosI6PY+pphrmpZbxXvCDdOoz1GnHDGkaQ6eSrRTkHnTaMV2FNigVmXHVH5e15ku35ETl4Rs+N8KP5/0xzEHGjBnD7NmzKS4uZvPmzfTv3/9MWKdOnfjll1/YsGEDU6ZM4ZlnnsHPz48pU6YwevRoNm7cyOjRo8863uDBg0lKSmLDhg2MGTOG11577UzY5s2bWbhwIatXr2bKlClkZXm2K5ld65fTfum9ZHq3odX93xIQ6LzpZ9yFl7c3cuU/aKIK2fbFM1bL0biAwiP5FKat5uqEZk51vBbR51q8RZG28munHdNV6JHyHkL37t3Zu3cvs2bN4sorrzwrrLCwkNtvv520tDRExKEJIjMzMxk9ejTZ2dmUlJQQGxt7JmzkyJEEBgYSGBjIxRdfzNq1a7n22mudfk7OIHtfKhELbqNQQgkZP5+wxs73aOkuOvQeytpfr6TPwbnsS72Pth17Wi1J40TSVs7jK9/n2NmsG+C8mcvjug8iZ35jvNN+wPTg4bFog1KROxdWHebXqPrwoIjqw2tgxIgRTJw4kWXLlpGXl3dm/3PPPcfFF1/MN998w969e7noootqPNbDDz/ME088wYgRI1i2bNlZ/ukrfj15qhvbouJS3przA4/hTcmYubRoGWO1pFoTN+Y1ZvzHn+2rjvBGR6vVaJyJ7PqBfELp0Mt51V1glG73NBlMl7wllJaW4Ovr59TjOxNd5eVB3HXXXTz//PMkJJw9qWFhYeGZRvqZM2ee2R8SEkJRUVGlx7JP8/HHH58VNn/+fIqLi8nLy2PZsmX07dvXiWfhHMptikdmbWBuQTxpY36lbafeVktyChHNWlM85Bm+TT3Bhv1HrJajcRJlpSXEH03it/BBLuksUtb9Zl4pG8OW/bk1R7YQqz02DheRVBFJF5GnKwl/XUQ2mssuESmwCyu3C6sXEya1bt2aRx999A/7n3rqKSZNmsSgQYMotxtxffHFF7N9+/ZKG+UnT57MTTfdxIUXXkhk5NnVRP369eOqq65iwIABPPfcc7Rs2dI1J1QL1rz/IC3TZ/HXEV0Z3Km11XKcyt2DY7mi0U5yZj+kuxHXE3al/EQox/HuNNwlx+/cdxif2y5lxV7P9uRo2fT1IuIN7AIuBTIxXPmOVUptryL+w0AvpdRd5vYxpdQ5ucbT09efG1ZdmzVfTqX/tr+xpulN9H9gmtvzdwdJs15iQOqrbLl4BglDr7dajqaWrHr/IRKzvqDkyXSCQ5u4JI9xb35LZ/YxqZKPTlfj6PT1VpZQ+gHpSqndSqkSYDYwspr4Y4FZblGmsYwtv8ynz9YX2RTYj8QJ71ktx2X0uu4xsoki4NeXdCmlHjClaAQvNnvdZcYE4P6AxTyZ/1dOHCt0WR61xUqD0grIsNvONPf9ARFpC8QCS+12B4hIiogkiUiVXZREZIIZLyUnJ8cZujUuIjN9K22X3k+Gdxva3TfHIxwLuQr/gEZkJDxIfFkaW375xmo5mlqwL+84O3NLiekxxKX5BHcehp+Uk5682KX51AYrDUplXYuqqn8bA8xTStlP2RptFsFuBt4QkbjKEiqlPlBKJSqlEqOioio9eEPyWuko7r4mJ0rK+PqrzynHC//b5hIS5rovPU+h59X3c4gIfFa9brUUTS3Ys2IuE33mMDSuZp/xtaF9n2GUKB9O7Fzi0nxqg5UGJRNoY7fdGqhqhN0YKlR3KaWyzN/dwDKg1/mICAgIIC8vTxsVO5RS5OXlERBw/lNvn2t+f/5qC28UDGbHDT/TMraTW/K1Gj//AFK7Ps4XxxNJ3uPZvXc0VRO662tu9FlJbNNwl+YTGBRCmn8XInPXujSf2mBlnUIyEC8iscABDKNxc8VIItIRaAysttvXGDihlDolIpHAIOC1imkdoXXr1mRmZqKrw84mICCA1q3d07tq5ex/kLNFmHj5dQxMaO+WPD2FviMf4LGdS8lctpuZsXV30GZDpbysjLjj69gVPpTmbpio9GizfvTb/xFFR48QEtrY5fmdK5YZFKVUmYg8BPwIeAPTlVLbRGQKkKKUOt0VeCwwW51dhOgMvC8iNoxS1itV9Q6rCV9f37NGkWvcy9aV3zFg58sENB5Cn6ETrZbjdgL9vLlnYCsOLP2AtC2K+IT+NSfSeAy7t6winuNI3EVuyU/63s2gtB68drCMIaFuyfKcsLTVUym1CFhUYd/zFbYnV5JuFZBQcb+mbnEwI52Wix/ggHdLOt07w2NH7Luam3s3xfeXL9i2JBsS5tScQOMx5G7+kXggpu+VNcZ1Bt06duCw7CFlbz5DOlTeJmwleqS8xhKKTx7n6Mdj8VOlyJjPXdrd0tMJaxzJ5qir6VGwhNyDGTUn0HgM2blH2ObdmcjmbWqO7ASC/X2YELGJ6C1vuyW/c0UbFI3bUUrx/SdT6VC2i7RBU4nuoCdJbHX5o/hJGWmL3rJaisZBikvLefrINXzd8yO35ntR4B6uOjqH0pJTbs3XEbRB0bidWWszeGJPb2Z3eY9elznHEVFdp018DzYG9qfD/jmcKq4b7l4bOuv25lNSZmNQvHs7U/jGDiRQStizdXXNkd2MNigat7Jj/a98sOBnhnRoxk03jrFajkfhNeB+9tqasjR5i9VSNA7gtfwVFvlNom+0e1vHo3teDED+juVuzdcR6u9QZI3HkXswg4gF43jXvwktRq+qc14XXU3ChSO5dF0TAjaWMHywarCdFOoKoYfX4OvrS0ijQLfmG9Uyhixphl9WSs2R3YwuoWjcQmnJKQ5PH0uIOobfde8QHuRvtSSPQ7y8uGNgDAcOZLI9dafVcjTVcKr4BHGndpIXWeN8iS4hM6Q7J0963szD2qBo3MK6aQ/RpWQL2xL/TlzCAKvleCwjEyL52X8ix3/8m9VSNNWwZ9MKAqQUv3aDLcl/c+Kr3HLiSQ4XFVuSf1Vog6JxOWsWfcyAw3NJajqaxGvutVqORxMSFMSuJhfRLX8JRYX5VsvRVMGRnUb7RUyvSyzKgK9KAAAgAElEQVTJv0e0MUp+c4ZnzTysDYrGpWzLKuS+1aHMCR5Hn/Ge2Xfe0wi/cAKN5BTbf3Rvd1SN46w63pK5ftfRpGmlE6S7nK4tgvnE7xX8k960JP+q0AZF4zIK83N47NOV+AeGcsm9U/H10+0mjhDf80J+825HZOoX2leKB1JuU8w81J6NnZ+0TEMjfz9a+RYRdjjJMg2VoQ2KxiWUl5Wx74Ox/Ov4s7x7cw+iQrQxcRTx8iK341jiyneza/sGq+VoKrBrz17CSw7QP8bayRnzwroRU7zToz46tEHRuIS1M56ke3EyJxJupVeMnkX3XOl8+XguL3+Tmbt8rZaiqUDh2i/41f9x+kdaPFK9VR9COU7m7m3W6rBDGxSN09nw4ydccGAmaxtfTf+brKsWqMuEhjUhoXsv/rspm+LS8poTaNyGX2YSWdKU5m0q9ennNiI7DgTg0PaVluqwRxsUjVPZt3MdHVb9iV0+Hehx74dWy6nTjOraiKm219i6+FOrpWhMlM1G9PHNZIV0t1oK0R1787PqTWqh57zGPUeJps5TVFzK09+ls9GrE2F3zMY/oJHVkuo0iR3b0cN7L75bPrdaisbkYEYakRRQ3qqv1VLw8fXlPy1e5KuirlZLOYM2KBqnYCsv58k5G0g+EoLPbV/TrLW11QH1AS9vb/a0vIquJ1LIPbjfajka4MDWXwCI6GTNgMaKdG0Zxt7sHMrLPaNa1FKDIiLDRSRVRNJF5OlKwu8QkRwR2Wgu4+3CbheRNHO53b3KNRVZ88kz3JA+ieeuiKN/uwir5dQbWg65Ax+xkf7TTKulaICfT3XmMdujtO1sfQkF4DJZQ4rXHWSme8aEopYZFBHxBt4BrgC6AGNFpEslUecopXqayzQzbRPgBaA/0A94wfQzr7GATUvn0n/v+0Q0bsy4QfFWy6lXtO3Um10+HYja/Y3VUjTAimwhq9UVHjOmqllMF7xEkZO21mopgLUllH5AulJqt1KqBJgNjHQw7eXAYqVUvlLqCLAYGO4inZpqyEzfSuwvj7LHJ5au985AvHQtqrM52PkuFhT3YscBPRWLlRSfPE7fg7MZ2uyk1VLOEN2xFyXKh7LMjVZLAaw1KK0Ae3+nmea+itwgIptFZJ6InPaz6WhaRGSCiKSISEpOTo4zdGtMjhcVUPbFWGx4EXjbbAKDQqyWVC/pNvxu3lE38s2mQ1ZLadDs3bqav3h/wsBGWVZLOYOvnz/7fdoSdGS71VIAaw1KZc4eVIXt74AYpVR3YAnw8TmkNXYq9YFSKlEplRgVFXXeYjVno5TijS8XE1h+lIxL/k3LmI5WS6q3NAnyY1h8OHnrv8XmIY2vDZGCXcZ4jzbdh1is5GzyQzvR+lS6R4yYt9KgZAJt7LZbA2eZfqVUnlLq9HDUD4E+jqbVuJa3fkrnw12N+O/QhSQMuc5qOfWe8VFb+WfZK6Su/Z/VUhosvtnryJKmRDaPtlrKWRTEjeSN0uvJLrDeP4qVBiUZiBeRWBHxA8YAC+wjiEgLu80RwA5z/UfgMhFpbDbGX2bu07iB9T/MRJb9nRt6teCuiz2nD3x9pvPQUZxUfhxdN9dqKQ2W1se2kh3czWoZfyAi4VI+Lr+cbQdPWC3FOhfASqkyEXkIwxB4A9OVUttEZAqQopRaADwiIiOAMiAfuMNMmy8if8MwSgBTlFK6xdINpG9aSafVTxEUEEvbEZ20m1o3ERQSzvqQAbTP/Ymy0hJ8fP2sltSgOJS9n0iVz56W1nhorI5OzUNpI4c5kLYJulxmqRZLfcorpRYBiyrse95ufRIwqYq004HpLhWoOYucg/sJ/mYcRRJMxPgvCQjUI+HdStfriVjzC1uTvqfbhY52iNQ4g3W5vkw89RGz+/WzWsofCPL34bPAqeSnRmNU1liH7uOpcYhTxSfImzaKMHWUous+9bh65IZA5yE3ckL5k79pUc2RNU5l/b4jlPk0olNbaxxq1UROUEdanNhltQxtUDQ1o5Ri2qy5xJamsWPAa7TvMchqSQ2SwKAQpsZN55G86yktt75HT0Oi5/ZXeKxJEn4+nvnKLIvqSnNyKcg9aKkOz7w6Go/iveW7+UdqFJ/2+4beV9xptZwGzQWJfSk4WcaK9FyrpTQYSktOMez4InoHWvuyro6gtr0AyExNriGma9EGRVMtKYums/V/M7m6ewvuuuJCq+U0eIZ0iOT5gNnIj89aLaXBsD91AwFSik/r3lZLqZIW8caIimP7N1uqQxsUTZVsT/qR7mv+xINBS5l6YwJeXrpHl9X4+3iTEFZC7/z/UnKq2Go5DYK8XasBaNZ5oMVKqiaieRue9PoTPzHAUh3aoGgqZf+ujbT84S4OeTWl5b1fEeCnXdF6Cr7dRhLKCVKTdOO8O1BZGzhKI1rFeu6YK/HyIqv5MNbmB1iqQxsUzR/IO5SJ96xR2PDC67avCY9sbrUkjR2dBo3ghPLnxOb5VktpEOQcK2NnYC+8vD37dXlBeCG9D31l6fQ8NY5DEZFE4EKgJXAS2Aos0QMJ6ycnS8qZ9+m/ud2Wz75r5tKxXWerJWkqENAomPUh/YnLW4atvBwvb2+rJdVbikvLeazoVu4Z0g7PG4FyNgO8t9PP6yMO7LubVu2sKU1VaXJN51brMQYWBgKpwGFgMLBYRD4WET0YoR5RVm7j4VkbeCXvQtZe+T0dEy+xWpKmCkq7juKHsj5s2p1ptZR6zc6DRZTZFN1bhVktpUbC2vYA4HD6Bss0VFdCCQIGKaUqnfxfRHoC8YD2TVoPsJWXs+Kd+8jK7snka65iSP8YqyVpqqHzxWO4dUUEd+w6Ti/t08xlnFw9je/9PiYsarHVUmqkdcfe2L4Tig9Y572xSoOilHoHQEQGKaVW2odVtk9Td1E2G2s/eICL8ufi0yGKwQNjrJakqYHQAF8GxzVhz5aVqCs6acdmLsInK4VmXgU0btrUaik1EhQSTqZXM/zzdtQc2UU4che+7eA+TR0l6ZNnGHBoNklRNzHo9r9bLUfjIPeGJzOteCK7t1k7mK0+E3l0OxmBdcdg5wTGEXki3bL8q2tDuUBEngSiROQJu2UyxuzAmnrAmrmvccHed0kJvZR+971fZx4cDcQPug6bEg6v/dJqKfWS40WFRJfv50Rkd6ulOExy56e5+uRfKS61pqdXdW8PPyAYo1osxG45CtzoemkaV/PN+gyObVnExsAB9Hjoc91bqI4R0aw1O/260vyA59fv10X2bUvCWxSBMX2tluIwrdrGc1Q1Iv3wMUvyr64NZTmwXERmKqX2uVGTxg3M37CfJ7/cwuDYKXxwa298/fytlqQ5D47GDqfLrqlkpm+ldXvPc/5Ul0nNK2NX+UAGdfXcEfIV6RThzUSfOeRtPQmt3P/dX12V13cisgB4W0QWVFyckbmIDBeRVBFJF5GnKwl/QkS2i8hmEflJRNrahZWLyEZzcYqehkLKwg+J+WYkF7f15f07BhHQKNhqSZrzpO2gUQBkrNbVXs7m56MtebXRRKLqkKuGtk0bM957Ef67l1iSf3Xdhqe6MmMR8QbeAS7F8BGfLCILlFLb7aJtABKVUidE5H7gNWC0GXZSKdXTlRrrI+sWzaDn2qfYFdCVt2/tT6Cfruaqy7Ro25Gngl/iwLGuXGC1mHpG5v7ddG9Vd4wJgI+vL3t9omlUaI1vlJqqvFxJPyBdKbUbQERmAyOBMwZFKfWzXfwk4FYXa6rXrFs0gx5rniDNrzNtH/4vjYI9f7CWpmaadR/GvJ/TOXK8hMZB2jWwMyjMz+Hrk3exgseButOGAnAkOJ7YwiRL8q6xS4+I7BGR3RUXJ+TdCsiw284091XF3cD3dtsBIpIiIkkicm1ViURkghkvJScnp3aK6zCrFkyj55rHSfPrTJuHFxIUEm61JI2TGBYfyqPeX7L916+sllJv2L91FQBhMT0sVnLulEd1JpICjuRkuz1vR3zKJ9qtBwA3AU2ckHdlc6GrSiOK3GrqGGq3O1oplSUi7YClIrJFKfXbHw6o1AfABwCJiYmVHr++M+3X3Xy0yo+XIy6n/wPTCAwKsVqSxokkRDellc/P7NtWAMPHWi2nXnBsz1oAorvVPe+kQW26czQtkIw9qTSOauHWvGssoSil8uyWA0qpNwBnTPKUCbSx224NZFWMJCLDgGeBEUqpU3a6sszf3cAyoJcTNNUrlM3Gf2e9y0sLt9EroSsDH5+ljUk9xMvbm93hA4kvWkNpyamaE2hqxP/wJjKlBWFNPH+EfEWiegyn+6lpbCqPcXvejlR59bZbEkXkPozxKLUlGYgXkVgR8QPGAGf11hKRXsD7GMbksN3+xiLib65HAoOwa3vRQHlZGWvencDVqU/zUtx23hrTy2P9YWtqj2+XqwwfKWv/Z7WUekHL4zs4GFw3Z9puFhZIiL8vuw65fyyKI1Ve/7RbLwP2AqNqm7FSqkxEHgJ+xBh5P10ptU1EpgApSqkFwD8wBld+KSIA+5VSI4DOwPsiYsMwiq9U6B3WoDl5vIid74xmwImVJDUbw6i7Jnq8LwdN7eg48GpKVj7Ksc3/hcHXWC2nTpNztJiXS0Zxafu6WekhIjwd/F+a7soFZrk17xoNilLqYldlrpRaBCyqsO95u/VhVaRbBSS4SlddJvdgBvnTrqdHaRpJnZ5iwFjte7wh0Cg4nOSgC8jKP4pSCvMDTHMebMkqZL5tMDcnWOtOtzbE++XR/sgKt+frSJVXuIg8IiL/EpG3Ti/uEKc5N9IPH+PZGQtpWprFpkH/1sakgbFz8Ns8cexWfsuxZtqN+sLBnUl09dpL1zrgA6UqbJGdaMJR8g8fcGu+jtSDLAJigC3AOrtF40GsWr2C695Zybry9mSMW0Ovy/SQnYbGJV0MV83LtmbUEFNTHV1T3+E/Ae8S7O9Ii4Bn0qiV4bExO32TW/N15IoFKKWecLkSzXlhKy9nzcynGJgxjdGhz3Hn3Q/SKjzQalkaC2gVHsiM0PdptjofLnF/dUd9QNlstDqZyp6w/rStObrH0qx9T1gGxzK3AFe6LV9HSiifisg9ItJCRJqcXlyuTFMjhUdy2Tz1Si7ImMba8CuZ+MD92pg0cAKaxtGxZCsFuQetllInycneRyQFlDevewMa7WnaMpYtKo7sIvdOY++IQSnB6G21mt+ru1JcKUpTMzuTf+L4WwPpeiKZNV2epe8jnxMQGGS1LI3FRPQeibco0ld9Y7WUOknW9tUAhMX1s1hJ7RAvL55v9jZzyp0xZNBxHDEoTwDtlVIxSqlYc2nnamGayim3Kd7+KY335i8FFL9dPZf+o57SjrE0ALTvcSF5hEGa9pFyPpzcv45yJbTt2t9qKbUmvmkwaYeL3JqnI2+hbcAJVwvR1EzW3lT+9fa/+OfiXdi63kjwE+vo1LfSntWaBoqXtze7w/rTvmgN5WVlVsupc3zKVTwZ/AqNgkKtllJrLiOJhaXjOZLzhwlIXIYjjfLlwEYR+Rmwn/rkEZep0pyFrbyctXNfpfvON7gLf9pft4Rr+3XQYw00lXKq+zheXtqWMRm59IptbrWcOoNSirXZNoZ0qLvjT+yJioykmRSwPX0TjaNauiVPRwzKt+aisYD9uzZy/Mv7GVC6nc2BfWl687tcFx1vtSyNB9PtgssZt8SLZmmF2qCcA4ezM7j55Cyim9xptRSn0CyuByyDoowtcMEVbsmzWoMiIt2VUh+LSIJSaotbFGkAOHaqjI9+TOLedddySvxI7vkSiSPu120lmhoJb+THJS3LKNnyDVz2B0eomirI2voLT/jOY2dY/ZixuVmrdhxTgXB4p9vyrKmEcpeIvIPhi+QxN+hp8NjKy1m+5DueSgkmp+gUzdpNZNiIm+lbh9yQaqzn7uBVXLDvPfIO3UpEs9ZWy6kTnNq/3miQ71L3G+TB6Ol1wDeaoKN/8OrhMqrzKf+CGZ4EeInI81XF1dQeZbOx4X+fseelRC5efTuDg7OZ/+Agxkx4mkhtTDTnSGSvqwHYnbSghpia0wTmbma/d3S9cvGwK2IYK8s6ui2/Kg2KUuqvwBJgDrBEKTXFbaoaEMpmY9PS2aS/2Jdeqx7E33aSlD6v8c+HbqZHG+1VUXN+xCUMJI8wJH2J1VLqBMpmo3XxLnJCulgtxakc7Dqel0+M5MjxErfkV1OFfH+l1APUNafKdYDi0nJmrd3PdW/8QNzyRwmyFbG2x99p/sxmEq+5Fy9vb6slauowRvfhAbr7sIMcPnSAQFWMalG3R8hXpH2zYPwo5bcs97g/r9agKKWeFZGmSqnn7PeLiPvKUPWMPdvWkPTufWx6cSiTvt5MqU8Iay6cQdSkLfS77mF8fP2slqipJ0iHSwnnGKmbk6yW4vFszPel26mP8O17h9VSnEqnoONs97+T8g1fuCU/R7oN/yoizyml5gKIyJMYjfT1q2zoQjL3pZG54guidn9LXPluWitvtgQP5MuxCSR2iNbjSTQuof2gG+m3ohFjcqP0w1oDWw8UIl7edG5T91z+Vkfzlm05gR/kuKenlyMG5SLgAxG5CWgG7ACcMtGNiAwH3sTw2DhNKfVKhXB/4BOgD5AHjFZK7TXDJmEYtnLgEaXUj87Q5AyKCvPZs+Fnfi2I5Nvd0Dn3f7zl92/SvNuT1PHPdBx2J72jWlgtU1PPCQ9vTJvoWH5OzeGJy3SlQnV03fIqz4YGEejnvpl53cHvPb3S3ZKfIx4bs0XkB2ASYAMmKaVq7cFHRLyBd4BLgUwgWUQWVHDlezdwRCnVXkTGAK8Co0WkC4YP+q5AS2CJiHRQSrl1as2ychu5x0rYm3UQ2fgZkrOTJoXbiS3bTXdRzC8bR9OYW0jsfTMHom8lvl1n9JBEjTsZ1SKHyHWvk3vwUyKbt7FajkeibDb6FS0mPXyw1VJcQkFwe+IKVrklrxoNiogsBrKBbkBrYLqI/KKUmljLvPsB6Uqp3WY+s4GRgL1BGQlMNtfnAf8Wo35oJDBbKXUK2CMi6ebxVtdSU6W883M6/VMex6/8OD7lp/CznSS8PI+vywbxUtktBFLMjoB/UEAwmf7tWdtyPMHxF/JEr6EEheieWhrr6B3ThPiNG0hOWkDktQ9aLccjOXRgN805Wu8a5E9ji+xARMEiCvMOERbRzKV5OVLl9Y5S6vTUKwUiMhCjtFJbWgH2ruUygYojis7EUUqViUghEGHuT6qQtlVlmYjIBGACQHT0+Y3nyCo4SWBxDj5io8zLn+N+keQHdKZF00H8PbYbsZFB5IRsIzKqJeF6JLvGg4hLGEjet2F4/bYE0AalMrK2r6Y5EBZXPzuzlvccxw2H+/NqeRCudmpcpUEREVEGZ83jpZQqA/5mH+c8866sJbrisaqK40haY6dSHwAfACQmJp6X1hevS4DrdE8ZTd3D6D58AfGFKygvK8Pbp+66tXUVp/avo0x5EVNPRshXZFC3OAZ1i3NLXtV9Tv8sIg+LyFmf9SLiJyKXiMjHwO21yDsTsK/UbQ1UnGf5TBwR8QHCgHwH02o0GsDL7D6ctmGZ1VI8kgPHFOt9exPQKNhqKXWe6gzKcIweVLNEJEtEtovIHiANGAu8rpSaWYu8k4F4EYkVET+MRvaK80Qs4HejdSOw1CwRLQDGiIi/iMQC8cDaWmjRaOot7S8YwXpbezbvPWS1FI9DKcVLx65mXqd/WS2lXlBl+VcpVQz8B/iPiPgCkcBJpVSBMzI220QeAn7E6DY8XSm1TUSmAClKqQXARxg+7dMxSiZjzLTbRGQuRgN+GfCgu3t4aTR1hbAmTXmxxduUHLQxymoxHkZWwUnyj5eQ0MrVrQsNA0crVG0Y7RahIhIKoJTaX9vMlVKLgEUV9j1vt14M3FRF2heBF2urQaNpCFzUIYp3Fm8h90gBkY11z8PTHF77JSv8p1AYMg+IsVpOnafGLkki8jBwCFgMLDSX/7pYl0ajcSKXNT/OJv8J7PnFPVNw1BWK96+nGUeIi+tgtZR6gSMllEeBjkqpPFeL0Wg0riG+UwJHpJE5+/ADVsvxGILytrLfJ5q4wCCrpdQLHBk0kQEUulqIRqNxHUb34f569mE7lM1Gm+JU8uvZlPVW4ohB2Q0sE5FJIvLE6cXVwjQajXOReKP7cPrGX6yW4hEczEijMUXYWva0Wkq9wRGDsh+j/cQPCLFbNBpNHSJuwAhsSsjftKjmyA2AndlFfFJ2KWEdh1otpd7gyOSQf3WHEI1G41oaRzbn3yGPsPNYZy6wWowHkFwQzAe2u9jaJdFqKfWG6nzKv2H+ficiCyou7pOo0WicRVnPW1l4MJR8N7mE9WSy9u2iU9NGBPhq76jOoroSyqfm71R3CNFoNK5naPvGbF+azJa1Xgy9+Aqr5ViGstmYnH0/qU0uwXD5pHEG1Y2UX2f+Lj+9T0QaA22UUpvdoE2j0TiZ7q0b86rfNNI3boMGbFCy9+2iJcegeXerpdQrHBnYuExEQkWkCbAJmCEieuIbjaYO4u3jw2+h/YkrTMJW3nBnK8resRKAiI66NcmZONLLK0wpdRS4HpihlOoDDHOtLI1G4zLaD6MJR/lti3u8+HkipfuSOaV8adu5fvpAsQpHDIqPiLQARqGnXNFo6jyx/a8BIHdDw32cQ/O3sNc3Dl8/f6ul1CscMShTMGYETldKJYtIO4wp7DUaTR0kollr0nzi8T64yWopllBWbmNq8UiS295jtZR6R40GRSn1pVKqu1LqAXN7t1LqBtdL02g0rmJxr3cYc/QhCk40vO7D6TnHWFraleBuDbdTgqvQDtA1mgZI/27x2JSwIj3XailuZ9+2NQz02kqPVnrCD2ejDYpG0wDp2aYxUwK+IGj5FKuluJ3wbZ/wnt+bxERql7/OxhKDIiJNRGSxiKSZv40ridNTRFaLyDYR2Swio+3CZorIHhHZaC56djeN5hzw9hK6hhwnIe97bOU2q+W4lYjCbewN6IR46e9pZ+PQFRWRS+x/ncDTwE9KqXjgJ3O7IieAcUqprhj+7d8QEXtXc39SSvU0l41O0qXRNBjK2w0jkgL2bEuyWorbKD5xjJiyPRyP0AMaXYGjJnpqhd/aMhL42Fz/GLi2YgSl1C6lVJq5ngUcBqKclL9G0+Bpd8EIAA43oO7De7euxkdsBMb0s1pKveRcy3zipHybKaWyAczfptVmKtIPY/r83+x2v2hWhb0uIlV2JheRCSKSIiIpOTk5ztCu0dQLIptHk+4dR1hmw/GPUpBulMZadx1ksZL6icsqEUVkiYhsrWQZeY7HaYExUeWdSqnTlb2TgE5AX6AJ8Oeq0iulPlBKJSqlEqOidAFHo7FnX+trSC5uydGTDaP78Fy5nNt9pxLZsq3VUuoljviUPy+UUlVOzyIih0SkhVIq2zQYh6uIFwosBP6ilDpT0Xu6dAOcEpEZwEQnStdoGgwhFz3KC6mraZqexxUJLayW43LWZx6nU3Rvq2XUW6zq5rAAuN1cvx2YXzGCiPgB3wCfKKW+rBDWwvwVjPaXrS5Vq9HUU3pHhxMa4MXGLfX/ESrMO8Rthe9zUZN8q6XUWxw1KMfM3yIn5fsKcKmIpAGXmtuISKKITDPjjAKGAHdU0j34cxHZAmwBIoG/O0mXRtOg8PH24r2QGdy9awLKVr+7D+/bsoK7fb6nW9hJq6XUWxyq8lJKDbH/rS1KqTzg/yrZnwKMN9c/Az6rIr2zui9rNA0en5iBNN3yI3t2JBPbtb/VclzGiT1rsSkhuptukHcVemSPRtPAielvdB8+tK5+dx8OyNlEhncrQsMjrJZSb9EGRaNp4DRt3Y49XjGEZC6vOXIdRdlstD2xjUOhekCjK9EGRaPRcLDZYOJPbeXY0SNWS3EJGRn7KFMCrbVDLVfiiAvghyqba0uj0dQfAvrewYTSJ1m1z1n9bjyLtbm+9D31H8IH3Wm1lHqNIyWU5kCyiMwVkeFmV12NRlOP6Nq9Dym+iSxLL7RaiktYv/8IIQG+tG8WXnNkzXnjiIOtvwDxwEfAHUCaiLwkInEu1qbRaNyEn48X17c5RqttH9TL7sMjtj3OM+FL8PLS38OuxKE2FKWUAg6aSxnQGJgnIq+5UJtGo3EjV4ft4cGyT9ifusFqKU6lqCCPfqUpxIRaraT+40gbyiMisg54DVgJJCil7gf6ANoVsEZTT4jub0yzd3DddxYrcS57N/+ClyiC2w+0Wkq9x5GBjZHA9UqpffY7lVI2EbnaNbI0Go27aR4dz16vNjTKWGa1FKdyLH0VNiXE9BhqtZR6jyNtKM9XNCZ2YTucL0mj0VjFwajBdCzewvGiAqulOI2gw+vZ692WkLAmVkup9+hxKBqN5gzB3a6gDG+2b1prtRSnYLMp1hc3J73ppVZLaRBog6LRaM4Q3384A9U0FuS2tFqKU0jPOcbk4ps52vcxq6U0CLRB0Wg0Z/D38ycxrgXLdh3G6NxZt9mUnoFgo09bPTbbHWiDotFozuKGZgd579hjZKZvtlpKrYlJnsKygD8RE9HIaikNAm1QNBrNWfToGE9Xr30cSF5gtZRa07xwE0caxaAn+HAPlhgUEWkiIotFJM38rbQ8KiLlds61FtjtjxWRNWb6OaZ3R41G4wRaxnZir1cbgvctsVpKrcjN2kcblUVxy/rr48XTsKqE8jTwk1IqHvjJ3K6Mk0qpnuYywm7/q8DrZvojwN2ulavRNCyym19Mp+LNFObnWC3lvNm30TCITbpof3zuwiqDMhL42Fz/GMMvvEOYk1NeAsw7n/QajaZmGve+Fh+xkbbqG6ulnDdlu1dwQvnTLuECq6U0GKwyKM2UUtkA5m/TKuIFiEiKiCSJyGmjEQEUKKXKzO1MoFVVGYnIBPMYKTk5dfdrS6NxJx16XcRCGcrKgw55CfdI5p7sx6wm9+Pjq2vE3YXL7hYRWYIx9X1Fnj2Hw0QrpbJEpDBGCRkAAA6KSURBVB2wVES2AEcriVdl/0al1AfABwCJiYl1vx+kRuMGvLy9+bXb31m4OZv7y2z4+dSt/jtHjpfwVV40sZcNs1pKg8Jld4lSaphSqlsly3zgkIi0ADB/D1dxjCzzdzewDOgF5ALhInLaGLYGslx1HhpNQ2VY52aEnDrIhi11r/vwti0p9JWd9GsbZrWUBoVVnx0LgNvN9duB+RUjiEhjEfE31yOBQcB2cyr9n4Ebq0uv0Whqx+DYEJb4/wlWvmm1lHPGZ8MnfOb3Mj1aBVstpUFhlUF5BbhURNKAS81tRCRRRKaZcToDKSKyCcOAvKKU2m6G/Rl4QkTSMdpUPnKreo2mARAQ2IjUoD60zf2lzjndishLId2/E/4BekCjO7GkxU0plQf8XyX7U4Dx5voqIKGK9LuBfq7UqNFooKz9cJpvXkX6ltW07zHIajkOUVSYT7vSdJLb3GW1lAZH3Wpp02g0bqXdoBuwKSEn5VurpTjM7nVL8BZFUAft/8TdaIOi0WiqJKJZa3b5dSYq6yerpTjMidSlFCtf4vvqHl7uRhsUjUZTLZt7Ps+Y4xPJPHLCaikO8WLxTfy12dsEBAZZLaXBoQ2KRqOpln79h5BLGD9sPWi1lBo5XFTM1kPFRHfRTaxWoA2KRqOplpjIIO6O2EbzVZOtllIj6Su/ZZLP5wxpG2i1lAaJNigajaZGhkXlc/XJ+RzMSLdaSrX47vyW0T7L6dSmqtmcNK5EGxSNRlMjrQeNBWDvr7MtVlI1ymYjumAtu4N74+1Td+cgq8tog6LRaGqkTXwPdnvFELZnodVSqiQjfTPNyKO0re4ubBXaoGg0Goc41GY4nUu3c/jAHqulVEr2hh8AaN3nKouVNFy0QdFoNA7RcuBYNtriSNq8w2oplbLvUD6p0o5W7TpbLaXBog2KRqNxiLYde/JU49f5bF+lHrstpbi0nBdyLuGz7p9YLaVBow2KRqNxmCsTWrBjXyaHPcxZ3f+3d+9BWtX3Hcffn72BCwhyKSqXQIVANiBYVwPSNnXRxpgWTBWjaS2DNtY0Gk0ziWiaTGYa46WZGtOoE+sldgpGa6w68YI36mWiyAZEwA0BwbiLUAi3ILrA7n77x/OQrnLZZfc8z2/J83nN7Ow+Z89zzuc8sPvZ3znnOWfR6g28v7eVupqhqaOUNBeKmXXazDGVvFr1D6x96rbUUT6g93Pf4rFe1zJ1dM8bPZUSF4qZddroUaNprBzF4LU95xZE0dbGiC0vsrfPcfSuqkwdp6S5UMzssGw94RzGtL7JWw31qaMA8PaqpRwfm9g92heDTM2FYmaHZWzdbFqijI0v3ps6CgAb6nOjpVFTz0mcxJIUiqSBkp6WtDr/eb8dn5JOl/Rau49mSefkv/djSevafW9y8bfCrDQNGjqcldW1jH7nMdpaW1PH4ejGZ1lbPoqhw8ekjlLyUo1Q5gLPRsRY4Nn84w+IiIURMTkiJgN1wHvAU+1m+dq+70fEa0VJbWYAbJt6DX+9+2peeWtb0hybd+7mtl11vDHm75PmsJxUhTIT2DdevhfoaKx6HvBERBwZN2Qw+z03Zeon2Vg1ioeXrk+aY8HKjfysdQpjT78oaQ7LSVUoQyNiA0D+c0eXBr0AuO9D066T9LqkmyX1OtgTJV0qqV5S/eYedu682ZGqd2U5F495l2nL/4ldO7cny7H91fmcNuhdPjq0b7IM9v8KViiSnpG04gAfMw9zOccBE4EF7SZfA4wHTgEGAlcf7PkRcUdE1EZE7ZAhQ7qwJWZ2IGePO5qZeoGVT92TZP1bN63nsi03ccWAl5GUJIN9UMEKJSLOiIgJB/h4BPjffFHsK4xNh1jU+cB/R8TedsveEDm7gXsA357NrMjG1U5nXdlHOKZhXpL1r3nxASrUxpBTZyVZv+0v1S6vR4HZ+a9nA4d6l9SFfGh3V7syErnjLysKkNHMDkFlZWz66IWMbVnNmmUvFX391b96mCYdywkTpxZ93XZgqQrlBuBMSauBM/OPkVQr6c59M0kaBYwAnv/Q8+dJWg4sBwYD3ylCZjP7kPGf+gLvRxVbnv9RUde78e3V1DQvo3HEDFTmt9P1FEluaxYRW4DpB5heD/xdu8dvAcMOMF9dIfOZWef0P2Ywzw05j2Wb2vj47hb69irOr5RlLz9NHWWM/LOLi7I+6xxXu5l1yzEzvsste2bw0JKmoqwvIripqYYvDL3P9z7pYVwoZtYtk0cM4OQR/Xjjfx6gtaWl4Ot7vXEbb27exadqXSY9jQvFzLpFEteOaeSG3dfx+rOFP+Or+eEr+XGvf+HsCccWfF12eFwoZtZtk6d/jvUaylGLbyciCraeHVs3c+KWJ6keOIz+1VUFW491jQvFzLqtvKKCpvEXM76lgV+++nTB1tPwxG0cpT0MOv1LBVuHdZ0LxcwyMfEzX2Qb/di78MaCLL+ttZXha+bTUFnDCSeeVpB1WPe4UMwsE9V9+7PqhIvp914jS1aty3z5yxY+wPDYyK5JczJftmXDhWJmmZl07lwuqPwBNz2/MdPlRgQ3Ngzm+5WXMOnPZ3f8BEvChWJmmTmquprL6saxbO07LK5flNlyX1z9G15Zv4c/OOMqKqsOenFxS8yFYmaZ+vwnRjKv+maGPj6HPbubu728aGuj5aHLmNX3dc49eb8LZ1gP4kIxs0z1qiin/I+/zMi29Sy5/7puL2/JE3dT1/wMsz7Wm14V5RkktEJxoZhZ5ibVnc/S6qmc+OaP2PDrVV1eTvN77zJs8Q28WT6ak2dekWFCKwQXipkVxNDzb6GNMrbNu6TLl2R57T++zrFs5v26f6a8Ism1bO0wuFDMrCCOHzWOhpO+yY73W7nrmWWH/fzlS37OqRvms2jgDCZM+8sCJLSsuVDMrGBqZ3yR+2t+yPXPb+TJFRs6/bxNO5u5bMH7XNfrKibM+bcCJrQsuVDMrGBUVsYN503mT4aVUf1fn2P5C4e6OWvOzh1buf6u+9n63l4+O/sr9Ok3oAhJLQtJCkXSLEkrJbVJqj3EfGdJWiVpjaS57aaPlrRI0mpJ90vyVeLMeqjeleXcPGsiI8p3MO7ZOSx+5LaDzrtp/To2/uBMvrntG/zw3DFMGNa/iEmtu1KNUFYAfwW8cLAZJJUDtwKfBmqACyXV5L99I3BzRIwFtgGXFDaumXXHoKHDGXj5M6zpVcMpS69h6U1n86vXXvrdlYl37tjKK/ddR69/n8awlkYaP/k9pk8ekzi1Ha5UtwBugNx9FA7hVGBNRKzNz/sTYKakBqAO+Hx+vnuBbwO3FyqvmXVf/4FD6PO153h5/reZuPZuXv7pt5j1yFyGVO3lieaLmKJW3qiaSL9ZtzJp7KTUca0LevJ5eMOAxnaPm4BPAIOA7RHR0m76Qd8+K+lS4FKAkSNHFiapmXVKRWUVU2d/l99u/yqtS5Zx9vYh7N7bykvNV3Hcx6dRUzs9dUTrhoIViqRngAPdUu0bEdHxkTk40PAlDjH9gCLiDuAOgNra2sLd+cfMOu3oAYM4q66Os343ZXLCNJaVghVKRJzRzUU0ASPaPR4OvAP8BhggqSI/Stk33czMEurJpw0vBsbmz+iqAi4AHo3cUbyFwHn5+WYDnRnxmJlZAaU6bfizkpqAqcBjkhbkpx8v6XGA/OjjcmAB0AA8EBEr84u4GvhHSWvIHVO5q9jbYGZmH6R9p+2Vgtra2qivr08dw8zsiCLpFxFx0PcM7tOTd3mZmdkRxIViZmaZcKGYmVkmXChmZpaJkjooL2kz8OsuPn0wuffAlLJSfw1KffvBr0Gpbv9HImJIRzOVVKF0h6T6zpzl8Pus1F+DUt9+8GtQ6tvfEe/yMjOzTLhQzMwsEy6UzrsjdYAeoNRfg1LffvBrUOrbf0g+hmJmZpnwCMXMzDLhQjEzs0y4UDpB0lmSVklaI2lu6jzFJGmEpIWSGiStlHRl6kypSCqXtFTSz1JnKTZJAyQ9KOmX+f8LU1NnKjZJX8n/DKyQdJ+k3qkz9TQulA5IKgduBT4N1AAXSqpJm6qoWoCvRsTHgCnAl0ps+9u7ktytFErRLcCTETEemESJvQ6ShgFfBmojYgJQTu4eTdaOC6VjpwJrImJtROwBfgLMTJypaCJiQ0QsyX+9k9wvkmFpUxWfpOHAZ4A7U2cpNklHA39K/r5DEbEnIranTZVEBXCUpAqgGt8pdj8ulI4NAxrbPW6iBH+hAkgaBZwELEqbJInvA18H2lIHSeAPgc3APfldfndK6pM6VDFFxHrge8DbwAZgR0Q8lTZVz+NC6ZgOMK3kzrWW1Bf4KXBVRPw2dZ5ikvQXwKaI+EXqLIlUAH8E3B4RJwG7gFI7lngMuT0To4HjgT6S/iZtqp7HhdKxJmBEu8fDKbGhrqRKcmUyLyIeSp0ngWnADElvkdvlWSfpP9NGKqomoCki9o1MHyRXMKXkDGBdRGyOiL3AQ8BpiTP1OC6Uji0GxkoaLamK3IG4RxNnKhpJIrfvvCEi/jV1nhQi4pqIGB4Ro8j9+z8XESXz12lEbAQaJY3LT5oOvJEwUgpvA1MkVed/JqZTYicmdEZF6gA9XUS0SLocWEDuzI67I2Jl4ljFNA24CFgu6bX8tGsj4vGEmaz4rgDm5f+oWgvMSZynqCJikaQHgSXkznxcii/Dsh9fesXMzDLhXV5mZpYJF4qZmWXChWJmZplwoZiZWSZcKGZmlgkXipmZZcKFYmZmmXChmCUk6RRJr0vqLalP/n4bE1LnMusKv7HRLDFJ3wF6A0eRu2bW9YkjmXWJC8UssfzlTBYDzcBpEdGaOJJZl3iXl1l6A4G+QD9yIxWzI5JHKGaJSXqU3GXxRwPHRcTliSOZdYmvNmyWkKS/BVoiYr6kcuDnkuoi4rnU2cwOl0coZmaWCR9DMTOzTLhQzMwsEy4UMzPLhAvFzMwy4UIxM7NMuFDMzCwTLhQzM8vE/wHSmz4xQyFzuAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.title('Two results, once calculated in NumPy, once in Matlab.')\n",
    "plt.plot(x, np.sin(x*mult), label='NumPy')\n",
    "plt.plot(x, y, '--', label='Matlab')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y = sin(x*mult)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's it, very lightweight. We only need the packages `subprocess` and `scipy.io`, and defined two small functions, `matlab` and `save4mat`. With this we can happily continue to code in our Python notebooks, but call that one bloody function or script that we have to use in Matlab."
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}