dfba.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\n",
    "\"Contains functions to perform dynamic FBA.\"\n",
    "\n",
    "from __future__ import absolute_import\n",
    "\n",
    "from optlang.symbolics import Zero\n",
    "\n",
    "from scipy.integrate import solve_ivp\n",
    "import pandas as pd\n",
    "\n",
    "def dynamic_fba(model_and_rxns, kinetic_eq):\n",
    "    \"\"\"\n",
    "    Perform dynamic FBA to simulate batch growth of an organism.\n",
    "    \n",
    "    Dynamic FBA finds application in industrial fermentation process\n",
    "    and modeling community growth of microorganisms, among others. It\n",
    "    is an important tool which enhances FBA by adding reaction and\n",
    "    metabolite dynamics. It can be be performed using mainly 3\n",
    "    techniques: (1) Static Optimization Approach (SOA), \n",
    "    (2) Dynamic Optimization Approach, and (3) Direct Approach.\n",
    "    This algorithm uses the Direct Approach which circumvents the\n",
    "    inaccuracy of SOA and the Non-Linear Programming (NLP) complexity\n",
    "    found in DOA. This approach employs Lexicographic Linear Programming\n",
    "    to provide unique flux distribution and also solves the LP\n",
    "    feasibility problem when integrating in DAE.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    model_and_rxns: dict(model: cobra.Model, biomass: cobra.Reaction.id,\n",
    "                         exchanges: tuple(cobra.Reaction.id))\n",
    "        The model(s) to perform dynamic FBA on and the metabolites to track.\n",
    "    kinetic_eq: func()\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    pandas.DataFrame\n",
    "        \n",
    "    References\n",
    "    ----------\n",
    "    .. [1] Gomez et al.: DFBAlab: a fast and reliable MATLAB code for dynamic flux balance analysis.\n",
    "           BMC Bioinformatics 2014 15:409.\n",
    "    \"\"\"\n",
    "    model = model_and_rxns[\"model\"]\n",
    "    biomass_rxn_id = model_and_rxns[\"biomass\"]\n",
    "    exchanges_rxn_id = model_and_rxns[\"exchanges\"]\n",
    "    rxn_list = [biomass_rxn_id]\n",
    "    rxn_list.append(exchanges_rxn_id)\n",
    "    # 1. LLP\n",
    "    y = lexicographic_lp(model, rxn_list)\n",
    "    # 2. ODE integration\n",
    "    # a. modify \"y\" so that if fits into the kinetic equation\n",
    "    # of the form \"kinetic_eq(t, y)\"\n",
    "    # b. set y0 to the initial concentrations of the metabolites\n",
    "    solution = solve_ivp(kinetic_eq, t_span, y0, method='LSODA')\n",
    "    return solution\n",
    "\n",
    "def lexicographic_lp(model, rxn_list):\n",
    "    \"\"\"\n",
    "    Perform Lexicographic Linear Programming.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    model: cobra.Model\n",
    "        The model to perform lexicographic LP on.\n",
    "    rxn_list: list(cobra.Reaction.id, tuple(cobra.Reaction.id))\n",
    "        The list containing the reactions to be considered as objectives.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    pandas.Series\n",
    "    \"\"\"\n",
    "    biomass_rxn_id, ex_rxn_ids = rxn_list\n",
    "    \n",
    "    with model:\n",
    "        # LP feasibility\n",
    "        obj_vars = []\n",
    "        prob = model.problem\n",
    "        for met in model.metabolites:\n",
    "            s = prob.Variable(\"s_\" + met.id, lb=0)\n",
    "            beta = prob.Variable(\"beta_\" + met.id, lb=0)\n",
    "            s_eq_beta_const = \\\n",
    "                prob.Constraint(s - beta,\n",
    "                                name=\"s_equal_beta_\" + met.id,\n",
    "                                ub=0.0, lb=0.0)\n",
    "            model.add_cons_vars([s, beta, s_eq_beta_const])\n",
    "            model.constraints[met.id].set_linear_coefficients({s: 1.0, beta: -1.0})\n",
    "            obj_vars.append(s)\n",
    "        model.objective = prob.Objective(Zero, sloppy=True, direction=\"min\")\n",
    "        model.objective.set_linear_coefficients({v: 1.0 for v in obj_vars})\n",
    "        sol_feasibility = model.slim_optimize()\n",
    "        feasibility_const = prob.Constraint(model.objective.expression,\n",
    "                                            name=\"fixed_feasibility\",\n",
    "                                            ub=sol_feasibility, lb=sol_feasibility)\n",
    "        model.add_cons_vars([feasibility_const])\n",
    "        \n",
    "        # Biomass\n",
    "        model.objective = model.reactions.get_by_id(biomass_rxn_id)\n",
    "        model.objective_direction = \"max\"\n",
    "        sol_biomass = model.slim_optimize()\n",
    "        biomass_const = prob.Constraint(model.objective.expression,\n",
    "                                        name=\"fixed_biomass\",\n",
    "                                        ub=sol_biomass, lb=sol_biomass)\n",
    "        model.add_cons_vars([biomass_const])\n",
    "        \n",
    "        # Exchanges\n",
    "        for rxn_id in ex_rxn_ids:\n",
    "            model.objective = model.reactions.get_by_id(rxn_id)\n",
    "            model.objective_direction = \"min\"\n",
    "            sol = model.slim_optimize()\n",
    "            exchange_const = prob.Constraint(model.objective.expression,\n",
    "                                             name=\"fixed_\" + rxn_id,\n",
    "                                             ub=sol, lb=sol)\n",
    "            model.add_cons_vars([exchange_const])\n",
    "\n",
    "        # Indexing list\n",
    "        loc_list = [biomass_rxn_id]\n",
    "        loc_list.extend(list(ex_rxn_ids))\n",
    "        \n",
    "        sol =  model.optimize()\n",
    "        return pd.Series(sol.fluxes.loc[loc_list])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cobra.test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = cobra.test.create_test_model('textbook')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = lexicographic_lp(model, ['Biomass_Ecoli_core', ('EX_glc__D_e',)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = np.array([0.001, 10.4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_max = np.array([1000, 10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 0.015"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mm_kin_eq(t, y, **kwargs):\n",
    "    return y_max * y / (k  + y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = solve_ivp(fun=lambda t, y: mm_kin_eq(y_max, k), t_span=(0, 36000), y0=y0, t_eval=[i for i in range(36000)], vectorized=True, method='LSODA')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  message: 'The solver successfully reached the end of the integration interval.'\n",
       "     nfev: 27\n",
       "     njev: 0\n",
       "      nlu: 0\n",
       "      sol: None\n",
       "   status: 0\n",
       "  success: True\n",
       "        t: array([    0,     1,     2, ..., 35997, 35998, 35999])\n",
       " t_events: None\n",
       "        y: array([[1.000000e-03, 5.000010e+02, 1.000001e+03, ..., 1.799850e+07,\n",
       "        1.799900e+07, 1.799950e+07],\n",
       "       [1.040000e+01, 1.540000e+01, 2.040000e+01, ..., 1.799954e+05,\n",
       "        1.800004e+05, 1.800054e+05]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(sol.y.transpose(), index=sol.t, columns=tuple(['Biomass', 'Glucose']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Biomass</th>\n",
       "      <th>Glucose</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.000000e-03</td>\n",
       "      <td>10.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5.000010e+02</td>\n",
       "      <td>15.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.000001e+03</td>\n",
       "      <td>20.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.500001e+03</td>\n",
       "      <td>25.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2.000001e+03</td>\n",
       "      <td>30.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2.500001e+03</td>\n",
       "      <td>35.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>3.000001e+03</td>\n",
       "      <td>40.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>3.500001e+03</td>\n",
       "      <td>45.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>4.000001e+03</td>\n",
       "      <td>50.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>4.500001e+03</td>\n",
       "      <td>55.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>5.000001e+03</td>\n",
       "      <td>60.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>5.500001e+03</td>\n",
       "      <td>65.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>6.000001e+03</td>\n",
       "      <td>70.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>6.500001e+03</td>\n",
       "      <td>75.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>7.000001e+03</td>\n",
       "      <td>80.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>7.500001e+03</td>\n",
       "      <td>85.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>8.000001e+03</td>\n",
       "      <td>90.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>8.500001e+03</td>\n",
       "      <td>95.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>9.000001e+03</td>\n",
       "      <td>100.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>9.500001e+03</td>\n",
       "      <td>105.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>1.000000e+04</td>\n",
       "      <td>110.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>1.050000e+04</td>\n",
       "      <td>115.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>1.100000e+04</td>\n",
       "      <td>120.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>1.150000e+04</td>\n",
       "      <td>125.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>1.200000e+04</td>\n",
       "      <td>130.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>1.250000e+04</td>\n",
       "      <td>135.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>1.300000e+04</td>\n",
       "      <td>140.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>1.350000e+04</td>\n",
       "      <td>145.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>1.400000e+04</td>\n",
       "      <td>150.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>1.450000e+04</td>\n",
       "      <td>155.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35970</th>\n",
       "      <td>1.798500e+07</td>\n",
       "      <td>179860.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35971</th>\n",
       "      <td>1.798550e+07</td>\n",
       "      <td>179865.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35972</th>\n",
       "      <td>1.798600e+07</td>\n",
       "      <td>179870.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35973</th>\n",
       "      <td>1.798650e+07</td>\n",
       "      <td>179875.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35974</th>\n",
       "      <td>1.798700e+07</td>\n",
       "      <td>179880.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35975</th>\n",
       "      <td>1.798750e+07</td>\n",
       "      <td>179885.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35976</th>\n",
       "      <td>1.798800e+07</td>\n",
       "      <td>179890.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35977</th>\n",
       "      <td>1.798850e+07</td>\n",
       "      <td>179895.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35978</th>\n",
       "      <td>1.798900e+07</td>\n",
       "      <td>179900.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35979</th>\n",
       "      <td>1.798950e+07</td>\n",
       "      <td>179905.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35980</th>\n",
       "      <td>1.799000e+07</td>\n",
       "      <td>179910.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35981</th>\n",
       "      <td>1.799050e+07</td>\n",
       "      <td>179915.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35982</th>\n",
       "      <td>1.799100e+07</td>\n",
       "      <td>179920.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35983</th>\n",
       "      <td>1.799150e+07</td>\n",
       "      <td>179925.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35984</th>\n",
       "      <td>1.799200e+07</td>\n",
       "      <td>179930.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35985</th>\n",
       "      <td>1.799250e+07</td>\n",
       "      <td>179935.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35986</th>\n",
       "      <td>1.799300e+07</td>\n",
       "      <td>179940.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35987</th>\n",
       "      <td>1.799350e+07</td>\n",
       "      <td>179945.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35988</th>\n",
       "      <td>1.799400e+07</td>\n",
       "      <td>179950.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35989</th>\n",
       "      <td>1.799450e+07</td>\n",
       "      <td>179955.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35990</th>\n",
       "      <td>1.799500e+07</td>\n",
       "      <td>179960.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35991</th>\n",
       "      <td>1.799550e+07</td>\n",
       "      <td>179965.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35992</th>\n",
       "      <td>1.799600e+07</td>\n",
       "      <td>179970.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35993</th>\n",
       "      <td>1.799650e+07</td>\n",
       "      <td>179975.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35994</th>\n",
       "      <td>1.799700e+07</td>\n",
       "      <td>179980.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35995</th>\n",
       "      <td>1.799750e+07</td>\n",
       "      <td>179985.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35996</th>\n",
       "      <td>1.799800e+07</td>\n",
       "      <td>179990.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35997</th>\n",
       "      <td>1.799850e+07</td>\n",
       "      <td>179995.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35998</th>\n",
       "      <td>1.799900e+07</td>\n",
       "      <td>180000.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35999</th>\n",
       "      <td>1.799950e+07</td>\n",
       "      <td>180005.4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>36000 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            Biomass   Glucose\n",
       "0      1.000000e-03      10.4\n",
       "1      5.000010e+02      15.4\n",
       "2      1.000001e+03      20.4\n",
       "3      1.500001e+03      25.4\n",
       "4      2.000001e+03      30.4\n",
       "5      2.500001e+03      35.4\n",
       "6      3.000001e+03      40.4\n",
       "7      3.500001e+03      45.4\n",
       "8      4.000001e+03      50.4\n",
       "9      4.500001e+03      55.4\n",
       "10     5.000001e+03      60.4\n",
       "11     5.500001e+03      65.4\n",
       "12     6.000001e+03      70.4\n",
       "13     6.500001e+03      75.4\n",
       "14     7.000001e+03      80.4\n",
       "15     7.500001e+03      85.4\n",
       "16     8.000001e+03      90.4\n",
       "17     8.500001e+03      95.4\n",
       "18     9.000001e+03     100.4\n",
       "19     9.500001e+03     105.4\n",
       "20     1.000000e+04     110.4\n",
       "21     1.050000e+04     115.4\n",
       "22     1.100000e+04     120.4\n",
       "23     1.150000e+04     125.4\n",
       "24     1.200000e+04     130.4\n",
       "25     1.250000e+04     135.4\n",
       "26     1.300000e+04     140.4\n",
       "27     1.350000e+04     145.4\n",
       "28     1.400000e+04     150.4\n",
       "29     1.450000e+04     155.4\n",
       "...             ...       ...\n",
       "35970  1.798500e+07  179860.4\n",
       "35971  1.798550e+07  179865.4\n",
       "35972  1.798600e+07  179870.4\n",
       "35973  1.798650e+07  179875.4\n",
       "35974  1.798700e+07  179880.4\n",
       "35975  1.798750e+07  179885.4\n",
       "35976  1.798800e+07  179890.4\n",
       "35977  1.798850e+07  179895.4\n",
       "35978  1.798900e+07  179900.4\n",
       "35979  1.798950e+07  179905.4\n",
       "35980  1.799000e+07  179910.4\n",
       "35981  1.799050e+07  179915.4\n",
       "35982  1.799100e+07  179920.4\n",
       "35983  1.799150e+07  179925.4\n",
       "35984  1.799200e+07  179930.4\n",
       "35985  1.799250e+07  179935.4\n",
       "35986  1.799300e+07  179940.4\n",
       "35987  1.799350e+07  179945.4\n",
       "35988  1.799400e+07  179950.4\n",
       "35989  1.799450e+07  179955.4\n",
       "35990  1.799500e+07  179960.4\n",
       "35991  1.799550e+07  179965.4\n",
       "35992  1.799600e+07  179970.4\n",
       "35993  1.799650e+07  179975.4\n",
       "35994  1.799700e+07  179980.4\n",
       "35995  1.799750e+07  179985.4\n",
       "35996  1.799800e+07  179990.4\n",
       "35997  1.799850e+07  179995.4\n",
       "35998  1.799900e+07  180000.4\n",
       "35999  1.799950e+07  180005.4\n",
       "\n",
       "[36000 rows x 2 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x11a259eb8>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEDCAYAAAAyZm/jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VAXWx/HvSaF3EjohoUNIQglEQFFRKYqioCviKlbUta27r4JiQWyAru6uuiIoIrsqu1QRVNBFRRdEQCGNFkJL6AklEELaef+YAQdMYBImmZnkfJ4nT+7c+pubmZObeyfniqpijDGm8gjwdgBjjDHlywq/McZUMlb4jTGmkrHCb4wxlYwVfmOMqWSs8BtjTCXjs4VfRKaLyH4RSXRj3jdEZJ3za7OIHC6PjMYY44/EVz/HLyL9gGPATFXtUoLlHga6qepdZRbOGGP8mM8e8avqciDTdZyItBGRL0VkrYh8LyIdi1j0FuCTcglpjDF+KMjbAUpoKnC/qm4RkTjgH0D/UxNFpBUQASzzUj5jjPF5flP4RaQW0AeYLSKnRlc9a7YRwBxVLSjPbMYY40/8pvDjOC11WFW7nmOeEcCD5ZTHGGP8ks+e4z+bqh4FtonITQDiEHNquvN8f31gpZciGmOMX/DZwi8in+Ao4h1EJE1E7gZuBe4WkfVAEjDUZZERwCz11Y8pGWOMj/DZj3MaY4wpGz57xG+MMaZs+OTF3ZCQEA0PD/d2DGOM8Rtr1649qKqh7szrk4U/PDycNWvWeDuGMcb4DRHZ4e68dqrHGGMqGSv8xhhTyVjhN8aYSsYnz/EXJS8vj7S0NHJycrwdxW9Uq1aNFi1aEBwc7O0oxhgf4jeFPy0tjdq1axMeHo5Lrx5TDFUlIyODtLQ0IiIivB3HGOND/OZUT05ODg0bNrSi7yYRoWHDhvYXkjHmN/ym8ANW9EvI9pcxpih+VfiNMaaiWr09kynfbS2XbVnhL4HAwEC6du1KTEwM3bt3Z8WKFQDs3r2bG2+80cvpjDH+6NjJfJ79NJGbpqzk41U7yc7NL/Nt+s3FXV9QvXp11q1bB8CSJUt48skn+e6772jWrBlz5szxcjpjjL/5bvMBnpqXwO4jJ7izbzj/N6ADNaqUfVm2I/5SOnr0KPXr1wdg+/btdOniuB98Tk4Od955J1FRUXTr1o1vvvkGgBkzZnD99ddz1VVXER4ezltvvcXrr79Ot27duOiii8jMdNxeeNq0afTs2ZOYmBiGDx9OdnY2ALNnz6ZLly7ExMTQr18/AJKSkujVqxddu3YlOjqaLVu2lPduMMaUwqHjufzpP+sYNf0nqlcJZM79fXju2khqVi2fY/HzbkVEpgNDgP2q2qWI6Y/j6JN/an2dgFBVzRSR7UAWUADkq2qsJ0I//1kSybuPemJVp3VuVofnro085zwnTpyga9eu5OTksGfPHpYt++2tfd9++21EhISEBDZu3MiAAQPYvHkzAImJifzyyy/k5OTQtm1bJk2axC+//MJjjz3GzJkz+eMf/8iwYcO49957AXj66ad5//33efjhh5kwYQJLliyhefPmHD58GIApU6bw6KOPcuutt5Kbm0tBgd1x0hhfpqp8kbiXZz9N5HB2Hg/3b8tD/dtSNSiwXHO48+tlBvAWMLOoiar6KvAqgIhcCzymqpkus1yuqgcvMKdPcD3Vs3LlSm6//XYSExPPmOeHH37g4YcfBqBjx460atXqdOG//PLLqV27NrVr16Zu3bpce+21AERFRREfHw84fjk8/fTTHD58mGPHjjFw4EAA+vbtyx133MHvfvc7hg0bBkDv3r156aWXSEtLY9iwYbRr167sd4IxplT2H83hmU8TWZK0j6jmdZl5Vxydm9XxSpbzFn5VXS4i4W6u7xbgkwsJ5I7zHZmXh969e3Pw4EEOHDjg9jJVq/56b/iAgIDTjwMCAsjPd1zQueOOO1iwYAExMTHMmDGDb7/9FnAc3a9atYrFixfTo0cP1q5dy8iRI4mLi2Px4sVcffXVvPvuu/Tv399zT9IYc8FUldlr03hxUTIn8wsZO7gj91wcQVCg9860e2zLIlIDGATMdRmtwFIRWSsio8+z/GgRWSMia0pSTL1l48aNFBQU0LBhwzPGX3LJJXz00UcAbN68mZ07d9KhQwe315uVlUXTpk3Jy8s7vR6ArVu3EhcXx4QJEwgNDWXXrl2kpqbSunVrHnnkEYYOHXr6rwZjjG/YlZnNbe//xBNz4unYpA5fPHoJ91/axqtFHzz7qZ5rgf+ddZrnYlVNF5FGwFcislFVlxe1sKpOBaYCxMbG+uT9IE+d4wfHb/EPP/yQwMAzz8394Q9/4IEHHiAqKoqgoCBmzJhxxpH++bzwwgvExcURGhpKXFwcWVlZADz++ONs2bIFVeWKK64gJiaGSZMm8c9//pPg4GCaNGnCU0895bkna4wptYJC5cMV23l1ySYCA4QXru/Crb3CCAjwjX+qdOueu85TPYuKurjrMs98YLaqflzM9PHAMVV97Xzbi42N1bNvxLJhwwY6dep03qzmTLbfjClfW/ZlMWZuPD/vPMxlHUJ5+YYomtWrXubbFZG17n6AxiNH/CJSF7gU+L3LuJpAgKpmOYcHABM8sT1jjPE1eQWFTPl2K28uS6Fm1UD+enNXhnZt5pOtU9z5OOcnwGVAiIikAc8BwQCqOsU52w3AUlU97rJoY2C+80kHAR+r6peei26MMb4hIe0Ij89Zz8a9WQyJbsr46yIJqeX+Kd7y5s6nem5xY54ZOD726TouFYgpbTBjjPF1OXkFvPH1ZqYtTyWkVlWm3taDAZFNvB3rvKxlgzHGlMKq1AzGzktg28Hj3NKrJWMHd6Judf+46ZEVfmOMKYGsnDwmfbmRf/24k7AGNfj4njj6tA3xdqwSscJvjDFu+mbjfp6an8C+ozncc3EEfxrQvlyaqnmaNWkrgX379jFy5Ehat25Njx496N27N/Pnz+fbb79lyJAh3o5njCkjmcdz+eOsX7hzxmpqVQ1i7gN9eHpIZ78s+mBH/G5TVa6//npGjRrFxx87/lVhx44dLFy48HSXTmNMxaKqLIrfw/iFSRw5kcejV7TjD5e3Kfemap5mR/xuWrZsGVWqVOH+++8/Pa5Vq1anG7KdMn78eF577df/UevSpQvbt28HYObMmURHRxMTE8Ntt90GOFo69+/fn+joaK644gp27twJFN2GuaCggMcff5yePXsSHR3Nu+++W5ZP2ZhKbd/RHO6duZaHP/mF5vWrs+iRi3nsqvZ+X/TBX4/4vxgLexM8u84mUTB4YrGTk5KS6N69e6lXn5SUxIsvvsiKFSsICQk53X//4YcfZtSoUYwaNYrp06fzyCOPsGDBgiLbML///vvUrVuX1atXc/LkSfr27cuAAQOIiIgodS5jzJlUlX+v3sVLn28gr6CQcVd34s6+4V7vr+NJFeeZlLMHH3yQmJgYevbs6db8y5Yt46abbiIkxHH1v0GDBoCjvfPIkSMBuO222/jhhx+AX9swT5s27XSf/aVLlzJz5ky6du1KXFwcGRkZdvMVYzxoR8Zxbn1vFWPnJRDZrA5fPtqPe/u1rlBFH/z1iP8cR+ZlJTIykrlzf208+vbbb3Pw4EFiY89sjREUFERhYeHpxzk5OaXaXlFtmFWVN99883SPfmOMZxQUKh/8bxuvLd1EcEAAL98QxYieLX2mqZqnVaxfY2Wof//+5OTk8M4775wed+q2iK7Cw8P5+eefAfj555/Ztm3b6eVnz55NRkYGwOlTPX369GHWrFkAfPTRR1xyySVA0W2YBw4cyDvvvENeXh7gaPt8/PhxjDGlt2lvFsPeWcGLizfQt00IS//Uj5FxvtNJsyz45xG/F4gICxYs4LHHHmPy5MmEhoZSs2ZNJk2adMZ8w4cPZ+bMmURGRhIXF0f79u0Bx18M48aN49JLLyUwMJBu3boxY8YM3nzzTe68805effVVQkND+eCDD4Ci2zBHR0ezfft2unfvjqoSGhrKggULyn1fGFMR5OYX8o9vU3j7mxRqVwvmbyO6cl2MbzZV8zS32jKXN2vL7Dm234z5rfW7DvPEnHg27ctiaNdmPDukMw19uKmaO8q9LbMxxviDE7kFvP7VJt7/YRuNalfj/VGxXNGpsbdjlTsr/MaYSmHF1oOMnZvAzsxsRsaFMXZwR+pU84+map7mV4VfVSvF+TdP8cXTeMaUt6M5ebzy+UY++WknrRrW4JN7L6J3m4bnX7AC85vCX61aNTIyMmjYsKEVfzeoKhkZGVSrVs3bUYzxmq+T9zFuQQIHsk4yul9rHruyPdWr+P9/3l4ovyn8LVq0IC0tjQMHDng7it+oVq0aLVq08HYMY8pdxrGTPP9ZMgvX76Zjk9pMvS2WmJb1vB3LZ/hN4Q8ODrbWBMaYc1JVFq7fzfiFSRw7mc+frmrP/Ze2oUqQ/cuSq/PuDRGZLiL7RSSxmOmXicgREVnn/HrWZdogEdkkIikiMtaTwY0xxtWeIye458M1PDprHa0a1mTxI5fwyBXtrOgXwZ0j/hnAW8DMc8zzvaqe0ZBeRAKBt4GrgDRgtYgsVNXkUmY1xpjfKCxUPlm9k1c+30hBofLMkM7c0SecwAr8n7cXyp2brS8XkfBSrLsXkOK86ToiMgsYCljhN8Z4xLaDxxk7N55V2zLp27Yhr9wQTVjDGt6O5fM8dY6/t4isB3YD/6eqSUBzYJfLPGlAXHErEJHRwGiAsLAwD8UyxlRE+QWFTP/fNv6ydDNVggKYNDyK38W2tE/8uckThf9noJWqHhORq4EFQLuSrkRVpwJTwdGywQO5jDEV0IY9RxkzN574tCNc1bkxL17fhcZ17GPLJXHBhV9Vj7oMfy4i/xCRECAdaOkyawvnOGOMKbGT+QW8vSyFf3y7lbrVg3lrZDeuiWpqR/mlcMGFX0SaAPtUVUWkF45PCmUAh4F2IhKBo+CPAEZe6PaMMZXPzzsPMWZOPFv2H2NYt+Y8M6Qz9WtW8XYsv3Xewi8inwCXASEikgY8BwQDqOoU4EbgARHJB04AI9TRKyBfRB4ClgCBwHTnuX9jjHFLdm4+ry3ZzAcrttG0TjU+uLMnl3do5O1Yfs9v2jIbYyqX/6UcZOy8eHZlnuC2i1rxxKAO1K6kTdXcYW2ZjTF+68iJPF5evIF/r9lFREhN/j36IuJaV+6map5mhd8Y4zOWJu3l6QWJZBzP5f5L2/DHK9tRLdiaqnmaFX5jjNcdyDrJ+M+SWBy/h05N6/D+qJ5Etajr7VgVlhV+Y4zXqCrzf0lnwqJksk8W8PjADozu15rgQOuvU5as8BtjvCL98AnGzU/g200H6B5Wj8k3RtO2UW1vx6oUrPAbY8pVYaHy0aodTPxiIwqMv7Yzt/W2pmrlyQq/MabcpB44xti5Cfy0PZNL2oXw8g1RtGxgTdXKmxV+Y0yZyy8oZNr323jj681UCwrg1RujubFHC2u34CVW+I0xZSpp9xHGzI0nMf0ogyKbMGFoJI2sqZpXWeE3xpSJnLwC3ly2hSnfpVK/RhXeubU7g6OaejuWwQq/MaYMrN2RyRNz4tl64DjDu7fgmSGdqFfDmqr5Civ8xhiPOX4yn1eXbOLDldtpVrc6H97Vi0vbh3o7ljmLFX5jjEcs33yAJ+clsPvICW6/qBWPD+pIrapWYnyR/VSMMRfkcHYuLy7ewJy1abQOrcns+3oTG97A27HMOVjhN8aU2hcJe3jm0yQOZefy4OVteLi/NVXzB1b4jTEltj8rh+c+TeKLxL1ENqvDh3f1JLKZNVXzF1b4jTFuU1XmrE3jxcUbOJFXwBODOnDvJdZUzd9Y4TfGuGVXZjZPzU/g+y0H6Rlen4nDo2kTWsvbsUwpuHPP3enAEGC/qnYpYvqtwBhAgCzgAVVd75y23TmuAMh397ZgxhjfUViozFy5nclLNiHAhKGR/D6uFQHWVM1vuXPEPwN4C5hZzPRtwKWqekhEBgNTgTiX6Zer6sELSmmM8YqU/ccYOzeeNTsOcWn7UF66oQst6ltTNX933sKvqstFJPwc01e4PPwRaHHhsYwx3pRXUMjU5an87est1KgayOu/i+GGbs2tqVoF4elz/HcDX7g8VmCpiCjwrqpO9fD2jDEelph+hCfmxJO85yjXRDVl/HWRhNau6u1YxoM8VvhF5HIchf9il9EXq2q6iDQCvhKRjaq6vJjlRwOjAcLCwjwVyxjjppy8Av723y1MXZ5Kg5pVmPL7Hgzq0sTbsUwZ8EjhF5Fo4D1gsKpmnBqvqunO7/tFZD7QCyiy8Dv/GpgKEBsbq57IZYxxz+rtmYyZE0/qweP8LrYF467uTN0awd6OZcrIBRd+EQkD5gG3qepml/E1gQBVzXIODwAmXOj2jDGec+xkPpO/3MjMlTtoUb86/7o7jovbhXg7lilj7nyc8xPgMiBERNKA54BgAFWdAjwLNAT+4bzwc+pjm42B+c5xQcDHqvplGTwHY0wpfLNpP+PmJbDnaA539Y3gzwPaU9OaqlUK7nyq55bzTL8HuKeI8alATOmjGWPKwqHjubywKJl5v6TTtlEt5tzfhx6t6ns7lilH9uvdmEpCVfk8YS/PLUzkcHYej/Rvy4P921I1yJqqVTZW+I2pBPYfzeHpBYksTd5HVPO6zLwrjs7N6ng7lvESK/zGVGCqyuw1abywOJnc/EKeHNyRuy+OIMiaqlVqVviNqaB2Zjiaqv2QcpBeEQ2YOCyK1tZUzWCF35gKp6BQmbFiO68t2URggPDi9V0Y2SvMmqqZ06zwG1OBbNmXxRNz4/ll52Eu7xDKSzdE0axedW/HMj7GCr8xFUBufiFTvtvKW8tSqFk1kL/e3JWhXZtZUzVTJCv8xvi5+LTDPDEnno17s7g2phnPXduZkFrWVM0Uzwq/MX7qRG4Bf/16M9O+TyW0dlWm3R7LVZ0bezuW8QNW+I3xQz+mZjB2bjzbM7K5pVdLxg7uRN3q1lTNuMcKvzF+JCsnj4lfbOSjVTsJa1CDj++Jo09ba6pmSsYKvzF+YtnGfYybn8i+ozncc3EEfx7QgepVrN2CKTkr/Mb4uMzjuUz4LIkF63bTvnEt/nFrH7qFWVM1U3pW+I3xUarKZ/F7GL8wiaycPB69oh0PXt6WKkHWbsFcGCv8xvigvUccTdW+3rCPmBZ1mXRjHB2bWFM14xlW+I3xIarKrNW7eHnxBvIKC3n6mk7c2TeCQGu3YDzICr8xPmJHxnHGzk1gZWoGvVs3ZOLwKFo1rOntWKYCssJvjJcVFCof/G8bry3dRHBAAK8Mi2JEz5bWbsGUGSv8xnjRpr2Opmrrdx3myk6NePH6KJrUrebtWKaCc+vjASIyXUT2i0hiMdNFRP4uIikiEi8i3V2mjRKRLc6vUZ4Kbow/y80v5I2vNjPkze/ZlZnN32/pxrTbY63om3Lh7hH/DOAtYGYx0wcD7ZxfccA7QJyINACeA2IBBdaKyEJVPXQhoY3xZ+t2HeaJOevZvO8YQ7s247lrI2lQs4q3Y5lKxK3Cr6rLRST8HLMMBWaqqgI/ikg9EWkKXAZ8paqZACLyFTAI+ORCQhvjj07kFvCXpZuY/r9tNKpdjfdHxXJFJ2uqZsqfp87xNwd2uTxOc44rbvxviMhoYDRAWFiYh2IZ4xtWbD3I2LkJ7MzM5ta4MMYO7kjtatZUzXiHz1zcVdWpwFSA2NhY9XIcYzziaE4er3y+gU9+2kV4wxrMGn0RF7Vu6O1YppLzVOFPB1q6PG7hHJeO43SP6/hvPbRNY3za18n7GLcggQNZJ7mvX2v+eGV7a6pmfIKnCv9C4CERmYXj4u4RVd0jIkuAl0XkVEepAcCTHtqmMT7p4LGTPP9ZMp+t303HJrWZdnss0S3qeTuWMae5VfhF5BMcR+4hIpKG45M6wQCqOgX4HLgaSAGygTud0zJF5AVgtXNVE05d6DWmolFVPl23m+c/S+LYyXz+dFV77r+0jTVVMz7H3U/13HKe6Qo8WMy06cD0kkczxn/sPnyCpxcksmzjfrqF1WPS8GjaN67t7VjGFMlnLu4a448KC5WPf9rJxC82UlCoPDukM6P6hFtTNePTrPAbU0rbDh5n7Nx4Vm3LpG/bhrxyQzRhDWt4O5Yx52WF35gSyi8o5P0ftvH6V5upEhTA5OHR3BTbwpqqGb9hhd+YEkjefZQxc+NJSD/CVZ0b8+L1XWhcx/rrGP9ihd8YN5zML+CtZSm88+1W6tUI5u2R3bk6qokd5Ru/ZIXfmPNYu+MQY+bGk7L/GMO6N+eZazpT35qqGT9mhd+YYmTn5vPqkk3MWLGdpnWq8cGdPbm8QyNvxzLmglnhN6YIP2w5yNh58aQdOsHtvVvxxKCO1KpqbxdTMdgr2RgXR7LzeOnzZP6zJo2IkJr8577e9Ipo4O1YxniUFX5jnL5M3MsznyaSeTyXBy5rw6NXtKNasDVVMxWPFX5T6R3IOsn4hUksTthD56Z1+OCOnnRpXtfbsYwpM1b4TaWlqsz7OZ0Ji5I5kVvA4wM7MLpfa4IDramaqdis8JtKKf3wCZ6al8B3mw/Qo1V9Jg2Ppm2jWt6OZUy5sMJvKpXCQuVfq3Yw6YuNKDD+2s7c3jucAGuqZioRK/ym0th64Bhj58azevshLmkXwss3RNGygTVVM5WPFX5T4eUVFDLt+1T++vUWqgUF8OqN0dzYw5qqmcrLCr+p0BLTjzBmbjxJu48yuEsTnh8aSaPa1lTNVG5W+E2FlJNXwJvLtjDlu1Tq16jCO7d2Z3BUU2/HMsYnuHvP3UHA34BA4D1VnXjW9DeAy50PawCNVLWec1oBkOCctlNVr/NEcGOKs2Z7Jk/MjSf1wHFu7NGCp6/pRL0a1lTNmFPOW/hFJBB4G7gKSANWi8hCVU0+NY+qPuYy/8NAN5dVnFDVrp6LbEzRjp90NFX7cOV2mtWtzsy7etGvfai3Yxnjc9w54u8FpKhqKoCIzAKGAsnFzH8L8Jxn4hnjnu82H+CpeQnsPnKCUb3DeXxgB2paUzVjiuTOO6M5sMvlcRoQV9SMItIKiACWuYyuJiJrgHxgoqouKGVWY37jcHYuLyzawNyf02gTWpPZ9/UmNtyaqhlzLp4+JBoBzFHVApdxrVQ1XURaA8tEJEFVt569oIiMBkYDhIWFeTiWqYi+SNjDM58mcSg7l4cub8tD/dtaUzVj3OBO4U8HWro8buEcV5QRwIOuI1Q13fk9VUS+xXH+/zeFX1WnAlMBYmNj1Y1cppLafzSHZz9N4sukvUQ2q8OHd/Ukspk1VTPGXe4U/tVAOxGJwFHwRwAjz55JRDoC9YGVLuPqA9mqelJEQoC+wGRPBDeVj6oyZ20aLyxKJie/kDGDOnLvJREEWVM1Y0rkvIVfVfNF5CFgCY6Pc05X1SQRmQCsUdWFzllHALNU1fVovRPwrogUAgE4zvEXd1HYmGLtyszmqfkJfL/lID3D6zNxeDRtQq2pmjGlIWfWad8QGxura9as8XYM4wMKCpWZK7fz6pJNCDB2cEdujWtlTdWMOYuIrFXVWHfmtc+7GZ+Vsj+LMXMTWLvjEJe2D+XlYVE0r1fd27GM8XtW+I3PySso5N3vtvL3/6ZQo2ogr/8uhhu6NbemasZ4iBV+41MS04/w+Jx4Nuw5yjXRTRl/bSShtat6O5YxFYoVfuMTcvIK+OvXW5j2fSoNalbh3dt6MDCyibdjGVMhWeE3XrcqNYOx8xLYdvA4N8e25KmrO1G3RrC3YxlTYVnhN16TlZPH5C838c8fd9CifnX+dXccF7cL8XYsYyo8K/zGK77ZtJ9x8xLYczSHu/pG8H8D21Ojir0cjSkP9k4z5erQ8VxeWJTMvF/SadeoFnMf6EP3sPrejmVMpWKF35QLVWVxwh6e+zSJIyfyeKR/Wx7s35aqQdZUzZjyZoXflLl9R3N4ekEiXyXvI6p5Xf51TxydmtbxdixjKi0r/KbMqCr/WbOLFxdvIDe/kCcHd+Tui62pmjHeZoXflImdGdmMnRfPiq0Z9IpowKTh0USE1PR2LGMMVviNhxUUKjNWbOe1JZsIDBBeuqELt/QMs6ZqxvgQK/zGYzbvy+KJOfGs23WY/h0b8dINXWha15qqGeNrrPCbC5abX8iU77by5rIt1KoaxN9GdOW6mGbWVM0YH2WF31yQ9bsOM2ZuPBv3ZnFtTDPGX9uZhrWsqZoxvswKvymVE7kFvPH1Zt77PpXQ2lWZdnssV3Vu7O1Yxhg3WOE3JbZyawZPzotne0Y2t/QK48mrO1KnmjVVM8ZfWOE3bjuak8fELzby8aqdtGpYg4/vjaNPG2uqZoy/ces/aURkkIhsEpEUERlbxPQ7ROSAiKxzft3jMm2UiGxxfo3yZHhTfpZt3MeA15cz66ed3HtJBF8+2s+KvjF+6rxH/CISCLwNXAWkAatFZKGqJp81679V9aGzlm0APAfEAgqsdS57yCPpTZnLOHaSCYuS+XTdbjo0rs2U23rQtWU9b8cyxlwAd0719AJSVDUVQERmAUOBswt/UQYCX6lqpnPZr4BBwCeli2vKi6qycP1unv8smaycPP54ZTv+cFlbqgRZuwVj/J07hb85sMvlcRoQV8R8w0WkH7AZeExVdxWzbPOiNiIio4HRAGFhYW7EMmVlz5ETPD0/kf9u3E9My3pMHh5Nhya1vR3LGOMhnrq4+xnwiaqeFJH7gA+B/iVZgapOBaYCxMbGqodymRIoLFRmrd7FK59vIK+wkKev6cSdfSMItHYLxlQo7hT+dKCly+MWznGnqWqGy8P3gMkuy1521rLfljSkKXvbDx5n7Lx4fkzNpHfrhkwcHkWrhtZUzZiKyJ3CvxpoJyIROAr5CGCk6wwi0lRV9zgfXgdscA4vAV4WkVO3WBoAPHnBqY3HFBQq03/Yxl++2kRwQAATh0Vxc8+W1m7BmArsvIVfVfNF5CEcRTwQmK6qSSIyAVijqguBR0TkOiAfyATucC6bKSIv4PjlATDh1IVe430b9x5lzJx41qcd4cpOjXjx+iia1K3m7VjGmDImqr53Oj02NlbXrFnj7RgV1sn8At4WO1LNAAAOzUlEQVT+Ziv/+CaFutWDGX9dJEOim9pRvjF+TETWqmqsO/Paf+5WMr/sPMSYufFs3neM67s249lrI2lQs4q3YxljypEV/koiOzefvyzdzPT/baNJnWpMvyOW/h2tqZoxlZEV/kpgRcpBxs5LYGdmNr+/KIwxgzpS25qqGVNpWeGvwI6cyOOVzzcwa/UuwhvWYNboi7iodUNvxzLGeJkV/gpqadJenl6QyMFjJ7nv0tY8dmV7qgUHejuWMcYHWOGvYA4eO8n4hUksit9Dxya1eW9ULNEtrKmaMeZXVvgrCFVlwbp0nv8smeyTBfz5qvbcd2kba6pmjPkNK/wVwO7DJxg3P4FvNh2gW5ijqVq7xtZUzRhTNCv8fqywUPnop51M+mIjBYXKs0M6M6pPuDVVM8ackxV+P5V64Bhj5yXw07ZMLm4bwivDomjZoIa3Yxlj/IAVfj+TX1DIez9s442vNlMlKIDJw6O5KbaFtVswxrjNCr8fSd59lCfmricx/SgDOjfmheu70LiONVUzxpSMFX4/cDK/gLeWpfDOt1upVyOYf9zancFdmthRvjGmVKzw+7i1OxxN1VL2H2NY9+Y8c01n6ltTNWPMBbDC76OOn8zntaWbmLFiO83qVmfGnT25rEMjb8cyxlQAVvh90PdbDvDkvATSDp3g9t6teGJQR2pVtR+VMcYzrJr4kCPZeby4OJnZa9NoHVKT/9zXm14RDbwdyxhTwVjh9xFfJu7lmU8TyTyeywOXteHRK9pZUzVjTJlwq/CLyCDgbzjuufueqk48a/qfgHtw3HP3AHCXqu5wTisAEpyz7lTV6zyUvULYn5XD+IVJfJ6wl85N6/DBHT3p0ryut2MZYyqw8xZ+EQkE3gauAtKA1SKyUFWTXWb7BYhV1WwReQCYDNzsnHZCVbt6OLffU1Xm/ZzOhEXJnMgr4PGBHRjdrzXBgdZUzRhTttw54u8FpKhqKoCIzAKGAqcLv6p+4zL/j8DvPRmyokk7lM1T8xNZvvkAPVrVZ9LwaNo2quXtWMaYSsKdwt8c2OXyOA2IO8f8dwNfuDyuJiJrcJwGmqiqC4paSERGA6MBwsLC3IjlfwoLlX/+uINJX24E4PnrIrntolYEWFM1Y0w58ujFXRH5PRALXOoyupWqpotIa2CZiCSo6tazl1XVqcBUgNjYWPVkLl+w9cAxxsyJZ82OQ1zSLoSXb7CmasYY73Cn8KcDLV0et3COO4OIXAmMAy5V1ZOnxqtquvN7qoh8C3QDflP4K6q8gkKmLk/lb//dQvXgQF67KYbh3ZtbuwVjjNe4U/hXA+1EJAJHwR8BjHSdQUS6Ae8Cg1R1v8v4+kC2qp4UkRCgL44Lv5VCYvoRxsyNJ2n3Ua6OasL46yJpVNuaqhljvOu8hV9V80XkIWAJjo9zTlfVJBGZAKxR1YXAq0AtYLbzSPbUxzY7Ae+KSCEQgOMcf3KRG6pAcvIK+Pt/t/Du8lTq16jClN93Z1CXpt6OZYwxAIiq751Oj42N1TVr1ng7Rqms3p7JmLnxpB44zk09WvD0NZ2pWyPY27GMMRWciKxV1Vh35rX/3PWQYyfzmfzlRmau3EHzetWZeVcv+rUP9XYsY4z5DSv8HvDd5gM8NS+B3UdOcEefcB4f2IGa1lTNGOOjrDpdgMPZuUxYlMy8n9NpE1qT2ff1JjbcmqoZY3ybFf5S+jxhD89+msjh7DweurwtD/Vva03VjDF+wQp/Ce0/msMznyayJGkfXZrX4cO7ehHZzJqqGWP8hxV+N6kqs9em8eKiZHLyCxkzqCP3XhJBkDVVM8b4GSv8btiVmc2T8xL4IeUgvcIbMHF4FK1DramaMcY/WeE/h4JCZebK7Uz+chMBAi8MjeTWOGuqZozxb1b4i5GyP4sn5sTz887DXNYhlJduiKJ5verejmWMMRfMCv9Z8goKefe7rfz9vynUqBrIGzfHcH1Xa6pmjKk4rPC7SEg7wuNz1rNxbxbXRDfl+esiCalV1duxjDHGo6zw42iq9sbXm5m2PJWQWlV597YeDIxs4u1YxhhTJip94V+VmsHYeQlsO3icm2Nb8tQ1nahb3ZqqGWMqrkpb+LNy8pj05Ub+9eNOWjaozkf3xNG3bYi3YxljTJmrlIX/m437GTc/gT1Hc7j74gj+PKA9NapUyl1hjKmEKlW1yzyeywuLkpn/SzrtGtVi7gN96B5W39uxjDGmXFWKwq+qLIrfw/iFSRw5kccjV7TjwcvbUDXImqoZYyqfCl/49x3NYdz8RL7esI/oFnX51z1xdGpax9uxjDHGayps4VdV/r16Fy99voHc/EKeurojd/W1pmrGGONW4ReRQcDfcNxs/T1VnXjW9KrATKAHkAHcrKrbndOeBO4GCoBHVHWJx9IXY2dGNmPnxbNiawZxEQ2YNDya8JCaZb1ZYzxLFbTQ5Xsh4DJ8xniKGV/U/HqO9biOx43tumb0RM7SPCf1cEbX+TyxL0uw3er1YdRnZf7SOm/hF5FA4G3gKiANWC0iC1U12WW2u4FDqtpWREYAk4CbRaQzMAKIBJoBX4tIe1Ut8PQTAUdTtQ/+t43Xlm4iKCCAl27owi09wzzfVK1MXkBl8UIvyZvXmznL4s1bhhkvqAiU4GdiSklAApxfLsNnjKeY8WfPL+dZj8s859tuQCBI0Lm3W61euewhd474ewEpqpoKICKzgKGAa+EfCox3Ds8B3hJHc5uhwCxVPQlsE5EU5/pWeib+r45k53Hw9Yu4LDebIVUDCKlVhaAfFVaWpLC5WQRMKXnwTeQ6b7HzixvrOfWGdGe77uQvwXY9/nydpzHd2m5JcuLGdkuQs0T7sjTPyZyPO4W/ObDL5XEaEFfcPKqaLyJHgIbO8T+etWzzojYiIqOB0QBhYWHuZD9DnepBbKkRTtPGQTSuVwMptzcvbr4pyuKF7m7R8HZOe0Ma40t85uKuqk4FpgLExsaW+O9cESH2T3M9nssYYyqaADfmSQdaujxu4RxX5DwiEgTUxXGR151ljTHGlCN3Cv9qoJ2IRIhIFRwXaxeeNc9CYJRz+EZgmaqqc/wIEakqIhFAO+Anz0Q3xhhTGuc91eM8Z/8QsATHxzmnq2qSiEwA1qjqQuB94J/Oi7eZOH454JzvPzguBOcDD5bVJ3qMMca4RxwH5r4lNjZW16xZ4+0YxhjjN0RkrarGujOvO6d6jDHGVCBW+I0xppKxwm+MMZWMFX5jjKlkfPLirogcAHaUcvEQ4KAH45QVy+l5/pLVcnqWv+SEss3aSlVD3ZnRJwv/hRCRNe5e2fYmy+l5/pLVcnqWv+QE38lqp3qMMaaSscJvjDGVTEUs/FO9HcBNltPz/CWr5fQsf8kJPpK1wp3jN8YYc24V8YjfGGPMOVjhN8aYSqbCFH4RGSQim0QkRUTGeinDdhFJEJF1IrLGOa6BiHwlIluc3+s7x4uI/N2ZN15EurusZ5Rz/i0iMqq47ZUw23QR2S8iiS7jPJZNRHo4n3uKc9lS3XKrmJzjRSTduV/XicjVLtOedG5zk4gMdBlf5OvB2V58lXP8v52txkuTs6WIfCMiySKSJCKPOsf71D49R05f3KfVROQnEVnvzPr8udYvjnbv/3aOXyUi4aV9Dh7KOUNEtrns067O8V57PxVLVf3+C0e76K1Aa6AKsB7o7IUc24GQs8ZNBsY6h8cCk5zDVwNf4Lgp4kXAKuf4BkCq83t953B9D2TrB3QHEssiG477LFzkXOYLYLAHc44H/q+IeTs7f9ZVgQjnayDwXK8H4D/ACOfwFOCBUuZsCnR3DtcGNjvz+NQ+PUdOX9ynAtRyDgcDq5zPv8j1A38ApjiHRwD/Lu1z8FDOGcCNRczvtfdTcV8V5Yj/9A3hVTUXOHVDeF8wFPjQOfwhcL3L+Jnq8CNQT0SaAgOBr1Q1U1UPAV8Bgy40hKoux3GvBI9nc06ro6o/quNVO9NlXZ7IWZyhwCxVPamq24AUHK+FIl8PzqOm/sCcIp5zSXPuUdWfncNZwAYc95P2qX16jpzF8eY+VVU95nwY7PzSc6zfdV/PAa5w5inRc/BgzuJ47f1UnIpS+Iu6Ify5XtxlRYGlIrJWHDePB2isqnucw3uBxs7h4jKX53PxVLbmzuGzx3vSQ84/k6efOn1SipwNgcOqmu/JnM5TDN1wHPn57D49Kyf44D4VkUARWQfsx1EIt55j/aczOacfceYp8/fW2TlV9dQ+fcm5T98Qkapn53QzT5m/nypK4fcVF6tqd2Aw8KCI9HOd6Pzt7ZOfn/XlbMA7QBugK7AH+It34/xKRGoBc4E/qupR12m+tE+LyOmT+1RVC1S1K477c/cCOno5UpHOzikiXYAnceTtieP0zRgvRjynilL4feKm7qqa7vy+H5iP44W7z/mnG87v+52zF5e5PJ+Lp7KlO4fLJLOq7nO+0QqBaTj2a2lyZuD4MzvorPGlIiLBOIrpR6o6zzna5/ZpUTl9dZ+eoqqHgW+A3udY/+lMzul1nXnK7b3lknOQ87SaqupJ4ANKv0/L9P10Krjff+G4d3Aqjgs5py7aRJZzhppAbZfhFTjOzb/KmRf7JjuHr+HMCz4/6a8XfLbhuNhT3zncwEMZwznzoqnHsvHbi1FXezBnU5fhx3CcvwWI5MyLeKk4LuAV+3oAZnPmhcI/lDKj4Dj3+tezxvvUPj1HTl/cp6FAPedwdeB7YEhx6wce5MyLu/8p7XPwUM6mLvv8r8BEX3g/FfkcPLkyb37huHK+Gcc5wXFe2H5r5wtpPZB0KgOOc47/BbYAX7v8YAV425k3AYh1WdddOC5IpQB3eijfJzj+pM/Dcc7wbk9mA2KBROcyb+H8r3AP5fynM0c8sJAzi9Y45zY34fLJh+JeD86f00/O/LOBqqXMeTGO0zjxwDrn19W+tk/PkdMX92k08IszUyLw7LnWD1RzPk5xTm9d2ufgoZzLnPs0EfgXv37yx2vvp+K+rGWDMcZUMhXlHL8xxhg3WeE3xphKxgq/McZUMlb4jTGmkrHCb4wxlYwVfmOMqWSs8BtjTCXz/y+93WiMzTb1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py35-cobra-cplex",
   "language": "python",
   "name": "py35-cobra-cplex"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}