Interactive plot of a SymPy expression (using Traits)

sympy_interactive_plot.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Interactive plot based on a SymPy expression"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Rationale : have a tool to interactively study an expression. First target : Bode diagram of a dynamical system (control theory)\n",
      "\n",
      "Pierre Haessig \u2014 October 7<sup>th</sup> 2013"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sympy\n",
      "from sympy import symbols\n",
      "sympy.init_printing()\n",
      "\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Laplace variable \"s\":\n",
      "s = symbols('s')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "A dummy example"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A first order system \n",
      "\n",
      "$$ H = \\frac{K}{1+ \\tau s} $$\n",
      "\n",
      "with the parameters $K$ (gain) and $\\tau$ (time constant) being sympy symbols that could be set interactively"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "K, tau  = symbols('K tau')\n",
      "H = K/(1+tau*s)\n",
      "print(type(H))\n",
      "H"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<class 'sympy.core.mul.Mul'>\n"
       ]
      },
      {
       "latex": [
        "$$\\frac{K}{s \\tau + 1}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAADkAAAAsBAMAAAAtLQ2eAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdqvNEDJUuyJEiWaZ\n3e/xv6KKAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABLElEQVQ4EWNgQABG/c8qDEz/TREiKCx+AwYG\ntgoUISTO/AUMDNuQ+KjM9RMYvBpQhZB49gysC5C4aMzPDHvQRJC4nL+4/iJx0ZhcHxecT0ATQ3CZ\nfzP0A/2EA/AXMLB/xSHHwDBfgIHzE07Z9QEMDPcdcEnfB0rIL8Auy6L8X5GB8f73BuzSI0P0Px7w\nYcgFgQk+F7d/wSLL0gARdD+ETZZNAKqFkRjZhdIBDPLA8IRoQtPLHMB+gGtpU0QHxEQ02X4GVocw\nBnaobQxosv5vJgBlgnDIctZ/B8qsBstyGRvbPjY2PgDigN3MkcAi38DAcBUkAAKoJjM/YOBnYOCA\nRzSqLHsAgyDQmAsQneh6OQQ3JjAwMEE8i2Qyj95fbZgOBA03GSGExOLcgMQBACcVXk+CkCZNAAAA\nAElFTkSuQmCC\n",
       "prompt_number": 4,
       "text": [
        "   K   \n",
        "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n",
        "s\u22c5\u03c4 + 1"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Extract all the symbols, using the `H.atoms` method"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "H.atoms?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "H.atoms()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\left\\{-1, 1, K, s, \\tau\\right\\}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAIMAAAAVBAMAAACXjEALAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZpkQzYnvq1QyRLvd\ndiJ+ofBJAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABzElEQVQ4EZ1UPUsDQRB98XLEuyQmWNmIERRL\nI2hlkYBYWJlGKyVip5WiCILB/APT2YmVpckvUH+AEO0sDF4vGMWIiAnnzu7dzuWjiNli572ZeW/3\nbpcFEJrEoMOqpUlqLA9qIHTWDolD70GLJU00ip82DhBxV3QlCKaImFuB1OunTxghkQWGb/xCe9wg\nGrRYnPEtGAHlHDDWrtRsvdMCcd8igJC/xIKjRe2gT4tVWGIfvUefFg2Msn7zsSRJ3HXdFiAtIlmu\nB7bPn2R/mU3dEi5FU5KsXRy/OcAEkTPlKvM9LcxWbraoykAdVpqwVUSFYjQHnNPlmL+nIU6e12YU\n/kA9S+00MtuXMgqTPQlOksBVyctRYCGjRAXRH7/Hvv31oCm0QDglpkiFoDdYyKj8APvbq48UjRdH\n4ZhcWv4L884rU2Aho7zoraZVU3gLCa+9LmPvQ81kqUgWClUFe8kpFi2hBqNBDc80qUMNXvCho+Y+\nYnRDNTIO3QLi1V9H5kdq40VvTwVy6LaQWVyrIGZGlNLMdojK0f0hMp30ygAjSmkW0g1qF+3vhSga\nupURaZjNsYV8L9D5avEajEjDTK8Be1e6DT+x6T+RNZ0G/gBaIHUs7KooXQAAAABJRU5ErkJggg==\n",
       "prompt_number": 5,
       "text": [
        "set([-1, 1, K, s, \u03c4])"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "params = H.atoms(sympy.Symbol)\n",
      "print(params)\n",
      "# Remove the Laplace variable\n",
      "params.remove(s)\n",
      "params"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "set([K, s, tau])\n"
       ]
      },
      {
       "latex": [
        "$$\\left\\{K, \\tau\\right\\}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAADkAAAAVBAMAAAAOWFv7AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZpkQzYnvq1QyRLvd\ndiJ+ofBJAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABLklEQVQoFXVSsUoDQRSc47JcNsmRYGWngtam\nsfYaCyuv0UqIrVaCNoLB/IHp7MTK0twXiH9wlhaK9wVGMSKiYX3vLbd7Qpxidt7Mvlve7gEIljAD\n+qnLbrgxIyNL77MfvBLFZ5MjRGaTjRLLLNQeczsB6resHHZZ2XSUAvMuELHj094V1op/0y1o6v6D\nSu8Ecy6KjTFTQNIoIbvxoX5cun158lIAC2ycD4nUNF0dcEXQA2S8NumoCxm69oZxwpZAH8py2gGu\nhyTbGZpfNiJWZAO1RaIoIxrdo/FJq0WLG+y56o5Uj+q8yx5jLOwnyql+ToH1hINHJjsR3WR4bPqI\n8+8CrXf2+0wulcLSjdf+y87rOGV75X1LL6yk8r6o/htBuY0u90B0/cFbXukVGvAXOYg9ziJ9AAgA\nAAAASUVORK5CYII=\n",
       "prompt_number": 6,
       "text": [
        "set([K, \u03c4])"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "H.subs?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "H_fixed = H.subs({K:2, tau: 0.001})\n",
      "H"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{K}{s \\tau + 1}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAADkAAAAsBAMAAAAtLQ2eAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdqvNEDJUuyJEiWaZ\n3e/xv6KKAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABLElEQVQ4EWNgQABG/c8qDEz/TREiKCx+AwYG\ntgoUISTO/AUMDNuQ+KjM9RMYvBpQhZB49gysC5C4aMzPDHvQRJC4nL+4/iJx0ZhcHxecT0ATQ3CZ\nfzP0A/2EA/AXMLB/xSHHwDBfgIHzE07Z9QEMDPcdcEnfB0rIL8Auy6L8X5GB8f73BuzSI0P0Px7w\nYcgFgQk+F7d/wSLL0gARdD+ETZZNAKqFkRjZhdIBDPLA8IRoQtPLHMB+gGtpU0QHxEQ02X4GVocw\nBnaobQxosv5vJgBlgnDIctZ/B8qsBstyGRvbPjY2PgDigN3MkcAi38DAcBUkAAKoJjM/YOBnYOCA\nRzSqLHsAgyDQmAsQneh6OQQ3JjAwMEE8i2Qyj95fbZgOBA03GSGExOLcgMQBACcVXk+CkCZNAAAA\nAElFTkSuQmCC\n",
       "prompt_number": 7,
       "text": [
        "   K   \n",
        "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n",
        "s\u22c5\u03c4 + 1"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Numerical Evaluation requires the `lambdify` function of SymPy (substitution fails miserably)\n",
      "\n",
      "see http://stackoverflow.com/questions/10678843/evaluate-sympy-expression-from-an-array-of-values"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Frequency vector\n",
      "freq = np.logspace(1,4) # Hz\n",
      "s_num = 2j * np.pi * freq\n",
      "f = sympy.lambdify(s, H_fixed, modules='numpy')\n",
      "H_num = np.abs(f(s_num))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plotting the Bode diagram"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.loglog(freq, H_num )\n",
      "plt.xlabel('freq (Hz)')\n",
      "plt.ylabel('gain')\n",
      "plt.title('Bode diagram of a first order system');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEfCAYAAABWPiGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVNX/P/DXgIgbEi6DIiBugYTiguUe9mk3yzQNUizU\nzAr9ZEW2/AxtUys/bvVxKbcozSyXlEQzw63SXPKjuRsgijqGSEIo2/n9cb9MoDPM3Jk7c++F1/Px\nmIece8898+Ye7xzmvO9iEEIIEBER2clD7QCIiEhfOHAQEZEsHDiIiEgWDhxERCQLBw4iIpKFAwcR\nEcnCgaOGycjIgIeHB8rKypxuKy0tDUFBQeZyREQEtm/f7nS7WrNr1y60a9cOPj4++Pbbb51uLz4+\nHo0aNUL37t2xc+dOhIWFKRClazz11FOYNGmS2mGQxtRSOwCyLSQkBCaTCZ6envDy8kLPnj0xf/58\nBAYGqh1aJYcPH1Y7BJd48803MX78eIwbN87ptnbs2IEtW7YgOzsbderUAQAcO3bMobYmT56M06dP\nIzk52em4rDEYDDAYDC5r3xkZGRlo3bo1SkpK4OHBv4HdiXtbBwwGAzZs2ICrV6/i/Pnz8Pf3V+RD\nTC9KS0tVff8zZ84gPDxckbYyMzMREhJiHjSqUlJSosh72svat1BHrhF2Z+y8htn9OHDojLe3NwYP\nHowjR46Yl+Xl5WHEiBEwGo0ICQnBu+++az6YysrK8PLLL6Np06Zo06YNUlJSKrWXl5eHUaNGISAg\nAIGBgZg0aZLVD5DCwkI89dRTaNSoEW677Tb8+uuvldaHhIRg69atAIA9e/agR48e8PPzQ0BAAMaN\nG4fi4mJz3c2bNyM0NBS33HILnn/+edx5551YtGgRAGDp0qXo1asXXnzxRTRp0gRTpkzBH3/8gbvu\nugtNmjRB06ZNMXz4cOTl5VV67w8//BAdO3aEj48PRo0ahYsXL+KBBx6Ar68v7rnnHly5csXqfv3k\nk0/Qrl07NG7cGI888gjOnz8PAGjTpg3++OMPDBgwAA0bNqz0O5SbNm0a2rZti4YNG+K2227D2rVr\nLb7HokWL8PTTT+Pnn3+Gj48PpkyZctN0X0hICN5//33z71FaWorp06cjMDAQDRs2RFhYGLZu3YrU\n1FRMnToVK1euhI+PDzp37mzxPY8ePYro6Gj4+fkhIiIC69evN6976qmn8Oyzz+LBBx9EgwYNkJaW\nhgMHDqBLly5o2LAhYmJicO3atUrtbdiwAZ06dYKfnx969eqFQ4cOWY3d0v+jCRMmwN/fH76+vujY\nsSN+//13/Prrr2jWrFmlAWD16tXo1KkTAOn/UlRUFHx9fdGsWTO8/PLLAIC+ffsCAG655Rb4+Phg\n9+7dAIDFixcjPDwcjRo1wv33348zZ86Y2/Xw8MC8efPQrl07NGzYEG+++SZOnz6NHj164JZbbkFM\nTIzFPqYbCNK8kJAQsWXLFiGEEAUFBWLEiBHiySefNK+Pi4sTAwcOFPn5+SIjI0PceuutYtGiRUII\nIebNmyfCwsLE2bNnxeXLl0V0dLTw8PAQpaWlQgghBg4cKMaOHSv+/vtvYTKZxO233y4WLFhgMY6J\nEyeKvn37itzcXJGVlSVuu+02ERQUVCnOH374QQghxL59+8Tu3btFaWmpyMjIEO3btxezZs0SQghx\n6dIl0bBhQ7FmzRpRWloqZs+eLby8vMwxL1myRNSqVUt89NFHorS0VBQWFopTp06JLVu2iKKiInHp\n0iXRt29f8cILL1R67x49egiTySTOnTsnjEaj6Ny5s/jtt9/EtWvXxF133SWmTJli8ff64YcfRJMm\nTcSBAwfE9evXxbhx40Tfvn0t/l6WrFq1Spw/f14IIcTKlStF/fr1zeUbLV26VPTu3dtc/vHHH0Vg\nYKC53LJlS9G5c2dx9uxZce3aNXHs2DERFBRkbi8zM1OcPn1aCCHE5MmTRVxcnNW4ioqKRJs2bcTU\nqVNFcXGx2Lp1q/Dx8RHHjx8XQgjx5JNPCl9fX/HTTz8JIYTIy8sTwcHBYtasWaKkpER8/fXXwsvL\nS0yaNEkIIcT+/fuF0WgUe/bsEWVlZWLZsmUiJCREFBUVWYz9RqmpqaJr164iLy9PCCHEsWPHzL9X\neHi42Lhxo7nuwIEDxX/+8x8hhBDdu3cXn3/+uRBC+v//yy+/CCGEyMjIEAaDwfx/WQgh1q5dK9q2\nbSuOHTsmSktLxTvvvCN69uxpXm8wGMTAgQPF1atXxe+//y5q164t+vXrJ9LT00VeXp4IDw8Xy5Yt\ns7pPScKBQwdatmwpGjRoIG655Rbh5eUlWrRoIQ4dOiSEEKKkpETUrl1bHD161Fx/wYIFIjo6Wggh\nRL9+/SoNBJs3bzYfbBcuXBDe3t6isLDQvH758uWiX79+FuNo3bq12LRpk7m8cOHCSh96VX3Azpw5\nUzz66KNCCCGWLVtW6WAWQoigoKBKA0dwcHCV+2TNmjWic+fOld57+fLl5vLgwYPFc889Zy7PnTtX\nDBw40GJbI0eOFBMnTjSX8/PzhZeXl8jMzLT5e1nSqVMnsW7dOovrlixZUuXAERISIpYsWWIunzx5\nUhiNRvOgWVFSUpIYPny41Ti2b98umjVrVmlZbGysmDx5shBCGjgq/gGybds2ERAQUKl+z549zQPH\n2LFjzT+XCw0NFdu3b7cY+422bt0qbr31VvHLL79U+rAXQohp06aJYcOGCSGEyMnJEfXq1RMXLlwQ\nQgjRt29fkZSUJC5dulRpm/T09JsGjvvvv9/8/0gIIUpLS0W9evXEmTNnhBDSwFE+UAohRNeuXcX7\n779vLr/00kuV/iAhyzhVpQMGgwHr1q1Dbm4url+/jrlz5+LOO++EyWTCn3/+ieLiYrRs2dJcPzg4\nGOfOnQMAnD9/vtJUSHBwsPnnzMxMFBcXo3nz5vDz84Ofnx/Gjh2LS5cuWYwjOzvbals3OnHiBB56\n6CE0b94cvr6+eOONN5CTk2Nu58bE/o3liu8DABcvXkRMTAwCAwPh6+uLuLg4c3vl/P39zT/XrVu3\nUrlOnTrIz8+3GOv58+cr7b/69eujcePG5n1oy2effYbOnTub9+Hhw4dvik2Oir9727ZtMWvWLEye\nPBn+/v6IjY01T6PZcmN/AUDLli2RnZ0NQPp/VXG/Z2dno0WLFjfVL5eZmYkZM2aYf08/Pz+cPXvW\n3N6Nsd+oX79+SEhIwPPPPw9/f38888wzuHr1KgBg2LBhWL9+Pf7++2989dVX6Nu3r7n/Fi1ahBMn\nTqB9+/a4/fbbb5purSgzMxP//ve/zfE1btwYACr1ZVX/T+rWrWv1/wn9gwOHzhgMBjz66KPw9PTE\nzp070aRJE3h5eSEjI8Nc58yZM+YPhObNm1ea4634c1BQELy9vZGTk4Pc3Fzk5uYiLy+v0rx1RVW1\ndaNnn30W4eHhOHXqFPLy8vDuu++a57wDAgJw9uxZc10hRKVy+e9Z0euvvw5PT08cPnwYeXl5SE5O\ntnlKsbAzaRoQEFBp/xUUFCAnJ+emD1FLMjMzMWbMGHz88ce4fPkycnNzERER4VTC9sbfPTY2Fjt2\n7EBmZiYMBgMmTpxosd6NAgICkJWVVSmWzMzMSr9XxTaaN29+02CZmZlp/jk4OBhvvPGG+f9Kbm4u\n8vPz8fjjj1uN/Ubjxo3D3r17ceTIEZw4cQIffPABAOkPh+7du2P16tX4/PPPERcXZ96mbdu2WL58\nOS5duoSJEyfiscceQ2FhocX3Cg4OxsKFCyvFWFBQgO7du1cZF8nDgUMnyg9+IYT520f79u3h6emJ\noUOH4o033kB+fj4yMzMxc+ZMDB8+HAAwdOhQzJkzB+fOnUNubi6mTZtmbrN58+a499578eKLL+Lq\n1asoKyvD6dOnrV6LMXToUEydOhVXrlzB2bNnMXfuXKvx5ufnw8fHB/Xq1cOxY8cwb94887oHH3wQ\nhw4dwrp161BSUoKPP/4YFy5cqPL3z8/PR/369dGwYUOcO3fO/IGjhNjYWCxZsgQHDx7E9evX8frr\nr6N79+5VfqMqV1BQAIPBgCZNmqCsrAxLlixR9LTkEydOYOvWrbh+/Tq8vb1Rp04deHp6AgCaNWuG\njIwMq4NU9+7dUa9ePbz//vsoLi5GWloaNmzYgJiYGAA3D6w9e/ZErVq1MGfOHBQXF2P16tWVToB4\n+umnMX/+fOzZswdCCBQUFCAlJcXuv9D37t2L3bt3o7i4GPXq1av0uwDAiBEjMH36dBw+fBiDBg0y\nL//888/N34J9fX1hMBjg4eGBpk2bwsPDA6dPnzbXHTt2LN577z3zySN5eXlYtWpVlXFV3A/ODPg1\niWYHjvT0dIwePRpDhgxROxRNGDBgAHx8fODr64tJkybhs88+Q/v27QEAc+fORf369dG6dWv06dMH\nw4YNQ3x8PADpYL/vvvsQGRmJqKgoDB48uNJfap999hmKiorMZ6EMGTLE6od4UlISWrZsiVatWuH+\n++/HiBEjrP6F+eGHH2L58uVo2LAhxowZg5iYGHPdJk2aYNWqVXjllVfQpEkTHD16FFFRUfD29gZg\n+dqBpKQk7N+/H76+vhgwYMBNv4clFddXdT3Cv/71L7z99tsYPHgwAgICkJ6eji+//LLKtsuFh4fj\npZdeQo8ePdCsWTMcPnwYvXv3rjKmG+Oo6ve4fv06XnvtNTRt2hTNmzfHn3/+ialTpwKA+dho3Lgx\noqKibtrWy8sL69evx8aNG9G0aVMkJCQgOTkZt956q8VYvLy8sHr1aixduhSNGzfGV199hcGDB5vX\nd+3aFZ988gkSEhLQqFEjtGvXDp999pnd13n89ddfGDNmDBo1aoSQkBA0adIEiYmJ5vWDBg3CmTNn\n8Oijj1Y6XXnTpk2IiIiAj48PJkyYgC+//BLe3t6oV68e3njjDfTq1Qt+fn7Ys2cPBg4ciIkTJyIm\nJga+vr7o0KEDNm3aVOW+tvf/Cf3DIDQ+xA4ZMsTmXwykb2VlZQgKCsLy5ctx5513qh0Oqahdu3ZY\nsGAB7rrrLrVDoSpo9hsHVW+bN2/GlStXcP36dbz33nsAwHnoGm716tUwGAwcNHTArQPHyJEj4e/v\njw4dOlRanpqairCwMLRr1w7Tp093Z0ikkp9//hlt27ZF06ZNkZKSgrVr15qnqqjmiY6OxnPPPYeP\nP/5Y7VDIDm6dqtqxYwcaNGiAESNGmM/cKS0tRWhoKLZs2YIWLVqgW7duWLFiBfz9/fH666/jhx9+\nwOjRo81nkhARkbrcepPDPn36VDrtEZBuJ9C2bVuEhIQAAGJiYrBu3Tq8+uqrmD9/vjvDIyIiO6h+\nd9xz585VumgoMDDQfM8ZW1q3bo309HRXhUZEVC21adMGp06dcnh71ZPjzpz6lp6eDiHdNkXVV1JS\nkuptydnOnrq26lhbL2e5kvtNC32nl/6Tu06rfafH/tPKsVfx2hdHqD5wtGjRAllZWeZyVlaW5p4z\nYUtkZKTqbcnZzp66tupYWy9neUFBgc04XE3JvnOmPXf2n9x1Wu07QH/9p5Vjz1mqT1VFRUXh5MmT\nyMjIQEBAAFauXIkVK1bYvX1iYiL69++P6OhomEwmAIDRaAQAt5V79eqlWHsVn/sgZ/tevXopGq+t\n9qyttxa/pfo+Pj4wmUxu7y8l9rfS7bmz/6ra3lL8Sv7/Vrqst/6zp31b+9va+vDwcIvHU8X6u3bt\nwsGDB+Est55VFRsbi23btiEnJwdGoxFvvfUW4uPjsXHjRrzwwgsoLS3FqFGj8Nprr9nVnsFggBvD\nJ4WlpaUhOjpa7TDIAew7fXP2s1PzV45XhQMHEZF8zn52qp7jqA7KvxKq2Zac7eypa6uOtfVyliu5\n3xyldAx66D+567Tad4D++k8rx56zVM9xOEsLOY5ySrSXm5urq3idiT83N9clv5+a+1sP/VfV9pbi\nd/b92H/y2i8nd72146lifV3mOJTGqSoiIvk4VUVERG7FgUMBzHE4tlwL8+R6myO3ty5zHO5tjzkO\nnWGOgzkOLe1vPfRfVdszx8Echz2Y4yAiqmGY4yAiIrfiwKEA5jgcW66FeXK9zZHbW5c5Dve2V9Ny\nHBw4iIhIFibHNVguJ2d7o9GoaDy22qtqvaX4LdUvr6PH/a10e+7uP2fjVzpeZ8r2xOvq9uTsD6Xj\nlXM8MTkOJseJiBzB5LgGMMfh2HItzJPrbY7c3rrMcbi3PeY4iIiIqsCpKiKiGoZTVURE5FY8q0qB\ncvkyJdrLzc1FaGio7O1vjMXZeG21Z229tfgt1T9+/Dj8/PxUPSvH0f2tdHvu7L+qtrcUv6X6cuJl\n/8lr/8Y27V1v7XiqWF+ps6ogdEwr4V+8eFH1tuRsZ09dW3WsrZezXMn95iilY9BD/8ldp9W+E0J/\n/aeVY8/Zz07mOIiIahjmOIiIyK04cCiA13E4tlwL1wLo7ToAe+vyOg73tsfrOIiIiKrAHAcRUQ3j\n7GcnT8dlmWWWWa4hZZ6OK3g6rqPbaeWUQC2c0qm30zntrcvTcd3bXk07HZc5DiIikoU5DiKiGobX\ncRARkVtx4FAAr+NwbLkWrgXQ23UA9tbldRzubY/XcRAREVWBOQ4iohqGOQ4iInIrDhwKYI7DseVa\nmCfX2xy5vXWZ43BvezUtx8ErxxUol1OivdzcXF3F60z8ubm5Lvn91Nzfeui/qra3FL+z78f+k9d+\nObnrrR1PFesrdeU4cxxERDUMcxxERORWHDgUwByHY8u1ME+utzlye+syx+He9mpajoMDBxERycIc\nBxFRDcMcBxERuRUHDgUwx+HYci3Mk+ttjtzeusxxuLc95jiIiIiqwBwHEVENwxwHERG5FQcOBTDH\n4dhyLcyT622O3N66zHG4t72aluPQ/b2qRowADAbpBVj+2cPjn58tLSv/2cPj5p9vLHt6Vv7XwwOo\nWxcoKZGW1aol/Vv+qlXr5peX1z8/164tlcv/vXYNKCqSyt7e/7zKfyciIrXpfuC4ejURoaH9ERYW\nDYPBBCGAsjIjhAA8PKRySYkRZWWAp6dULi6W1nt6SiPx9evSei8vaX1hobS+dm2pXFAgrff2lsp/\n/SWV69c3obQUOHNGKjdoYEJZGfDnn0aUlgK33CKtv3jRiJISoFEjqXz2rFRu2tSEkhIgPd2I4mIg\nIEBaf+yYNHi0bGlCcTHwv/8Z4eUFdOxogpeXtH3dukBYmAne3sDly0bUrWtEcLAJdetKv1+DBkCT\nJlLZ21sq169vQoMGQLNmxv/bDybUrm35pmxGo7HKm7ZVtb6crfrlddS+aZ6leN3dnq397c727Ilf\n6XidKdsTr6vbk7M/lI5XzvHEmxyi5iTHhZAGkuvXpW8k164BhYXA33/f/G9BAZCff/Pr6lUgLw+4\ncqXyy8MD8PMDmjT559W06T8/G41AQMA/r3r11N4bROQsZz87OXAooOIor1ZbcrYrryuENAhdvgz8\n+Wfl199/m3DmjBEXLgDnz0uv7Gxp2iwgAOja1YRatYwICQFatZJeISFArVomNG9+cxyW4lNyvzlK\n6Rjc2X+O1pG7Tqt954o4XN1/zvZdVevlLHf2s1P3U1XkOINBys+0aCG9KjKZpG8bFQkhfUvJzgYy\nM6XBJD0d2LxZ+jcjQ/q2kp8P3Hor0L699AoLA/z9b26PiPSJ3zhIUdeuSYPK8ePA0aNSvuboUelV\nuzYQHg506gR06SK92reXThIgIvfhVJV+w69RhJC+oRw5Avz2G7B/v/TKygIiIqRBpGtXoFcv6RsK\nzyIjch1eAKgBvI7D9vJLl0wICADuvht4+WVg+XJg+3YTLlwAZsyQvnls2wY8+KA0pTVwoLR8926g\nuNj27+QovV0HYG9dXsfh3vZ4HQeRG/n4AL17S69yZ88CO3dKr+Rk4PRp4I47gPvvl1633cZvJERq\n4lQVad6VK9K3kU2bgNRU6dTk++6TBpG775ZOJyYi+zHHod/wyQFCAKdOSQNIaiqwY4eUbH/sMWDQ\nICAwUO0IibSPOQ4NYI7DseWO/K4GA9CuHTBuHJCSIp02/MorUqK9Y0egZ08pN5KRYV97epsjt7cu\ncxzubY85DiIdqVMHeOgh6VVUBGzdCnzzDdCtm3RBYmwsMGyYdB0JESmDU1VULZWUAGlpwOefA2vX\nAn37Ak8+KQ0w3t5qR0ekLuY49Bs+uUl+vvQtZOlS4PBhICZGGkS6duXZWVQzMcehAcxxOLbcXfPk\nDRpIA8WPPwJ79kg3bxw6FIiKApYuNaGwULn30kP/McehfHs1LcfBgYNqlFatgKQk6cysd96RzsoK\nDgYSE4E//lA7OiJ90OxUVUFBAZ577jl4e3sjOjoaTzzxxE11OFVFSvjjD2DePGkq6/bbgYQE6ToR\nD/5ZRdVUtc1xJCcno1GjRujfvz9iYmLw5Zdf3lSHAwcpqbAQ+PJLYM4c6QytxETgiSekmzMSVSfV\nNsdx7tw5BAUFAQA8PT1VjqZqzHE4tlwL8+QVY6hbF4iPl64JmT1bup9W69bAhx8Cf/0lvz1H41Ci\nLnMc7m2POQ4XGjlyJPz9/dGhQ4dKy1NTUxEWFoZ27dph+vTpAIDAwEBkZWUBAMrKytwZJtVwBoN0\nK5PNm4H166WBpFUr4NVXpTv8EtV0bp2q2rFjBxo0aIARI0bg0KFDAIDS0lKEhoZiy5YtaNGiBbp1\n64YVK1agZcuWSEhIQJ06ddCnTx/ExsbeHDynqshN0tOBmTOl60KGDZMGkRsffkWkF7qaqurTpw/8\nbrgj3Z49e9C2bVuEhITAy8sLMTExWLduHerVq4fFixfjv//9r8VBg8idWrWSch9Hj0pXq3foIN32\n5Nw5tSMjcj/VbzlSMZcBSFNUu3fvtnv7xMRE1K9fHwAQGRmJXr16mZ+vWz635+py+TIl2svNzUVo\naKjs7W+Mxdl4bbVnbb21+C3VP378OPz8/NzeX87sb4MB+OADIxITgfnzTRg8GIiKMuLVV4HatfXR\nf1Vtbyl+S/XlxKul/nNVe/buD3vav7FNe9dbO54A4MiRI0hJSQEA8+elU4Sbpaeni4iICHP566+/\nFqNHjzaXk5OTRUJCgl1tqRC+RRcvXlS9LTnb2VPXVh1r6+UsV3K/OcrZGC5eFCIxUYhGjYRISBDi\n6FHt95/cdVrtOyGUj8PVx59Wjj1nPztVP6uqRYsW5iQ4AGRlZSFQZ/fGLh/h1WxLznb21LVVx9p6\nOcuV3G+OcjYGoxF4/31pCsvLC+jVy4j/9/+AvDzXxeFs/8ldp9W+A5SPw9XHn1aOPWepPnBERUXh\n5MmTyMjIQFFREVauXImHH35Y7bCIZDEagf/8BzhwAMjOlm79/uGHUPR2JkRa4dYcR2xsLLZt24ac\nnBwEBQXhrbfeQnx8PD766CPcd999KC0txahRo9C+fXu720xMTET//v0RHR2tyhxrOeY4akaOw1a5\nsDAXixeH4sgRYPZsEx54ABg+3IinngIuX9ZG/1W1PXMc1TPHUV5/165dOHjwIJzm1ESXyrQSPnMc\nji3Xwjy5q+fIf/5ZiOhoIcLChNiwQYiyMufjYI7jH8xxOLbc2c9Ozd5yxB68joP0QAjpaYWJidK1\nHzNmAJGRakdFNZmzn52qn47rLC1MVbHMsq3yQw8BXbqYsGEDcN99RvTvD7z0kglNmmgjPpZrRplT\nVYJTVY5up5Wvy1qY7lBjquPKFSEmTpRO4Z08WYj8fE5VOYpTVY4td/azU/WzqohqGl9fYNo0YN8+\n4NgxoH176VnpnHUlvWCOg0hlO3cC48dLTyqcPRvo3FntiKi609W9qojoZr17A7/+CsTFAQ88ADzz\nDHDpktpREVnH5LgC5fJlejiP3N54bbVnbX1Nv47D0fZyckx45BHgrruAuXONCA8H3npLWhYQoGz/\nVbU9r+PgdRx2cSpDojKthM/kuGPLtZBg1Wpy9fffhbj7biEiIoTYtk3+ezA57t72alpynDkOIo0S\nAvjmG2DCBCA6GvjgA6BZM7WjouqAOQ6iaspgAB57TLqBYosW0jNAZs8GSkrUjoxqOg4cCqg4l6hW\nW3K2s6eurTrW1stZruR+c5TSMbii/xo0kE7f3b5depTtwIEm7NjheHty12m17wB99J/ceu449pyl\n+4EjMTERaWlpAKQdVHEn6bGcm5urqXhcGX9ubq6u4lW7/xo3NuGLL0yIiwOeeEK68vzECW3vD/af\n+9q353has2YNJk+eDGcxx0GkQ1evAklJ0jPQ33sPGDkS8ND9n4HkLs5+dnLgINKxAweAZ58FPD2B\n+fOlPAiRLUyOa0DFr4NqtSVnO3vq2qpjbb2c5UruN0cpHYO7+69zZ+Cnn4ARI6RrQBITgfx85jjc\n3Z6922nl2HMWBw4infPwkK42P3wYuHABCA8Hdu1SOyqqzjhVRVTN/PADMHYs0LEjMGeOdCovUUU1\nfqqqup1VxTLLzpb/9S/g0CGgWzcThgwx4eOPgdJS7cTHsnplpc6q0sY9OxyklfB5yxHHlmvhthV6\nu2WFvXXL6/z+uxC9eglxxx1CHDxoe3vecoS3HLGH7r9xEJF14eHShYOjRwN33w1MnAhcu6Z2VKR3\nzHEQ1RAXLwIvvCDdwn3+fGkgoZqJ13HoN3wiVaSkAM89B/TrB8yYATRurHZE5G41PjmuBRUTUGq1\nJWc7e+raqmNtvZzlSu43Rykdgx76r1s3Ew4flh5hGxEBLF/+z2Nr7e0nLfQdoL/+08qx5yybA8fO\nnTtxzz33oF27dmjVqhVatWqF1q1bKx4IEbmPj490p921a6UbKPbvD2Rmqh0V6YXNJwCOGjUKs2bN\nQpcuXeDp6emOmGTRwhMAXfFEQbnbK/UEQnvbs/WEMnvql9fR4/5Wuj139195uVUrYN8+Iz74AIiJ\nMWHECGDMGOkWJu6M1937W+n25OwPpeOVczwp9QRAmzmOO+64A7t373b6jVyBOQ4i5Rw/Djz9NFBc\nDHz6KXDbbWpHRK7i8hxHv379kJiYiJ9//hn79+83v+gfzHE4tlwL8+R6myO3t25VdaytCw0FvvrK\nhCeflJ7P2z3tAAAXDUlEQVQ4OHkycP26dvsO0F//aeXYc5bNqapffvkFBoMBe/furbT8xx9/VDwY\nIlKXh4d0u5KHHpLOvOrSBfjoI6DCbAgRT8clIsuEAFatAv79b2DoUODdd6WnEZL+uew6juTkZMTF\nxWHGjBkwGAzm5UIIGAwGvPjiiw6/qVI4cBC5Xk4O8OKL0hXoCxcC99yjdkTkLJflOP7++28AwNWr\nVy2+6B/McTi2XAvz5HqbI7e3riM5DmvrSktNWLYMmDdPunXJyJHA6dPq9x2gv/7TyrHnLKs5jmee\neQYAlLmTIhHp3v33S8/8eP11ID4emDABePRRtaMiNdjMcRQWFmLRokU4cuQICgsLzdNWixcvdkuA\nVeFUFZE6du6Uvn106CAlz/391Y6I5HD56bhxcXG4ePEiUlNTER0djaysLDTQUIaMz+NgmWX3l3v3\nBrZsMaFTJxM6dgSSk4GLF7UTH8uWy257HkdkZKQQQogOHToIIYQoKioSt99+u1P3cleKHeG7BZ/H\n4dhyLTzTQW/Pc7C3rtxnblS1ztayffuEiIwU4oEHhDhzxmZoitJb/2nl2HP2s9PmN47atWsDAHx9\nfXHo0CFcuXIFly5dcn7EIqJqoUsX6VbtPXtKPy9YAJSVqR0VuZLNHMenn36KQYMG4dChQ4iPj8fV\nq1fx9ttvY+zYse6K0SrmOIi05cgR6ayrunWl25a0aaN2RGSJy5/HUX4dR8Vqvr6+iIqKQqdOnRx+\nYyVw4CDSntJS6c67770HvPEGMH68dNNE0g6XJ8f37duH+fPnIzs7G9nZ2Vi4cCFSU1Px9NNPY/r0\n6Q6/cXVSMQGlVltytrOnrq061tbLWa7kfnOU0jHoof/krpPbd56e0gWDv/wCrFsnJdKPHrURsIP0\n1n9aOfacZXPgyMrKwv79+zFjxgzMmDED+/btg8lkwrZt27B06VLFAyKi6qFtW2DrVmDECKBvX2Dq\nVKCkRO2oSAk2p6rCwsLwv//9z5wkv379Ojp27Ijjx4+jc+fOOHDggFsCtYRTVUT6kJkp3bL98mVg\n8WKgY0e1I6rZnP3stHl33GHDhuGOO+7AwIEDIYTA+vXr8cQTT6CgoADh4eEOvzER1RwtWwKbNgFL\nlgB33y3deff114H/+3uUdMbmVNWkSZOwcOFC+Pr6ws/PDwsWLEBSUhLq16+PL774wh0xah5zHI4t\nZ47Dse20nuOwxmCQzrg6cADYtw+IipL+dYbe+k8rx56zbH7jAIBu3bqhW7duir85EdU8LVoA334L\nfPEF8OCDwKhRwJtvAnXqqB0Z2YvP4yAi1Vy4ADz/vHTW1eLFQPfuakdUM7j8Og4t48BBpH/lD4wa\nPx4YPhx4+23pAkJyHZcnx7UuMTER/fv3R3R0tHkuz/h/z7l0V7l8mRLt5ebmIjQ0VPb2N8bibLy2\n2rO23lr8luofP34cfn5+bu8vJfa30u25s/+q2t5S/Jbqy4nXnvLQoUCnTibMnQtERhqxaBEQGlr9\n+s+e9i3tb3vWWzueKtbftWsXDh48CKc5dacrlWklfN7k0LHlvMmhY9tp7SaHSlu9WojmzYUYP16I\n/Pyq6+qt/7Ry7Dn72cmpKiLSnJwc4IUXgJ9+AhYtAqKj1Y6oemGOQ7/hE5ENGzYAY8cCDz8MTJ8O\n+PioHVH14PJ7VZFtvI7DseW8jsOx7fR6HYcjHnoIOHQIKCyUnja4ZYtr4+B1HPbhwEFEmubnJ11x\nPm+edAHhM88Af/2ldlQ1G6eqiEg38vKAl18GNm8GPvkEuPdetSPSJ+Y49Bs+ETno+++lmybefTcw\nYwbg66t2RPrCHIcGMMfh2HLmOBzbriblOKy55x4p99G0qQkREcB33ynTLnMc9uHAQUS65OMDTJgA\nLFsGJCQA8fFAbq7aUdUMnKoiIt3LzwdefRVYuxaYP186G4usY45Dv+ETkcLS0qS77fbqBcyaBTRq\npHZE2sQchwYwx+HYcrXnyV0Rgx76rzrkOMrdGEd0NPC//wG33CJd9/Htt861p/R2Wjn2nMWBg4iq\nlfr1gTlzgBUrgBdfBIYNk25hQsrhVBURVVt//y09ovarr4CPPwYefVTtiLSBOQ79hk9EbrJzp3TV\nedeuwNy5QJMmakekLuY4NIA5DseWa2GenDkO2+u02neA/XH07g389hsQECDlPr75xrn2HN1OK8ee\nszQ7cKSnp2P06NEYMmSI2qEQUTVQr550lfk330jTV48/Dly6pHZU+qT5qaohQ4Zg1apVFtdxqoqI\nHFFYCLz5JvD551Iivab9fcqpKiIimerWBT74AFizBpg0CRg6FNDI7JsuuHzgGDlyJPz9/dGhQ4dK\ny1NTUxEWFoZ27dph+vTpAIDk5GRMmDAB2dnZrg5LUcxxOLZcC/PkzHHYXqfVvgOcj6N7d+DAAaBV\nK6BjR2DlSuY47OHygSM+Ph6pqamVlpWWliIhIQGpqak4cuQIVqxYgaNHjyIuLg4zZ85EQEAALl++\njLFjx+K3334zDyxEREqrW1d6uuC6dcDSpdK0lUbGRc1yS44jIyMDAwYMwKFDhwAAP//8M6ZMmWIe\nUKZNmwYAePXVV2W1yxwHESnp2jVg8mRpAJk9W5rCMhjUjkp5zn521lIwFrudO3cOQUFB5nJgYCB2\n797tUFuJiYmoX78+ACAyMhK9evWC0WgE8M9XNJZZZplle8vTphkxaBCQlGTCjz8CU6YY4e+vnfgc\nKaelpSElJQUAzJ+XThFukJ6eLiIiIszlr7/+WowePdpcTk5OFgkJCbLbdVP4Nl28eFH1tuRsZ09d\nW3WsrZezXMn95iilY9BD/8ldp9W+E8K1/VdYKMSrrwphNAqxYoUQZWXOx6GVY8/Zz05Vzqpq0aIF\nsrKyzOWsrCwEBgaqEQoRkUV16gBTpwLr1wNvvw089hhw8aLaUWmDKlNVUVFROHnyJDIyMhAQEICV\nK1dixYoVDrWVmJiI/v37Izo6WlNfDZ0pl5OzvdFoVDQeW+1Vtd5S/Jbql9fR4/5Wuj1395+z8Ssd\nrzNle+J1pr3bbzdi3z5gxgwTBg0CEhKMiIkBLl1ybH8oHa+c42nXrl04ePAgnOXy5HhsbCy2bduG\nnJwcGI1GvPXWW4iPj8fGjRvxwgsvoLS0FKNGjcJrr70mu20mx4nInX79FXjqKeDWW4F584BmzdSO\nyDFOf3Y6NdGlMq2EzxyHY8u1ME/OHIftdVrtOyHU6b9r14R4/XUp9/HFF1LugzkOIiKyytsbePdd\nICUFeO89YNAg4PJltaNyL1VyHEpijoM5Dq3PkWux/5yNvyblOKzVDw42ITUVmDfPiDvvNGLmTBP+\n9S/A3585Dk1jjoOItGDvXiA+HmjTRsp9NG+udkRV400ONeDGvyTUaEvOdvbUtVXH2no5y5Xcb45S\nOgY99J/cdVrtO0A7/RccbMLevUBEBNCpE/DFF4Clz2WtHHvO4sBBRKQAb2/gnXeA774Dpk0DBg4E\nzp9XOyrXYI5Dg+VyWp4jZ45D3/3nbPzMcVjfH127GrF3LzBzpgmDBwPPPWfEsGH/XPehdLzMccjE\nHAcRadm+fVLuo1UrYP587eQ+mOPQAOY4HFuuhXlyrcyRM8fhGK33X9euUuK8Y0cgMhL4/HOTxdyH\nnBiY4yAiquZq15budZWaCnz5JfDII/rPfXCqiojITYqKpAT6/PnAjBnA8OHqPO+jxk9VJSYmIi0t\nDYD0lazi1zKWWWaZZS2Vr1wx4a23pG8fX31lwtNPm1D+pGx3vP+aNWswefJkOM2pG5aoTCvh815V\nji3Xwv2OeK8q2+u02ndC6K//Kta7fl2IN98UomlTIZYt++d5H7xXFRERWVS7NjBlCrBpkzRt9fDD\nMH/70DrmOIiIVFZUJN04cd484MMPgbg41+Y+nP3s5MBBRKQRBw5Iz/sIDgYWLAACAlzzPjU+Oa4F\nFRNQarUlZzt76tqqY229nOVK7jdHKR2DHvpP7jqt9h2gv/6zVa9zZ+C770zo2lW659WyZTff80qJ\nY89Zuh84qttZVbm5uZqKx5Xx5+bm6ipe9p/2ynrrP3va9/ICJk8G1q414ZtvTBgwADh3TpnjSamz\nqjhVRUSkUUVF0sOi/vtf4IMPgBEjlMl9MMeh3/CJiOzy229S7iMwUMp9tGjhXHvMcWhAxa+DarUl\nZzt76tqqY229nOVK7jdHKR2DHvpP7jqt9h2gv/5ztO86dQL27AGiooAhQ0yK5D6cwYGDiEgHateW\nch8ffgjMnAlz7kMNnKoiItIZZ3MfzHHoN3wiIqc4mvtw9rOTTwBUoFy+TIn2cnNzERoaKnv7G2Nx\nNl5b7Vlbby1+S/WPHz8OPz8/VZ8g5+j+Vro9d/ZfVdtbit9SfTnxsv/ktX9jm1WtDwgA9uwxYupU\n4IUXjuPhh/0wfLgRBoPl+ko9AVAbdwl0kFbC500OHVuuhRvl6e0mefbW5U0O3dueIzc5dLSOtfW7\ndl0UkZFCPPigEGfPVl3f2c9OTlUREVUTRUXA1KnAxx8D778PPPmk5dwHcxz6DZ+IyCV++0161nlA\nALBw4c25D17HoQG8jsOx5Vq4FkBv1wHYW5fXcbi3PVdfx2HP+orLy6/7uP126bqPpUtvvu7DGRw4\niIiqIS8vIClJuu5j9mzgoYeUu+6DU1VERNVccbF03Ud57iM+njkOtcMgItKFgwel6z5++405DtUx\nx+HYci3Mk+ttjtzeusxxuLc9reU4rC2PjJRyH87iwEFEVIN4eTnfBq8c12C5nJztlbpy3d72bF3Z\nak/98jp63N9Kt+fu/nM2fqXjdaZsT7yubk/O/lA6XjnHk1JXjjPHQURUw/A6Dg1gjsOx5VqYJ9fb\nHLm9dZnjcG97eslxKIUDBxERycKpKiKiGoZTVURE5FYcOBTAHIdjy7UwT663OXJ76zLH4d72mOMg\nIiKqAnMcREQ1DHMcRETkVhw4FMAch2PLtTBPrrc5cnvrMsfh3vaY4yAiIqoCcxxERDWMs5+dvMkh\nyyyzzHINKSt1k0MIHdNK+BcvXlS9LTnb2VPXVh1r6+UsV3K/OUrpGPTQf3LXabXvhNBf/2nl2HP2\ns5M5DiIikoU5DiKiGobXcRARkVtx4FAAr+NwbLkWrgXQ23UA9tbldRzubY/XcRAREVWBOQ4iohqG\nOQ4iInIrDhwKYI7DseVamCfX2xy5vXWZ43Bve8xxEBERVYE5DiKiGoY5DiIicisOHApgjsOx5VqY\nJ9fbHLm9dZnjcG97zHEQERFVgTkOIqIapto+j2PdunVISUnBX3/9hVGjRuGee+5ROyQiIoKGp6oe\neeQRLFy4EPPnz8fKlSvVDqdKzHE4tlwL8+R6myO3ty5zHO5tjzkOjXnnnXeQkJCgdhjkArt27VI7\nBHIQ+65mc/nAMXLkSPj7+6NDhw6VlqempiIsLAzt2rXD9OnTAQDJycmYMGECsrOzIYTAxIkT8cAD\nD6BTp06uDtMp5Y9mVLMtOdvZU9dWHWvr5SxX5BGWTlKy75xpz539J3edVvsO0F//aeXYc5bLB474\n+HikpqZWWlZaWoqEhASkpqbiyJEjWLFiBY4ePYq4uDjMnDkTAQEBmDt3Ln744Qd8/fXXWLBggavD\ndEpaWprqbcnZzp66tupYWy93udqUjksP/Sd3nVb7DtBf/1WXY8/lA0efPn3g5+dXadmePXvQtm1b\nhISEwMvLCzExMVi3bl2lOuPHj8fevXsxb948PPPMM64O0ykpKSmqtyVnO3vq2qpjbb2c5QUFBTbj\ncDUl+86Z9tzZf3LXabXvAP31n1aOPac59cRyO6Wnp4uIiAhzedWqVWL06NHmcnJyskhISJDdbps2\nbQQAvvjiiy++ZLzatGnj1Ge6KqfjGgwGRdo5deqUIu0QEZH9VDmrqkWLFsjKyjKXs7KyEBgYqEYo\nREQkkyoDR1RUFE6ePImMjAwUFRVh5cqVePjhh9UIhYiIZHL5wBEbG4uePXvixIkTCAoKwpIlS1Cr\nVi189NFHuO+++xAeHo7HH38c7du3d3UoRESkAF3fq4qIiNxP81eOy5Geno7Ro0djyJAhaodCDli3\nbh3GjBmDmJgYfP/992qHQzIcO3YMzz77LIYOHYpFixapHQ45oKCgAN26dbPr9N1q+Y1jyJAhWLVq\nldphkIOuXLmCl19+GZ9++qnaoZBMZWVliImJwVdffaV2KCRTUlISfHx80L59e/Tv37/KutXqGwdV\nD7w/mT6tX78e/fv3R0xMjNqhkEzff/89wsPD0bRpU7vqa37gkHOvK9IeOf0ndHR/sppA7rE3YMAA\nbNy4EcuWLXN3qGSBnP7btm0bfvnlFyxfvhyffPKJ7Wd1OHX5oBts375d7N+/v9KV5yUlJaJNmzYi\nPT1dFBUVicjISHHkyBGRk5MjnnnmGdG2bVsxbdo0FaOmcnL6b86cOaJr165i7NixYv78+SpGTULI\n67u0tDQxfvx4MWbMGDFz5kwVo6Zycvqv3NKlS0VKSorNtjU/cAhx8y1LfvrpJ3HfffeZy1OnThVT\np05VIzSyA/tPv9h3+uaq/tP8VJUl586dQ1BQkLkcGBiIc+fOqRgRycH+0y/2nb4p1X+6HDiUutcV\nqYP9p1/sO31Tqv90OXDwXlf6xv7TL/advinVf7ocOHivK31j/+kX+07fFOs/RTIwLhQTEyOaN28u\nateuLQIDA8XixYuFEEJ899134tZbbxVt2rQR7733nspRkjXsP/1i3+mbK/uvWl45TkRErqPLqSoi\nIlIPBw4iIpKFAwcREcnCgYOIiGThwEFERLJw4CAiIlk4cBARkSwcOIgAzJkzB+Hh4YiLi3NJ+3ff\nfTeuXr0KAGjQoEGldUuXLsW4ceOsbvvtt9/i7bffdklcRI7gwEEEYN68ediyZQuSk5MrLS8pKXG6\n7a1btyI0NBQ+Pj4Abr7RnK0bzw0YMADffPMNiouLnY6FSAm11A6ASG1jx47FH3/8gfvvvx8jR45E\nXl4eTp06hfT0dLRs2RKzZ8/G2LFjcebMGQDArFmz0LNnT+Tk5CA2NhbZ2dno0aMHvv/+e+zfvx+N\nGjWq1P7y5csxdOhQq+9f8eYNnTp1Mg8kx48fx6ZNm9CnTx/06NEDmzdvtvksaCK3UOzGKEQ6FhIS\nInJycoQQQiQlJYmoqChx7do1IYQQsbGxYufOnUIIITIzM0X79u2FEEKMGzdOvP3220IIIVJSUoTB\nYDC3UVFYWFil5Z6enqJTp07mV3BwsBg3blylbb799lvRt29fUVJSIoQQYvHixeKVV15R+Lcmcgy/\ncRDdwGAw4OGHH4a3tzcAYMuWLTh69Kh5/dWrV1FQUIAdO3ZgzZo1AIAHH3wQfn5+FtvLzs6u9C2k\nbt26OHDggLm8bNky7N2711w+efIkXnnlFaSlpcHT0xMAEBAQgNTUVOV+SSIncOAgsqBevXrmn4UQ\n2L17N2rXrn1TPaHAPUIrtpGfn4/HH38cn376Kfz9/c3Ly8rK+BAl0gwmx4lsuPfeezFnzhxz+eDB\ngwCAvn37Yvny5QCAjRs3Ijc31+L2AQEByMnJseu9Ro4cifj4ePTq1avS8vPnz6Nly5aOhE+kOA4c\nRKj6TKc5c+Zg7969iIyMxG233YYFCxYAAJKSkrB9+3ZERERgzZo1CA4Otth27969K01FWXovg8GA\nM2fO4JtvvsHixYvRuXNndO7cGfv37wcA7NmzB3379lXkdyVyFp/HQaSQVq1aYd++fTedVZWWloaV\nK1di3rx5DrVbVlaGLl26YO/evahVi7PLpD5+4yBSiLUcRHR0NE6ePGm+AFCuDRs24LHHHuOgQZrB\nbxxERCQLv3EQEZEsHDiIiEgWDhxERCQLBw4iIpKFAwcREcnCgYOIiGT5/7s876F69u0CAAAAAElF\nTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x4b812d0>"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Building a Traits object around a SymPy expression"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import traits.api\n",
      "from traits.api import HasTraits, Instance, Float"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A Traits class which can contain a Sympy expression\n",
      "\n",
      " * see http://docs.enthought.com/traits/traits_user_manual/defining.html#this-and-self for the `Instance` traits\n",
      " * see http://docs.enthought.com/traits/traits_user_manual/notification.html#specially-named-notification-handlers for name-based notification"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class ExprTrait(HasTraits):\n",
      "    expr = Instance(sympy.Expr)\n",
      "    \n",
      "    def _expr_changed(self, old, new):\n",
      "        print('new expression: {:s}'.format(str(new)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 54
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m = ExprTrait(expr = H)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "new expression: K/(s*tau + 1)\n"
       ]
      }
     ],
     "prompt_number": 55
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m.expr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{K}{s \\tau + 1}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAADkAAAAsBAMAAAAtLQ2eAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdqvNEDJUuyJEiWaZ\n3e/xv6KKAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABLElEQVQ4EWNgQABG/c8qDEz/TREiKCx+AwYG\ntgoUISTO/AUMDNuQ+KjM9RMYvBpQhZB49gysC5C4aMzPDHvQRJC4nL+4/iJx0ZhcHxecT0ATQ3CZ\nfzP0A/2EA/AXMLB/xSHHwDBfgIHzE07Z9QEMDPcdcEnfB0rIL8Auy6L8X5GB8f73BuzSI0P0Px7w\nYcgFgQk+F7d/wSLL0gARdD+ETZZNAKqFkRjZhdIBDPLA8IRoQtPLHMB+gGtpU0QHxEQ02X4GVocw\nBnaobQxosv5vJgBlgnDIctZ/B8qsBstyGRvbPjY2PgDigN3MkcAi38DAcBUkAAKoJjM/YOBnYOCA\nRzSqLHsAgyDQmAsQneh6OQQ3JjAwMEE8i2Qyj95fbZgOBA03GSGExOLcgMQBACcVXk+CkCZNAAAA\nAElFTkSuQmCC\n",
       "prompt_number": 56,
       "text": [
        "   K   \n",
        "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n",
        "s\u22c5\u03c4 + 1"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Adding traits corresponding to each parameter**\n",
      "\n",
      "see http://docs.enthought.com/traits/traits_user_manual/advanced.html#per-object-trait-attributes for dynamic initialization of traits using  `add_trait`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Add traits manually:\n",
      "for p in params:\n",
      "    print(p)\n",
      "    m.add_trait(p.name, Float(1))\n",
      "\n",
      "# for a mysterious reason, the newly added traits need to be visited:\n",
      "m.K\n",
      "m.tau\n",
      "\n",
      "# List the traits\n",
      "m.print_traits()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "K\n",
        "tau\n",
        "K:    1\n",
        "expr: K/(s*tau + 1)\n",
        "tau:  1\n"
       ]
      }
     ],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Look at all the traits:\n",
      "m.get()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\begin{Bmatrix}K : 1, & expr : \\frac{K}{s \\tau + 1}, & tau : 1\\end{Bmatrix}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAARkAAAAcCAMAAABidDsiAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAMolUmc0idhCr\nu+9m3UQgGr1z5AAAAAlwSFlzAAAOxAAADsQBlSsOGwAABGlJREFUWAnNWYm2oyAMRUA2geH/v3YS\nEAhqfTp9dco5rWHLcglJahn78hZdcgJ0nJMzX67q0+rpBSUqq54WnOVNHI/lO5t1oJfQjykX5dRl\nyfn/HEjX4IRaABQdTxb89pSa8Sxy0+mLb7BJ0ThfVX3mufBVjp27QG1Tmi0ihQTxq74EouGDSIoU\neT9EqsUv08QoVwGhyDCW1mHhXlwx5Z70MZ48i+nFEZ2AY8IL9Y/3DEY1QBqBm2IqruRfxGTjpH0S\nmVmCUsEe23MyKtLJ5HZqY1QDpBG4gSeMdoafIP5oXEp4QiKdqINa75usN2I/dThCjWqANAK3zIj1\nVEPQj0wOF/zi4FQwWdBzbrVwasKe1c/IYJjxWz1Gb6ZM9iLKiOFaSzxoo7VVUXPIvZE7NUHXAGVN\nhBVt9yigDTMeEnIRS9qq1NccUNFC5WyRvbJcTwEyhuEy95eePQah1KjmKq7l7xJmeAqjNLMMrkmZ\njAtbL2L8U3BsBj4+aBaAA2c2wN3gDig3C+ZydYt7NgIam38n1jATF9DD4TWQAC88dA8/o1BqVEXG\nYC1VG0+aM5kvdx3aPSmT3WQeUAsgYLCY5HBEHCK2UEwJNmMk9WkSOabG+xnnWNzBaAkzZsFKSMLJ\nK8/yEPWCYRs1ypcUo3L0r6tmSJJM1cxdR8cnZTLO1J5dvNYe/RY/azkCZEacpz8M6re6dvdM19pu\n3zBQwozMJoLLoh7ZAagXDBsGo3iAQ/Nj3k/5ZrnTAmJgMrCvnTRk2VR9MmZF8RLdyqmV643nivyC\n1phiTBY+vTRsNMrMoLOYSe0dixVxNG2j0shkM5m7a01Upoo26Dwc45fBYHo3pxZO6/dE4uIw0TsR\nw4kxWQ88BZAuUTgJM31xpgajVkwUiSp8xTSUm7bZvXYHJodL1h8eiLsp2mhEJt9bj9mh/TI53H44\nWK/fJO0mQRws5xDwxWSyZVDrQ5RjDv0Ywky+5Pst1KhWc5OoFNbQrYdTH6M4mIpmwuEPGYsK0zhj\nIN9G0A1TAwRAaHhVpwWusNqEmY0AyqrR3VHEisyJAtJCuQBHAR8PyRlrG4zDEeqiVueMQikyNTex\nShg7p5SrAB1Soj+baOaXdknBoqUnjqWl9rjEWK2N5pkGmKCGyUVO3J46FdCwGIk9MicKKJeLJRDP\nJ2XzSwwDlY2YrO4JkQjtRqHQCkgnRlV+7kVw0hsth5lL66EIhLSv4HygrdfoABl2U4FLsmHR+8iQ\n0H1F6FAenG2QkWGUlZOxJl9cXHyEzE0FzmTSubeRMfcUK/GQavCKdk4AIohJh4OQNc6wmwq8Ercb\nfxuZe8CIkDAgXmkTvDXL4SCWcGlnaAt+ZZENmXsKXJFc1ryNzHVR91Zijs2JgElie3efhsw9ttdX\nt3d6vP2yu775kyuxNC8v8JGqjSDzaX1D+2FfXvxWFf7702u//llAk3tFRslQSouP6WnxrcXaRP4b\nsPa+81mR+bR28el/Kt42qHn425yuMfgLSzUftkwmX2IAAAAASUVORK5CYII=\n",
       "prompt_number": 60,
       "text": [
        "\u23a7               K           \u23ab\n",
        "\u23a8K: 1, expr: \u2500\u2500\u2500\u2500\u2500\u2500\u2500, tau: 1\u23ac\n",
        "\u23a9            s\u22c5\u03c4 + 1        \u23ad"
       ]
      }
     ],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Launch the configuration GUI\n",
      "m.configure_traits()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 48,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "A complete TransFuncTrait class:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def symbols_set(expr):\n",
      "    if expr is None:\n",
      "        return set()\n",
      "    params = expr.atoms(sympy.Symbol)\n",
      "    # Remove the Laplace variable\n",
      "    params.remove(symbols('s'))\n",
      "    return params\n",
      "\n",
      "\n",
      "class TransFuncTrait(HasTraits):\n",
      "    expr = Instance(sympy.Expr)\n",
      "    f_min = Float(10)\n",
      "    f_max = Float(1e3)\n",
      "    \n",
      "    def _expr_changed(self, old_expr, new_expr):\n",
      "        print('new expression: {:s}'.format(str(new_expr)))\n",
      "        old_sym = symbols_set(old_expr)\n",
      "        new_sym = symbols_set(new_expr)\n",
      "        print('params: {} -> {}'.format(str(old_sym), str(new_sym)))\n",
      "        # Remove the old traits:\n",
      "        for p in old_sym.difference(new_sym):\n",
      "            self.remove_trait(p.name)\n",
      "        for p in new_sym.difference(old_sym):\n",
      "            self.add_trait(p.name, Float(1.))\n",
      "        \n",
      "        # store an evaluable tranfser function:\n",
      "        self._params = list(new_sym)\n",
      "        self._eval_tf = sympy.lambdify(self._params+[s], new_expr)\n",
      "    \n",
      "    def eval_tf(self, s_eval):\n",
      "        '''evaluate the tranfser function at s=`s_eval`'''\n",
      "        params_val = [self.get(p.name)[p.name] for p in self._params]\n",
      "        args = params_val + [s_eval]\n",
      "        #print(args)\n",
      "        return self._eval_tf(*args)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 207
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tf = TransFuncTrait()\n",
      "tf.expr = H_fixed\n",
      "tf.expr = H"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "new expression: 2/(0.001*s + 1)\n",
        "params: set([]) -> set([])\n",
        "new expression: K/(s*tau + 1)\n",
        "params: set([]) -> set([K, tau])\n"
       ]
      }
     ],
     "prompt_number": 208
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Visit the traits:\n",
      "tf.K = 10\n",
      "tf.tau = 4.5\n",
      "#tf.configure_traits()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 209
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Evaluate the transfer function\n",
      "print(tf._eval_tf(10,4.5,2))\n",
      "tf.eval_tf(2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.0\n"
       ]
      },
      {
       "latex": [
        "$$1.0$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAABkAAAAPBAMAAADjSHnWAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAzRAiu5mrdu/dZoky\nVEQKohj3AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAeklEQVQIHWNgYBBiAAGm6NIDQMrkE4jDwHiB\naTIDg0oYhLeCgWEjUJAdwnNjYLivAOd9Y2B4fwHGY/oK5BXAeMy/GRjOb8DOYwLKIVQyAPXdh5vC\nMIWBYT3ChgoGhkCo7foCDJwXmHoYGFgzfmYxcCYwMG0vOwCUQwYADrco9JpgipkAAAAASUVORK5C\nYII=\n",
       "prompt_number": 211,
       "text": [
        "1.0"
       ]
      }
     ],
     "prompt_number": 211
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "*To be continued:*\n",
      "\n",
      "* automatic **numerical** evaluation of the transfer function based on fmin, fmax and npts traits\n",
      "* **plot** the Bode diagram\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}