Old Notebook

gistfile1.json

{
    "worksheets": [
        {
            "cells": [
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "Julia and its Ecosystem\n",
                        "-------------------------------------\n",
                        "-------------------------------------\n",
                        "\n",
                        "[Julia](http://julialang.org/) is a new language specifically designed for numerical computing.\n",
                        "\n",
                        "Syntax wise, if you are familiar with MATLAB, Python and R, you'll feel right at home.  Under the hood though,\n",
                        "Julia tries to have the best of both worlds by combining expressive syntax and all the convenience that dynamic\n",
                        "languages offer, with a very fast JIT based on LLVM.\n",
                        "\n",
                        "Julia right now cannot compete with R and Python as far as breadth and quality of packages in the scientific\n",
                        "ecosystem, although that's quickly changing.  In the meanwhile though,  calling C code from Julia is a breeze \n",
                        "using [ccall](http://julia.readthedocs.org/en/latest/manual/calling-c-and-fortran-code/), and there is an excellent\n",
                        "interace between Python and Julia ([PyCall](https://github.com/stevengj/PyCall.jl))\n",
                        "\n",
                        "To show just how easy it is to mix and match some of the best packages that Python has to offer in Julia, let's\n",
                        "examine a simple example.  Note - this entire post is written in iPython, using iJulia, so you can just download\n",
                        "it and modify it at all."
                    ]
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "Ordinary Differential Equations\n",
                        "----------------------------\n",
                        "Let's define a very simple differential equation describing the concentration of a protein $y$ \n",
                        "which is synthesized with a zeroth order rate $k_1$ and is degraded with a 1st order rate $k_2$:\n",
                        "\n",
                        "$$\\frac{dy}{dt} = k_1 - k_2 * y $$\n",
                        "\n",
                        "For very simple equations like this one, we can just use SymPy to solve them analytically.\n",
                        "SymPy is available in Julia through the excellent [Julia SymPy](https://github.com/jverzani/SymPy.jl)\n",
                        "package. "
                    ]
                },
                {
                    "cell_type": "code",
                    "metadata": {},
                    "outputs": [
                        {
                            "latex": [
                                "$$k_{1} - k_{2} y{\\left (t \\right )}$$"
                            ],
                            "prompt_number": 1,
                            "text": [
                                "k1 - k2*y(t)"
                            ],
                            "metadata": {}
                        }
                    ],
                    "input": [
                        "using SymPy\n",
                        "k1 = Sym(\"k1\")\n",
                        "k2 = Sym(\"k2\")\n",
                        "t = Sym(\"t\")\n",
                        "y = SymFunction(\"y\")\n",
                        "# These are all symbolic variables now\n",
                        "\n",
                        "standard_diff_eqn = k1 - k2*y(t)\n",
                        "# Note that SymPy pretty printing using LaTex is enabled by default"
                    ],
                    "language": "python",
                    "prompt_number": 1
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "SymPy requires us to set the equation to zero before we can solve it, so we have to re-arrange it a bit."
                    ]
                },
                {
                    "cell_type": "code",
                    "metadata": {},
                    "outputs": [
                        {
                            "latex": [
                                "$$y{\\left (t \\right )} = \\frac{1}{k_{2}} \\left(k_{1} + e^{k_{2} \\left(C_{1} - t\\right)}\\right)$$"
                            ],
                            "prompt_number": 2,
                            "text": [
                                "             k2*(C1 - t)\n",
                                "       k1 + e           \n",
                                "y(t) = -----------------\n",
                                "               k2       "
                            ],
                            "metadata": {}
                        }
                    ],
                    "input": [
                        "diff_eqn = k1 - k2*y(t) - diff(y(t), t)\n",
                        "sol = dsolve(diff_eqn, y(t))"
                    ],
                    "language": "python",
                    "prompt_number": 2
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "This solution has an constant $C_1$ which is determined by initial conditions.\n",
                        "\n",
                        "Let's set the initial condition $y(0) = 0$ and find a fully determined solution."
                    ]
                },
                {
                    "cell_type": "code",
                    "metadata": {},
                    "outputs": [
                        {
                            "latex": [
                                "$$\\frac{1}{k_{2}} \\left(k_{1} + e^{k_{2} \\left(C_{1} - t\\right)}\\right)$$"
                            ],
                            "prompt_number": 3,
                            "text": [
                                "      k2*(C1 - t)\n",
                                "k1 + e           \n",
                                "-----------------\n",
                                "        k2       "
                            ],
                            "metadata": {}
                        }
                    ],
                    "input": [
                        "rhs = sol[:rhs] # Get the right hand side of the solution."
                    ],
                    "language": "python",
                    "prompt_number": 3
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "The syntax to access the methods of Python objects through PyCall is unfortunately a\n",
                        "little ugly, and we have to type sol[:rhs] rather than sol.rhs.  This is currently unavoidable\n",
                        "   unfortunately."
                    ]
                },
                {
                    "cell_type": "code",
                    "metadata": {},
                    "outputs": [
                        {
                            "latex": [
                                "$$\\frac{1}{k_{2}} \\left(k_{1} + e^{k_{2} \\left(- t + \\frac{1}{k_{2}} \\log{\\left (- k_{1} \\right )}\\right)}\\right)$$"
                            ],
                            "prompt_number": 4,
                            "text": [
                                "         /     log(-k1)\\\n",
                                "      k2*|-t + --------|\n",
                                "         \\        k2   /\n",
                                "k1 + e                  \n",
                                "------------------------\n",
                                "           k2           "
                            ],
                            "metadata": {}
                        }
                    ],
                    "input": [
                        "init = solve(subs(rhs, t, 0), :C1)\n",
                        "def_sol = subs(rhs, Sym(\"C1\"), init[1])"
                    ],
                    "language": "python",
                    "prompt_number": 4
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "We can simplify this a bit:"
                    ]
                },
                {
                    "cell_type": "code",
                    "metadata": {},
                    "outputs": [
                        {
                            "latex": [
                                "$$\\frac{k_{1}}{k_{2}} - \\frac{k_{1}}{k_{2}} e^{- k_{2} t}$$"
                            ],
                            "prompt_number": 5,
                            "text": [
                                "         -k2*t\n",
                                "k1   k1*e     \n",
                                "-- - ---------\n",
                                "k2       k2   "
                            ],
                            "metadata": {}
                        }
                    ],
                    "input": [
                        "def_sol = simplify(def_sol)"
                    ],
                    "language": "python",
                    "prompt_number": 5
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "This is a symbolic expression.  To evaluate it, we need to turn it into a Julia function.\n",
                        "\n",
                        "SymPy has very advanced functionality to do this through the lambdify module, but it's overkill\n",
                        "for what we want to do.  As a quick solution, we can simply replace the parameters with fixed values, and write\n",
                        "a function that replaces the $t$ symbol with a value, and then convert the output of that substitution\n",
                        "into a julian float."
                    ]
                },
                {
                    "cell_type": "code",
                    "metadata": {},
                    "outputs": [],
                    "input": [
                        "k1 = 0.01  # Arbitrary Parameter Values\n",
                        "k2 = 0.025\n",
                        "\n",
                        "bound_sol = subs(subs(def_sol, Sym(\"k1\"), 0.01), Sym(\"k2\"), 0.025)\n",
                        "\n",
                        "function eval_analytical_solution(v)\n",
                        "    return float(subs(bound_sol, t, v))\n",
                        "end;"
                    ],
                    "language": "python",
                    "prompt_number": 6
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "Let's evaluate the function now between 0 and 100.  We can take advantage of the @vectorize_1arg macro\n",
                        "to do it quickly."
                    ]
                },
                {
                    "cell_type": "code",
                    "metadata": {},
                    "outputs": [],
                    "input": [
                        "vect_eval_analytical_solution = @vectorize_1arg FloatingPoint eval_analytical_solution\n",
                        "\n",
                        "t_steps = [0:1.0:100]\n",
                        "analytical_result = vect_eval_analytical_solution(t_steps);"
                    ],
                    "language": "python",
                    "prompt_number": 7
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "Let's compare the analytical solution to a numerical integration.  We can use the [Julia wrappers](https://github.com/JuliaLang/Sundials.jl) around\n",
                        "[Sundials](http://computation.llnl.gov/casc/sundials/main.html), one of the most advanced ODE libraries.\n",
                        "\n",
                        "Currently, in the Julia world, there is a big push to improve [ODE.jil](https://github.com/JuliaLang/ODE.jl)\n",
                        "but the package is currently in a state of flux, and I'd suggest just using Sundials.jil for now."
                    ]
                },
                {
                    "cell_type": "code",
                    "metadata": {},
                    "outputs": [],
                    "input": [
                        "using Sundials\n",
                        "function f(t, y, ydot)\n",
                        "    ydot[1] = 0.01 - 0.025*y[1]\n",
                        "end\n",
                        "\n",
                        "y0 = [0.0] # Initial Conditions\n",
                        "num_result = Sundials.cvode(f, [0.0], t_steps);"
                    ],
                    "language": "python",
                    "prompt_number": 8
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "Julia has two excellent and highly developed plotting libraries ([Winston](https://github.com/nolta/Winston.jl),\n",
                        "[Gadfly](http://dcjones.github.io/Gadfly.jl/)) as well as a very advanced wrapper around [matplotlib](https://github.com/stevengj/PyPlot.jl).\n",
                        "\n",
                        "I'll just use Winston for this simple plot - the syntax is very straight forward."
                    ]
                },
                {
                    "cell_type": "code",
                    "metadata": {},
                    "outputs": [
                        {
                            "output_type": "display_data",
                            "png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEsCAIAAADfNCTgAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVxN+f8H8Ne9t5VIQ4oWSWJQCCO7ypJUosVSyJ41g7EMY50x9j1EVGRLlkSLpexLihRFUmjfS3vde8/vjzvffqah0HJu9X4+5o/bueee8+revOYs93wOh2EYEEII+VFctgMQQkj9RjVKCCHVQjVKCCHVQjVKCCHVQjVKCCHVQjVKCCHVQjVKCCHVQjVKCCHVQjVKCCHVQjVKCCHVQjVKCCHVQjVKCCHVQjVKCCHVIsF2gAbozp07ZWVlbKcgpC40bdq0X79+bKdgGW2N1jAPDw8nJ6fK51m0aFHlM1y9ejUlJaWSGRiGOXbsWOULcXFxqXyGmzdvOjs7V3MhV69eTU1NrWQGhmGOHz9e+UKqfEPevHlz7969yuepMqqPj0/lUfl8/uLFi6u5ltevXz948KCaC/nrr7+ioqIqmUEgELi6ulZzLVFRUVVGrfKj+eOPP548eVL5PA0ebY3WPFVV1WHDhlUyg7y8fOUzvHv3bsCAAe3bt//aDAzDBAYGVr6QmzdvVj5DcnIyn8+v5kJiYmIGDBigoaHxtRmEQmFQUFA13xB5efmPHz9WM+rbt28HDhzYrl27r81QWlpaZZIq19KsWbOkpKRqLsTFxaVPnz76+vpfm4HP59+7d6+aa5GTk0tJSanmR3P06NFKnm0kqEZrmLy8fKtWrSqfR1FRsfIZtLS05OTkKpmBw+Ho6upWvpDu3btXPkPz5s1btmxZzYVoaWk1bdq0khm+JWqVb0jLli2rHF+8+lG5XG7r1q2ruZaWLVtyuVXs5FW5kNatW8vKylYyA5fL1dHRqeZaWrVqJSFRRQNU+dEoKipWHrVRYEiNCgkJ8fT0rHyeFStW1E2Yyj19+vT8+fNsp2AYsXlD+Hz+77//znYKhmEYJyenDx8+sJ2CYb7hozlw4MDHjx/rJozYomOjhBBSLbRTz4LJkyezHQEAtLS02rZty3YKQGzeEB6PZ2try3YKABg1alSVe9N1Q0w+GjHHYeiWdjUqNDR027Zt7du319bWnj59OttxCKmWkhJIS3/5qX379iUlJUVHR+/du1dNTa1uc4kXqtEaFhoaGhsba21tzXYQQn5ESQmki3NjfKKeBmQpvbxl/3LZlT6beiweCmtrcDj/nd/Jycnc3LyR1yjt1BNCgPj4E055Zlnulp7jNQsi/PnDkvGLFAxLITXs0aabj4b18PLCuXNfbFJCNUpII5UUmVNw42Gif8T6+8P65AftwDJ5/J4L+SDoiWYohsxybLuIccNw88b54T0necPCgt3M4onO1BPSiKTGFQoPOS/sfveC0jyDrmmDF+uZ+i+4k99rB5YByIV8+ZytkDEcN3ZgmR08muPTcNx4fvw5e8HFGm2NEtLQffr0wT/qzqkE2ad37ZP/7o8ONzH4AAZX/qIMtNJCDAPORqydhNMesJt91y6Yod36L6AaJaRhehGYqXPXaa67vkJy5LkSiw/4hYcxfEjcRGUXd37uEOa2QI4QXA/YaSL2wpgTHM6mWs1cT1GNEtKA5ORkXH/mdIDReeEx8ZOzFqwi0QUYIXqS//3/3nPQAoAmYoM4hupzz9Zw2oaCapSQeu9jVEFbv2P7TrR4ESUVVtolHLqAEYBIdKn+wjURGyQxXH2TA74+VEojRzVKSP0kEDAvwl13ZPa9s21Y0ok2GPQcPWtkwRwwqkjIRMudEivc5Redm3hZfd5V/PxzjSy8QaIaJaQ+KSsD7+G9kJNRbpfkBVm5RzBbEkPKIJkC5eostjk+fULzUfBr2ay0b8esCTbCqLZGg8bvdJCSAn6rqfANFdUoIfVDxLWP2p5/Trw8XvCp4A7Gl385qQySP7ZAHgRcCOfikK/02Jume26r2E6Y3VG6c3vweAAG1Vjwho9qlBAxJhCURr8/vfa16i1382x3Raz5CPVqLlIX4VmcluOa3xwzNDdHU2/cnJHbNNSkpXdOrZHAjRLVKCFiKTLyruu7Jyej76Z2ugpTHowF4P1whyohNR9yrhKzLrVfsn9elGC4sWKXqeXfAP3K2CPkW1GNirH371FYiI4dIfmDe22k3hFm58LFZd4hHaW4x7uwJB9moukC8H5gaVqI0eR90FQu3G73Iqr7hD6WbtZSUkDvGo1MqEZrwYcPH549e9a8eXMtLa0fXISXFxYvRmIiAMjIYOVKrFz51QHL/uf58+erVq1iGKa4uHjKlCnu7u5Hjhzp3LkzgKSkJHNz85CQkKdPn65evVpCQkIgEGzZsqVnz57v37/X09PT19fncDgKCgp///23mpra+/fve/bs2atXLwCysrI+Pj4/+IuQb/bJOyjN5coYPwdtQYfLGFH+Zc/vpYyUNLR2kXCI6DphycxPrcf0k1RV4nDM+tRsXOD169eFhYWV33uxkaAarXlpaWkxMTHKyso/WKOenpgwAeX38ykpwfr1ePsWHh6VvCg3N9fW1tbX11dDQ4NhmPv370+cOPHs2bPr168HcO7cufHjx+fk5NjZ2d24cUNdXT0uLs7Y2Dg0NBSArq6ur68vgOvXr1taWj5+/BhA9+7db968+SP5yffIeZ8j3L7TzbPJngw7BksToBqJTuXP8iAQbYdyIRRWOgJGC+QYcwKy5TUu2l4I1Z8/aNJhVHVLqGqKj4/Pzs7Oysqq1bXUD+zew6Th+ZZ7MVWhXTuGy2WAf/3H4TBv3lTyIg8Pj8WLF38+JT09vWvXrqLHffv2jY+PP3nypKOjY/kMc+bM8fLyiouLGzJkSPnEgQMHRkVFVZhIap5QGH7wXnq3ob05IWr4WOHTBpg2SAIYDoRSKBFN4UJQYZ6myAcYBxza2MHt6R/ewui3JSV1/XvQvZgYhqGtUTGTnIwPH74wnWHw8CG0tb/2usTERFVV1c+ntGrVSkND49mzZ/Ly8tLS0qqqqomJiZ8PryuaItpzL6eiopKUlKSpqRkeHm5sbAygQ4cOTk5O1fy1yP/LzLzz62XmvJdpsRcP3p/Q/PMnZVCsgOxktMlBi7ZISkLbMkhywDDgfL412g+PPnDbB/ZZETlkrsWK8ZyfFETTper0NyH/oBoVM3z+V58qK6vkdaqqqk+fPq0wcdKkSWfOnJGXl584ceJ/54mPjx85cmSFlyQmJqqoqADQ1dX19/f/zvSkMhlBEVFLXbyfq+/EUmnYlfzvDLkkysogqYJELoTxUCuCrBJSU6Ek+kIog/8fUkkN8cMkbqvryK+fl5Zq3FtJ1b3Tl1dF6hbbm8MNTXV36gUCRlGR4XAq7uMBzLNnlbwuJyenS5cu79+/ZxhGKBTev3+fYZj8/PzOnTt37949IyODYZjs7GxtbW3RnXvfvXunra2dl5f3+f57QEBAnz59BAIB7dTXpNzcx5tvhSsP74HnPPDLP09JlMohD2B6IUQZyQCjho8/IbPCxy56yQpsWap5MWbnZSYhge3f519op56hnXqxw+VixQosW1ZxoqEhelZ2xbS8vLyHh8ecOXMYhiktLZ08efKAAQOaNm2qp6eXm5vbsmVLAC1atDh58uT06dN5PJ5QKDx79qycnFxGRkZ4ePioUaNEZ+q9vLy4XC6AFy9eDBv2z4hqfn5+kvSlqx+Ql3d34XnOKY/RfO8C+JfvlasgUbSxaQLfQBiGopcenqVDMR7/uqORElI5YALUZz03XjV1zSQ07vsdiTW2e7yhqYFTTEIhs3UrIyv7/yeX7OyYzMwaCkjqQtnbuHvD1h/hzOZC8PkWaGdEiU4NzcVBCZQBjD4eVdj8bI7crng5v+WZT45r3gZnsf2rVIG2RhnaGhVHHA6WL8esWXj+HIWF0NFBu3ZsZyLfKjMg5NPSDate2V3AagF4oiObHDC9EBqC3gLwJuLMGUw8icm9EfIY+o/x/6PPGSAoldfWZfSlbtP6yI2ZwOGgGXu/CPl2VKPiSkEBhoZshyDf4YOzf5O//xjx4WgUvMpPH6khXghuIlQ64Y0Csm9guGhiPNTKC5QHwUpsea089OSKaOnJPbgtV7L2O5AfQjVKSPUwzKfjXvFbTxu+dS7BzfKBl/riSRh6JKHtWmw8gtmnYGsGHy6Eb9Gx/KUaeO8ofVjRtK/tihHo0wcYwNLvQKqFapSQHyUQwMkpYHfkpPebJTAkDa1Fk/Xx+DH009B6F5b8it0bsdYMPolQ8fnfBfIAbHGqm3LGxHXa7ayWoVUrln4BUjOoRgn5fgIB3Nw8t8TujxkVgtnFkAEgiyINvI/CzzqI6I4XzpizA8sG4EEQDC7jn9u7cyF04iwo7DPE8VA3Xg+d2r5ek9QNqlHxVVqKtDT8+9IkwrayMpw9u3VldkIS9zA2iG4S1xpp+ZBjwFmP9Vux4ihmzcJRJaTGoX0c2otep4TUo00c863sJ25YDg0NNn8FUtPof4ZiqrQUlpYYMABxcd/xKltbWzMzs6rnAwC8fPnS3t6+wsS0tLTLly9/7dlvWUKDJRQiMHC/6tZ1U+JWJi06gAV8SMgjF0BXvPKEDQeMPdymwp0L4VHMSoWS6HUDcd+j9ZJbB9+YZbpNdDemDm14aGtUHJWUwMoK165BTg5DhyIoCJqaVb+qoKAgPj5eVlY2MzNT9H37HyCqUQsLi27durm5uf3YQhqgN2/Cxm08GdlrF9aIJojGXvoLq5+ijzumSqF0LTauxBZH7C1/0QpsVe7QdKZTT7mhf1c5ziGpv6hGa56Xl1doaKi2tvb06dN/4OVFRRgzBrdu4fhx6OvD0BD9+yMwEF2quleut7f3uHHj5OXlz58/7+DgEBISsnjxYn19/ZiYmClTpowbNy4rK8vW1lZKSqqgoMDd3V30quXLlw8ePNjU1DQzM9PCwkJfX//+/fv29vYmJia+vr5ubm5Lly6Niori8XgrVqzo0qXLf5fQsDGv35Q6/rbyhtFexqP88nYZFN+C0SScdsReN9inQskfxtf/N0KoBPgHOfNL9PotcDdF167sZa9d+/btS0pKio6ONjc3ZzsL29j+/n9DU82rmIqLmdGjGQ6HOXz4nynPnzMtWzLq6sy7d1W8dsyYMSkpKbm5uSNHjmQY5unTp4MHD2YYpqioaMCAAQzDlJWV8fl8hmEuXLiwZs2aiIiIqVOnvnnzxtLSkmGYvXv3Ojs7iyYyDCN64O3t7eDgUL6KLy7hh39ZcZeSIly4aKHU4c8HshuKoL+xEmD6IPgZerZGKg98GRSVX4C0k7P0vKkbExXFdvo6QlcxMXQVk7iJjMTt29DQgI3NP1O6d8fIkTh9Gj4+cHT86gszMzODg4OnTp0KICIiIiEhAUDHjh0ByMjIlJSUACgsLFyxYkV+fn5eXp6Cwj9Dq2lra+fk5GRkZJw7d87Pz+/jx4+fL/bly5cDBvz/lxm/uIQGqKyMOef524ys96WDL8BSNG02joSg9wMMmA+nbVi+HNts4JmBVkJwBeC1w4dh3MDpk4r7L58KHR1245M6RqeYxEvPnrhxAxkZGDAAKSkQCjFzJk6fxs6dlXUoAE9Pz/Xr1/v7+/v7++/atevMmTOfP8swDIBTp0517dpVNDqJaIrIlClTli5d2rFjx+bNm0tJSfE/G6yva9euDx8+LP/xa0toUJ4/v9zOcfHkzJ2lC0UdOhIBXAjj0P4KzHUQMRFnwqELIAZaQnC744U+50ng/AsuMUP7n5xLHdoI0dao2OnXD5cuwcwMo0ejVy8cP441a7BkSRWvOnPmzMmTJ0WPhw8fPmLECAMDgwrzDB061NHRMT4+Pjc39/Pp1tbWixYtEt1wqX379qmpqVZWVpaWlgDMzc2DgoJMTU2lpKSWLl36tSU0EHFxJQuWbvTtvRkHRRMUkP0JzZdhhw08Z+HoNLga4dYz6HnADkAXREpzSm/YubdcOx9aVX1CpAFj+aBCg1MDIzwxDMMwAQGMjAwDMKtWVX9hVSgoKOjbt2+tr0ac5eUxa9ceaLpcEWnlh0G18eYj1LTxRhaFNzDMHN7lT8khbzTXt3jC1OJXMWxHZxkdG2Xo2KjYGjECly7h/n38+Wftrig0NHTRokVLly6t3dWIs0ePLo9xvZ7e8zAcROfiJ+OkItJ3YclS7LyOEb8g2AS+orHoAVjg8pFezs2P7JDWc2MzNhEbVKPiy9gYxsa1vpZevXo9ePCg1lcjnrKzU20WPrxZOB7nRC1pAt87GHIDwwNhKABvLxzfomP5xfJTcMJKPXjk6alSA/xYzU3EC51iIo2SUIgjR150mdjt5h5rnC+DZBMUAhiCO9cwOg/NDBEoukIpDD0A2MFjndzO45cUzGL3Sg2o8Vu+k/qNapQ0Pu/fp5rOmDWHa5RyKgOtBOA1Q94z6PVA2ApsDcYvjtibAuWNWAugB8KmS55035W5PsWBZ2EGHo/t9ETs0E49aVxKtu559cfZyWXHI9EFQC+E9kLoEcyeCvcLsDSDz3Js40IIgAeBJffSKevLEju3QkWF7eBEfNHWKGk0oqJKR4+1XtmhV9njSHRRQioXwhho2cNtNo48Qd+xuPQeGgCE4JrAN6qj+akHGhJnPahDSeWoRr8gLCysR48e7du3NzExycrKqvCskZGRqqqqiorKiBEjKlzzQ8QUw+Ds2aM6+0x955WPnTwXh45gdh6aGcN/EO7JoDgcuoVoYoCg7VKrL7rmdnx1WUK/N7vBSb1ANfoFM2bM2LRpU1xcnI6OzoYNGyo86+rqmpCQkJCQ0KNHj99//52VhOQ7fPxY0m+o08T7cwQHb2A4B8weLJZH7gasK4bMHiz+hOaTcbIYMlqIGYmAq4uuL4t3lLafCLqnNPk2VKMVxcbGJicni0btnD179vnz5yvMoK6ujv9dttC6dWsWIpJv5+kZ0X9OhyenF2I/A05LZDLgPIb+NYxuhryF2L8Hi0Uz6iDiYZeZ/te5Tfb+DfpYyfegU0wVxcfHq6mpiR6rq6unpqby+XwJiX+9USYmJo8ePVJTU7t3795/lxAREdGsWTNZWVllZWUAysrK8vLydZCc/EtZWZH9XLfTUutwIh2KANrhwyWMHYtLZzGBAWcErnvBKhaaikg/wbXv5dBH8UAQOBy2c4u7rKys9PR0AElJSSUlJTExMWwnYh/VaGWYr4y+4evry+fzV65c6eDgUGEQEADBwcGJiYkKCgpdu3YFMHDgQKrRunb7duH0BeZxe2/BCIAx/LOh8AR9Z+PIZViYwPccxgOQQqkx/P/Uu6hzaic6d2Y7dP2QmJgYEhICICwsLD8/n04PAHRN/X+8e/euTZs2oscxMTHlj/8rIiKiXbt2FSbW1DX15MedPn2UO/tnRAIMD3yAUUFCMPr0xlPRlfKiiVIo8ZafzGzbxpSUsJ24HqNr6hmGoWOjFWlqaiorK1+5cgWAs7OzlZWVaLq/v39mZmZmZmZ0dDQAPp/v6urauzedyRUn6eklg4YdnRQ0R3goCj8DcMHM4biRCBULXF6CXRww0dAWgLcc26K0x5g//QO//QYpKbZzk3qO7R4XR6GhoTo6OqqqqiNHjszIyBBN7NChw7179z58+KCnp9emTRtVVVUbG5uUlJQKr6WtUdbcuZPcaUgXvOJAyAO/D4IBRh45gTAYgYDywZn0ELqa8xdz5AjD57OduCGgrVGGRnj6Ij09vfDw8AoTyw+lh4aG1nkiUpXNmwNX35qPw6/RGQAPguXYdgq2l2Fhgct98UQ011Dcvqr/Z9N9f6MPXRdPagzVKKnn0tOxevVml9arcQtAb4QoI+UqTG1xygN2aWj9EP1vYHgrZDhzHIxnqzc5dINOx5OaRTVK6rO0tMQeo/9InueKaaIJz9HTBTM5YHxgZgcPBhwOmJ8Rdbr96u6nV0Bfn928pEGiU0yk3nJ3T+xsZJB8WtShK7DVBL4C8GbCZRhuyiG/FFJ8SOzmLHk1fmP3h4eoQ0ktoRol9dPu3SfsAwdne79Fx5+QBWAnlk7A2dG4JgDPEXvzITcAD05zbB09B+DsWSgrs52YNFhUo6S+ycjA6NGuSyKmwTUWmgC2YfkAPOBDYjqOiyoVwEDc9+u1ZkLIMvzvK2uE1BI6NkrqlZycQlObqU/mXcS4ZshricxYaM7FIWfM4YC5j4EnMVkaJScwxWQEX+6yL2Rl2U5MGj7aGiX1h7d3Yfuupk/+8IKVENxSSK3Fxo54WwbJ2TiSCiUAiki/IDXRxttWLuACdSipG1SjpJ54+DDU6u9OOU+CYKCPxy2RWQTZuTj0Ozar4yMfEm/R0Qw+CSp9R/svhLk523FJI0I1SsSeQIB580IHOg7n+yVAFUAsNPfCsRUyRE2ajDYyKF6MPV4mx6XeRMDAgO3EpHGhGiXijWGwadPBQ5yBzL1sKEzAWQnw09B6CXZtx288CIohw4PAmzt299IEKa/TaNqU7cSk0aFTTDXPy8srNDRUW1t7+vTpbGep51JSYG5+8GmfBTjAgAMgAjq78euv2J2G1jPhIgBPH4+3clcNvvQr7cjXsX379iUlJUVHR5s3+neew3xlSE3yY0JDQ2NjY62trdkOUv8VFAjNLRYEjjsMB0mUjcUl0SChuggfjWt/YxWAfnjkr/Nb88Pb0L8/23EbKScnJ3Nz8/KRzhsn2qknYunpU2HnLjMDJx3CXAacUkiVQHo6jgMIh+7fWMUBsw4b/NVmNb97lTqUsItqlIifzMzcURNGJhxzxbSueKWNaACXYfEJzY1wCwAHzH4sXD8/vXnYXbRowXZc0thRjRIx4+mZq9VrZObpmxgGIBJdLHGhE94A8ILVLRi1ROYtGM2fVogDB/DTT2zHJYRqlIiVZ8/iJy4fnHPlCfoOxW1JlDHgbMFKG3gqIBvAT8i6ITfW4MxsHDvGdlZC/kE1SsTGoUMJQ2wNhTfDoQvgFbquxl+iJt2EP7KhMBS3wzg9e7o5YsIEGjOUiA+qUSIeXF1fz9tnmH8lBlq9EAogHYoHsGAZdoieH4rbV7utVAu5BEtLVoMSUhHVKGEbw8DFJcrxsAGC3qIjgJbINIEvgAy02oKVABbgwLUWtk0vnoSeHstpCfkPqlHCtjVrImftNszzLoKsHPIBXMcIAKImZcBZiP37xj9oEh2Gjh1ZjkrIl1CNElb5+r7ads0QgSlQLobMPBwUNakvTO5hEBfC/Vi49+fDHNfjUFRkOyshX0Y1StizYMGr0csN+dfToQigBNJ7sHgeDjZDHoACND2O6QvmM5z792jIOyLO6Jp6wpLbt1863TFCYBpa90ZIUxTcwZBSSO3BYmWklEDaA3bWJgU44MZ2UEKqQFujhA2bNr0c/qsRbqWhNYAQ9C6F1BDcAVAKqRQoe3InWC9rh3Pn2A5KSNWoRkmdO3UqZu0JI35AGloPwAMeBAAeoV8RZCVR1gx5l2ExZlE7bN8OOTm2sxJSNarRmpednZ2QkJCRkcF2ELF06lTMwr0GCEpDay6ET9HHEhdETRqMX6RR4ss1HTVfE1u3sh2UVCEtLS0hISEvL4/tIOyjY6M179mzZ6Wlperq6jQOY0Vubh+nrTPAg0SoqCAxCW0F4HljzHDc8IexHPKvYfTAXePg6Mh2UFK1oKCg9PT0169fsx2EfTTeaA2j8Ua/KjExtfuIoZkX3qEDHxK9EJqMNkloywEjBLcZ8nxhMnCiGk6cgAT9373eoPFGQTv1pI5cupTWoZ9h5vnX6MyHhBZiQtGrNdKUkSLqUD+MGug6A6dPU4eSeof+ZEnty8wsnOJgWXIhEl0UkK2KhJfopoWYMPSQQqkkyk5gygDrtrC3ZzsoIT+CtkZJLYuMLDCbMDr/3H0MbIskBpx0KHbDyxhoNUWBENyzmGBxYBjOnmU7KCE/iGqU1KaMDGaoge2jBbcxVAmpSWirjo8lkE6HYhMUFqCpC2fWuLEM5s8Hl/4USX1Ff7uk1giFWL9+afpKb4xpg+QMtOqBsHDotsOHNLQuRJO9nMVTZkjCzY3toIRUCx0bJbUjPx8jRvz2aNxu/ApAFQlqiA9Fr254GYYeAPbCceGkTBz1YDsoIdVFW6Okdqxf7/1YaSeWAlBB4jPoMeD0QuhLdANgA895Pwdh9262UxJSA6hGSS3w8blyMMGGOceDwAw+iVBRRcIz6CWjDQAbeJ6y8ZYIC6Gx70jDQDv1pKbt3fts8QlrPCqDpAT479BhLC5dwtgWyImH2i8IPtV0tsTWMEhJsR2UkJpBW6OkRuXnJ67cPwFnJcAHIIuit+j4Bp0G4EEOWvyMKO9uqyVu34SGBttBCakxVKOk5rx6lTjAxqDY7y06KiNlBbbmoZmoSR9gQBdEBkkMV758GL17sx2UkJpENUpqSEEBYzTMMXyG6LZ0sdD0gpWoScsgKYPi/ZxFSrtWoEMHtoMSUsOoRknNYI4dd0z9/QIsjeEviTIAMdDygpUa4mVQfIkzzjB0OxYuZDsmITWPTjHVPC8vr9DQUG1t7enTp7OdpU4wDBYuXOakuR9LAKRA+RDmzsPBUkjFQAvAKY6dsb0yevZkOyipSfv27UtKSoqOjqYBIWmgvBrWGAfKO3Xqjd2mIbiTCiXRhB4IW4ADDjjMh4QdPNymBPIOO9Ft6RokGigPtFNPqqug4PVWb9Fo9pvxuzRKAIShxzps4ENiEk67jfXmuR+nDiUNGNUoqYaYmCItHaOIPclow4Djh1GumCaFUgCi8e3dZebw1v/BdkpCahfVKKmGOXM2pDokoa3op3sYdAhz/8JqAPLIvaC1QuLxfejqshqRkFpHNUp+1O3bKwNHbGWWayHmFGxF37e/h0ErsUUeuQEY2XenDbp3ZzslIbWOapT8kEOH/A23bcUKADHQuo4RJzFZ1KQC8ObDqa+tFhr9CVzSSFCNku+Xnp7+6+bfmG0yKOaAAeCOqTcw3AaeAGzguWGnHDxoBDzSWND3Rsl3EgozNh8ZVnLtJbotwS5tRM/FIQac45gOwIZz/pTSUom50WynJKTu0NYo+WJqKZMAAB6jSURBVB5CIUxNN+yRD4cugD1YLAH+IcwVbZPyIFjS7KjE+TP09SbSqFCNku9x7twtv9JjmCGPXABCcGfjSCmkOuGNBPgeUtP7vjuNgQPZTklInaIaJd8sJ+fmlhAz+BRD5hpGz4QLACG4i7AvBlonOVMm2MugVSu2UxJS16hGvyAsLKxHjx7t27c3MTHJysr6/Km0tDRTU1NlZWU1NTV7e/vi4mK2Qta11FTo6CwOn14EWQac8Tj3G7aLmhRAO3wYMzQX27ezm5EQVlCNfsGMGTM2bdoUFxeno6OzYcOGz58SCoUODg4JCQkxMTGZmZk7d+5kK2RdW7t2d6LNK3QV7c4nQsUQgQYIkgBfE7GBQzbIBl5D8+ZspySEBVSjFcXGxiYnJ5uZmQGYPXv2+fPnP39WWVnZ1NRUQkJCWlp60KBB8fHxLMWsWxkZuz0UlzA71fExBlr98RBAIlRscUodH4NgoL7Mhu2IhLCGvvBUUXx8fPlwNerq6qmpqXw+X0Ki4htVVFTk7u6+Y8eO/y4hKyvrw4cPkpKSTZo0AdCkSROpen3foefPk4ZOWl4YAeAj1FdiyzWMHo1rD9EfQH/OI/VtC2FqynZKUkdKSkqKiooA5OfnCwSC3NxcthOxj2q0Ml8bRZDP50+aNMnU1HTUqFH/ffbYsWOXLl366aefdHV1AZiYmOjW5+vKmclT/spbxIeEFEpLIXUMMxhwJuH0Q/QfhptHNmdg2TK2M5K68/z589u3bwMIDQ3Ny8vLzs6ePHky26HYxpB/e/fuXZs2bUSPY2Jiyh+XEwgEkyZNcnBw+OLLQ0JCPD09azdiHRKGR8yDE8DMxNEH6N8cuQAj+m8YbhQ2ackkJrKdkbDpwIEDHz9+ZDsFy+jYaEWamprKyspXrlwB4OzsbGVlJZru7++fmZnJMMyMGTMkJSUPHjzIasw6ceuW54C9BzEPgBesJFHmD+Pm+CR6cpLsJdlzbmjbltWIhLCPavQLXFxc1qxZo6amFh4evm7dOtHEBQsWREVFPX361M3N7fr162pqaqqqqvb29qwmrU1ZWQJLG+98I9FPOWgxAtclwLeCF4DpOD71/iw6JEoI6NjoF+np6YWHh1eYGBMTI3rANI7brgi8LtnlOp3FhF1YEgOtg5iXgxaDcbcYMtM5rkd/3s3tWfEtIqRxoholX/L0qeNy6bOYACAa2k6YzwHjhPnFkJEA/48W+7keJ8DhsJ2SELFANUr+IzqaP8ToQ9Fp0U+H4QDgABbcwZCX6LZY0knj/W36pj0h5ejYKKmIv2OPbfGxqzBdhh3KSAFwGA66CH+Jbss4O7dPCKUOJeRzVKPkX/gfkyadNvVkrBWRvgj7gmAgatJIdOmOF9sUt2PTJrYzEiJeqEbJZ8LD0a1bWEFHAOlQHIrbTVAYCEPRbZONlcM4UZFo147tlISIFzo2Sv6f8Nels/N3vUVHPTx7Br1YaBogaCDul0B6BbZumV+Gn35iOyMhYoe2Rsk/hAVFM4LsXBn7yTj5FH1m4SiAWGiewJT5cNrSZi/mz2c7IyHiiGqUAAASE18OmX+GmQDAC1a3YHQYDqImBdBKQw4PH0JBgdWIhIgpqlECMAzGjAl8rlAKKQBFkB0D75sYpoRUAHPgvG7vT9DQYDkkIeKKjo0SIDh4d+jgJdhpiQtTcMIGnkWQHY1rfEjM5Bw7pLqZM+IN2xEJEV9UozXPy8srNDRUW1t7+vTpbGf5BkLhsQ0JS7ALwEDcN8eVC7Ach4ulkOKAsW/pw7l4ATIybKckYmffvn1JSUnR0dHm5uZsZ2EZ1WjNs7Kysra2ZjvFt2EYWFvL+UkClgCWYBcPgoXY74DD+7DIEIE97+xBFw22UxJxtGjRIgBOTk5sB2Ef1WjjduXK/Ytps+DbCW8E4MVAyxF7n6DvWUww4gT6aC6W7RLBdkRCxB3VaKN23y1mFPykUBoEAwF4BgiKgdYp2HbGax+ehezhi2wHJKQeoDP1jdjFi1t9dfIhl4WfDmCBKhKCYNAB7wAUcJvxg+5h2DC2IxJSD9DWaGN19+4n6xkZQj/RT5vxOw8CS1zIhbwa4oP6rGw28BS7AQmpL2hrtJHK27h7FOMbil6nYGuAIACb8Ed/PJRFURAMOvza2M+9EvLtqEYbo/w8xuTOiodMPzP4TMLpqzAVNWkhmlhxL3bYMAXjx7OdkZB6g2q08UlNlTQanMtvCuAixq3C301QeBHjWiCHA0ZrUBusXct2RELqEzo22vjMmnUstMdLdGuJzEy03IKVAvBeo3Mu5A9gwby5g9nOR0g9Q1ujjUxGxsGr6guE+8zgE4f2o3ENwHb85gMzJ86CeTr3YGnJdkRC6hmq0cbFzyl2AbOfAWcQ7jVD3gVYipoUALeHLm7dggTtoBDyfahGGxNX196bxjRFAYDl2LYfC6VRYgcPAANx3/bPn6GoyHZEQuof2vRoNJKTc+auGi3wKYPkINy7h0GO2BsOXXdM7Y+Hvgq2coY0jBMhP4K2RhuL3PPXjUu8X6C7F6xuYpg5rjDguGBmO3zw45g0O7KThnEi5MdQjdY8gUBQVlYmEAjYDvIZoXD1YbUn6MuDQBolUig9D2vRUdE4tE8+cxtWVmxHJPUMn88vKysTCoVsB2Ef1WjNu3LlysaNGz08PNgO8j9lZRg1yiDqIP43sv0NDBeCKwSXC+Fh2V872XRnOyKpf44cObJx48a7d++yHYR9dGy05o0dO1a8xht1c7t+HZNxsjteNEXBQ/QfA+8eCHuCvs4ch5nTJMHhsB2R1D/z5s0DjTcKgGq0MbjummiBy9qIvgUjKZQaw/8h+j9CP0fsndkzFJtusB2QkPqNduobOEH2p6XPbIsgWwTZEkg3Q54/jHsjBMAtmdGCR8F063lCqolqtEF7/57X7WfDEj8A0dA2QFAiVF6h61t01EKMn8l+nhSP7YiE1Hu0U9+grVy5Jnn+Piyyg8dj6EdDezDuZqKlItKDYKA6353tfIQ0BFSjDZdQuPqC3mZm+SwcPQyHJLQV3SOEAyaAZ6K6xRGGhmxHJKQhoJ36BuvN7eSdfEcALZDDhVB0jxAlpDLgnOm1A8uWsR2QkAaCarSBevhQa5zuADwAsB2/Lcc2AFdhmobWY+C9fV4c2/kIaTioRhuisjLB+EmTPzkFwnAKTkihdDt+M4b/PBw05Vzz/MlB0mI02xEJaTioRhsg/sNgu4QtZ5gJW7HCHVPPYbwkygIwUgPvzzebLnXaDfLybGckpOGgGm2AnE/InsWE5vhkjisALHD5DCZywMSh/Y1VgRg5ku2AhDQoVKMNTkCA/ZVxrZH2Cc0NEBSFnxlwgmDAgLMMO0ztWrCdj5CGhmq0YYmPL7aYYJ15OAOt+uJJCpQNETgFJ5wwfxl31/Z+F6GqynZEQhoaqtEGpdjDa1zx6QBmhDPmPMCAiTiTAmUP2M3Ase3KO+FO37cnpObR1+9rnre3d2RkZIcOHezs7Op41fYuA/3QxxrnZ8IFwElMBnAGE69xTLMeT/hJrWkd5yENmLOzc0pKSmRkpLm5OdtZWEY1WvNMTEzGjh3L49X55eohIYtKtp/H2fOw3omlS7GTA0YSZQAWt3D7SW1FXechDZq9vb1QKDxy5AjbQdhHNVrzJCUlZWVl63qtERFlg412Fx8XgtsUBcuwQwjua3Q+gSlbOKtW2ObXdR7S0ElLSwOQoFvJUo02GGV/bZtY7HqBGbcLSwbh3nDcEF25tBUrlqufw9qnbAckpMGiU0wNQVkZJvpMusCM24oVv2J3b4TcwHAFZAPgtlPDixd052RCag/VaEPAz/qUVyYD4BaMiiEDoBdCDREIILt5O7pmiZBaRTVa/0VGyup27FoWBuA6RozFpSLIzoHzBVj+js1/zaRRSAipXXRstP6bO3dl+tLd+HUuDr1DB38Yd8TbRKis5mz+U+UQpkWynY+QBo5qtJ5LT199d+RWLHfAYSfML4G0BS4HYGRHvF3b7wZO3kGzZmxHJKSBo536+u3WlYLN+F0R6RuxlgNGBsWXYaGL8LfouEd1BzQ12Q5ISMNHNVqfffpk9OhPW5xKh+Jw3MhAKwBHMSsCOha4vHgmfVeUkLpAO/X12cSJAX7CC7BURcILdB+Gm9Y4/wc2WXC8z7VZLDkkmu18hDQKtDX6BWFhYT169Gjfvr2JiUlWVlaFZxcuXKiurs7hcDIyMliJ94/gYH9foQVzqTtevES3ddjwAt3X4E9zXDnXfJbkKTdISbEZj5BGg2r0C2bMmLFp06a4uDgdHZ0NGzZUeNbGxiY4OFh0JRyL/N1SxuJSN7z0h7E8ctdj/TpsABCH9rm3QjB0KLvxCGk8qEYrio2NTU5ONjMzAzB79uzz589XmGHQoEHKyspsRPuXT/lcPiRy0CIH/4zEPBIBcsgvgmyJDH3fnpC6QzVaUXx8vJqamuixurp6amoqn89nN9IXHDhgcM5BHR9joGWIwI9Qf4R+xvBvg+QgjWkqXWmIe0LqDp1iqgzDMD/wqg0bNuzatUtBQUFXVxeAlZVV7969azLWkycZizYOY26mQukPbPobq4bgThZ+UkJqEAxUNm2pyXUR8m+3b9/29/cHEBoamp+fX1BQQOONUo1WpKamFh8fL3ocHx+vpKT0vUOBrVu3ztrauhai/SPriNdw3IiFpi9MBuOuLsIn4owAvBtNLVSc/kKdjxVNGpWhQ4cO/ezIu5OTE3tZxAXt1FekqamprKx85coVAM7OzlZWVqLp/v7+mZmZrEYDAIaB6YVpYUz3k5g8GHcBWMHLA3YA7PkupROnsh2QkEaHavQLXFxc1qxZo6amFh4evm7dOtHEBQsWREVFAZg5c6aqqmpJSYmurm7d785wykpXq7hJoXQ1/kqFEoBCNDmC2VwIf+9wjr7jREjdo536L9DT0wsPD68wMSYmRvTAxcWlzhP9T3o6hg79EDm0DJKR6GKAoKswnQmXuxjsjqm2s2r0CCwh5NvQ1mi9smrVwSiDBThgiquXMDYWmp3x+p8O7fUac+awnY+QxohqtD5xOq2wgNlvBh8vWFng8kWM40KogkSjsc3x4AHq/gZQhBCq0Xrkwzv+0qI/m6LABTOlUArABL6Lsecj1Fe/mw62r6oipNGiGq032hVGneBNy4ecKa6Krlw6iwk7sXQw7u6b9pztdIQ0XlSj9YSPT3HPfscFUyXAD0WvEbjujDl28BiAB76ylk3Hm7Kdj5DGi87U1wdlZcUz5lsIL96EkQfsiiEzA8eeos8Q3Lkma9XUwxlt2rAdkZDGi2q0Hih6+Nwi3eUWjNwxdQLOAhCANxtHmuOTpLcXhg9hOyAhjRrt1NcDF65KX8eIAXgwEWdEU6bCvQ+e+sDM51ErdrMRQqhG6wE79buO2HsXg6fBVQAeHxJ28HiCvts5yy2n05h4hLCMdurFnqPj231+XgjiQXACUxhwSiDtCZvt3BXLjMKgqsp2PkIaO6rRmhcQEPDx40cNDQ1LS8vqLis09O1+fwMElUIqGL9sxu8nMRnADixb2j0Qrj41EJeQH3Ly5Mm0tLSwsDAaKI9qtOb179/fzMxMqiaGCXnr/tCACSyF1C0Y6SDiDCZOwmkvWHFlpRESAi4dkyGsMTU15fP5QqGQ7SDso3+HNa9Zs2aKiory8jVw1HLjjX6JUFmI/TqIACCJsh1YJoPi9UUrklPpsyNsUlBQUFRUbNKkCdtB2Ef/FMXakSGn+uHRBqw7gSkAEqEyAtelUOqnNpu+KkqImKCdenHF58PW9pln4ktsFIA3Da7pUHTGnFQo+WFU/4UWbOcjhPyDtkbF1ZEjDzwTR8GvLZJ8YfITspZhRyqU/Dmj+k/rhCVL2M5HCPkH1aiYun8kUtShQTAYBb+bGNYKGc2Q13pMPxw/Dh6P7YCEkH9QjYqjoiLYvFybD7nTmNQGyQC648U+LEqEyuS7s37odqWEkNpCNSqOZGXh2WW9HPIn4XQy2gB4ge6LsE8N8SdHneZw2M5HCPkM1aj4KS2FgwPz8hUDzht0MkCQH0YNw01ZFAXBoMMsQ7bzEUL+hWpU/Kxaddc5yoS5pooEG3i+QScT+MqiKIhj2GH1BAyh8ZwIES/0hScxU1Z299Cr0bimhvhAGCoiXQqlHrCzgleHw79h9my28xFCKqKtUfHy/EaGSdGFFsgJhKEyUngQuMHeApd349d9t7qynY4Q8gVUo+KlSw+pobidhLZ+GCWaEoWfH2CAOj6a9k1nNxsh5Itop16cMIz0zWubePv9BKNmwoUHQW+E/HNyiTdM0yaQ7XyEkC+gGhUny5dH7vAdhcCfkJUNhek43gI5TVEQBAPNdfY0tCgh4olqtOZ5eXmFhoZqa2tPnz79O14WH/9qV4ARArkQ3sGQ5+g5GSezoRAkMVzzzFZYWdVaXkJ+xL59+5KSkqKjo2m8UarRmmdlZWVtbf29r3p1KsxIeIMHQSAMO+FNJ7zhgLGDxwL+nmu9u8jVRlBCqmHRokUAnJyc2A7CPjrFJC4evpJPQ2t9PNZCjGjKUNxWRUI4dOPelLKbjRBSCapRcTFrcvE2LL+IcfZwE4CXAmUDBGVDIaDpOB2j1mynI4R8Fe3Ui4eLF4unz7sJdwAesCuFVAR0EqESgJG/LBsJCfqYCBFftDUqBmJiiifYW3w6cRPD1mFDc3zyhE0CVP04JvpzumPNGrbzEUIqQzXKvmJndwu+103G6ASmrMd6P4xqhrw2SO5g2A6HD9OmKCFijmqUfdYeYwKYEfuxcBJOA+iPh1dg/hYdje6tLypiOxwhpCpUo+yb8fNDKZSewJRcyAPgQ+IwHBhw7Fv7ysqyHY4QUhWqUbaVlVmoPJ2Fo4+hPxIBmWhpi1PnMH47Z/my6VlshyOEVI2Ou7EqJwcGBm5hPQ5hLoAn6KuJ2E9ovp2zfJmGF359xnY+QkjVaGuUVevWub3oOQPHjHArCAbyyP2E5r9j87LJqXj2DC1asJ2PEFI1qlE2ubpzZzAuw3DTG2OG4nYARsoj9xzGx7fuRR1KSH1BNcqalBQszP1TCqXHMEMWRQD64slKbHmHDr95D2Q7HSHkW1GNskZZGZeU53LAmOJqBloB8MaYddjQF0+cJ99nOx0h5FvRKaaaFxQUlJmZqaamNnr06MrmCw7uKx3WEW9foLsRbjli71wc6onnAVJm8pOD6yosIT/Iy8srIyMjNDSUBsqjGq15urq6xsbGspV/5/PcubyJs0fBPwo/y6A4HLozcEwfj/1lLOSP7YGGRh1lJeRH9e/fv7S0NDMzk+0g7KMarXktW7bUqLwHBYL8Rb+bMNeeovdZTJBB8ThcZMA5x50oH/kI7dvXUVBCqqFt27YAWtC5UDo2yoqCJy9Hp7k+hv4p2I7DRRP4noc1ABvh2U8xaWynI4R8H6pRFty6zbuPgd3xwgw+oinG8P8FwcH4xfeWNLvZCCHfi2qUBebTWjpzHJ6j5zhcLIZMGSTH49wDDNiDxRPmyLOdjhDyfejYaJ1jGNy4oScZIVea74dRY3FJBsXeGLOH8+sio0g6MEpIvUM1Wudmznx2PGw4boh+8ocxB8xeOC7s/wwnzrMbjRDyA2invm4FB4s6tBnyXqD7GvwJoB/n8fx+z3DvHtq0YTsfIeS70dZonQp1fzkcN+WRGwQDDbzfhD84YDYxf8wOnnmE4XA5bOcjhHw/2hqtU2dDO2ZDwQS+GngvmmIHj2bI8xGYJLznsxqNEPKDaGu0Tm2bG5f/5NVBzGuGvC1YGQ1tAwTJouiWmr26pi/b6QghP4K2Rr9bWFhYjx492rdvb2JikpX1PQPUC4WczAwT+HLAbMWK2ThigCA+JAJh2O23UbWWlxBSu6hGv9uMGTM2bdoUFxeno6OzYcOGb31ZSQkMDa8uDbKCFxdCAEcxSwBeIHdY11XmWLCgFhMTQmoT1ej3iY2NTU5ONjMzAzB79uzz57/5K0qHDvncaW6JC93x4iW6/YwoAH9iTdctk7F5Mzh0domQ+opq9PvEx8erqamJHqurq6empvL533RqyMc1wwpePRB2HSM643UQDLogcjH23PGneygTUr/RKaYfxzDMF6dv2LBh165dCgoKurq6AKysrHR1ey95PbsUUluxogVyACghdSPWWsFracjEpwxtjJJ64/bt2/7+/gBCQ0Pz8/MLCgpovFEw5Hu8e/euTZs2oscxMTHlj8uFhIR4enr+94XxxjO18FYeOY+gzwA3MEwWhbqc8LTJS2o9NCG15sCBAx8/fmQ7Bctop/77aGpqKisrX7lyBYCzs7OVldU3vlB1rtl1jGiBnNG4dghzx+KSBt4HYKSig2Vt5iWE1Dqq0e/m4uKyZs0aNTW18PDwdevWfevLzM3bL7cO4hg2Q948HFRDfKDkSOW9q9C/f22GJYTUOjo2+t309PTCw8N/5JVbt7a3tr597Oxcb2NXa1/leTfQqVNNpyOE1DWq0brVu7dG795+hwB0ZzsKIaRm0E49C/z8/NiOAAAJCQkRERFspwDE5g1hGCYgIIDtFADw6NGjnJwctlMAYvPRiDmqURbcuXOH7QgAkJKS8ubNG7ZTAGLzhgiFwrt377KdAgCeP3/+6dMntlMAYvPRiDmqUUIIqRaq0Rr27t2758+fVz7P48ePK5/h9OnT8fHxlc+zbdu2ymfYunVr5TPExMRUGbXKhZw6dSohIaGSGRiG2b59e+ULqfINCQ8PF33luxJVRvXw8EhMTKxkBqFQWGWSKtcSFhZW5ZGBKhcSHBycmppayQwCgWDHjh3VXMvz58+vX79e+TxVviFPnjxJS2vst7OlGq1hpaWlxcXFlc9TUlJS+Qx5eXmVX2PKMEx2dnblC6lyhrKysiqTVLmQGolaZYySkpKCgoLK56mRqFV+dt/yuxQWFlZzIcXFxQKBoJIZGIap8uBpjUT9lo/mG6+HbsA4zFeuaCQ/xt/f/+DBg7/88ksl83h6etrY2FQyQ1RUlKqqarNmzSqZJzg4uPK1VDlDdHR0RkZG/0q/uFrlQiIjI9XU1CqJyjBMSEhInz59KllIlW9Ienp6fn5++0rv9/ctUdXV1eXk5L42g1AovHjxYuWXVFS5lrS0tMLCQg0NjeosJCAgoEePHkpKSpVEDQ0NrfxdrZGoVX40d+7c2bt3b5cuXSqZp8GjGq15BQUF9K6SRoLL5TZp0oTtFCyjGiWEkGqhY6OEEFItVKOEEFItVKOEEFItVKOEEFItVKOEEFItVKN16sdvzlxtaWlppqamysrKampq9vb25d8zZzHSggULOBxO+dfI6z5JQUHBtGnTFBUVVVRUyi8KYuUNOXHiRJcuXbp169anT58nT57UcZKFCxeqq6tzOJyMjIzyiV9cO4t/LWKNpVH3Gyk9Pb0rV64wDLN8+fJFixbV5aqTk5N9fHzKysqKi4tNTU3//PNPdiPduXPHzs6Ox+NlZ2ezlWTOnDlTpkwRXTIUFxfHVoyioiIpKamYmBiGYTw8PPT19es4yd27d5OTk6WlpdPT08snfnHtLP4BizOq0bpT5X2c6szWrVvnzJnDYqTCwkJ9ff20tLTyGq37JJ8+fWrSpElGRsbnE1l5Q/h8vqqq6osXLxiGOXz4sKWlJStJPq/RL65dfP6AxQ0N21x3vnhzZgmJuv4IioqK3N3dRfuwbEVas2bN3LlzFRUVy6fUfZLY2FglJaUNGzbcunVLWVl59+7durq6rLwhPB7P1dV14MCB8vLyPB7v3r17YPuv5YtrF5M/YDFEx0bZwbB08Rifz580aZKpqemoUaPYivTkyZNXr15NmTLlazPUTRI+nx8XF9e7d+9Xr17NmjXrv5fS19kbkp2dPWXKlNu3b8fHx69evXrSpElsJfmiL66d3Ujihmq07qipqZUPfxcfH6+kpFTH/ycXCoVTp05VVlYuH0KNlUj37t0LDw/X0NDQ0NAQCAS6urrR0dF1n0S0YTVx4kQA48ePj42Nzc/PZ+UNiYiIUFRU1NPTA2Bra/vgwQOBQMDuX8sX1876H7DYohqtOz98c+YawTDMjBkzJCUlDx48yG6kZcuWJSUlvX///v379zweLzw8XFtbu+6TtG7desiQIaKbZAQEBKipqcnJybHyhnTs2PH9+/fR0dEArly50qlTJx6Px+5fyxfXzm4kscbeYdnGKDQ0VEdHR1VVdeTIkRVObtQ20ddo2rRpo6KioqKiMnXqVNYjMQzz+Zn6uk/y+vVrfX39jh07ir5mxFYMhmE8PDw6d+6sra3dp0+f4ODgOk4yY8YMFRUV0Z+HmZlZJWtn969FbNEIT4QQUi20U08IIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdVCNUoIIdXyf8nSrbhsRxncAAAAAElFTkSuQmCC",
                            "text": [
                                "FramedPlot(...)"
                            ],
                            "metadata": {}
                        },
                        {
                            "output_type": "stream",
                            "stream": "stderr",
                            "text": [
                                "Warning: imported binding for transpose overwritten in module __anon__\n"
                            ]
                        }
                    ],
                    "input": [
                        "using Winston\n",
                        "\n",
                        "p = FramedPlot(aspect_ratio=1)\n",
                        "\n",
                        "p1 = Points(t_steps, num_result, color = \"red\", kind=\"filled circle\")\n",
                        "setattr(p1, label=\"CVODE\")\n",
                        "add(p, p1)\n",
                        "\n",
                        "p2 = Points(t_steps, analytical_result, color = \"blue\", alpha=0)\n",
                        "setattr(p2, label=\"Analytical\")\n",
                        "add(p, p2)\n",
                        "\n",
                        "l = Legend(0.1, .9, {p1, p2})\n",
                        "\n",
                        "add(p, l)"
                    ],
                    "language": "python",
                    "prompt_number": 9
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "Closing Thoughts\n",
                        "----------------\n",
                        "Julia is an excellent language, and even though it's very new, it still feels quite polished.  For this very simple\n",
                        "example, we used a pure Julia package (Winston) a pure Python package (SymPy) and an advanced C library (Sundials)\n",
                        "with very little overhead and almost no boiler plate code.\n",
                        "\n",
                        "I'll definitely use it more, and try to show off some its cooler features (meta-programming) later on."
                    ]
                }
            ]
        }
    ],
    "cells": [],
    "metadata": {
        "language": "Julia",
        "name": "",
        "signature": "sha256:1286369f8c280382cb1bbe934370c9060ac0c2e09dfdc5e17a6d39322ffcdd4f"
    },
    "nbformat": 3,
    "nbformat_minor": 0
}