market_simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Market Simulation\n",
    "\n",
    "This notebook creates one file.\n",
    "For a given number of market situations, each market situation is simulated with the same starting condition multiple times for a specified time horizon.\n",
    "In the `demand_learning_data` file we save explanatory variables for market situations and whether there was a purchase in the current situation. We also create three graphs displaying a selected time period from our simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** We create the output files for three different pricing setups:**\n",
    "* _Setup 1 (S1):_  For the watched company we chose a random price between 0 and 15. For all the other companies, we chose a price between 5 and 15.\n",
    "* _Setup 2 (S2):_ For each company we use the same pricing strategy:  \n",
    "```min_competitor_price = min(prices_companies)\n",
    "if min_competitor_price < 5:\n",
    "    10\n",
    "else:\n",
    "    max(4.9, (min_competitor_price - 0.5))\n",
    "```\n",
    "* _Setup 3 (S3):_ For the first three companies we use the approach from Setup 2, except with different thresholds. The last two companies have a fixed price of either 11 or 13."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import random\n",
    "import numpy as np\n",
    "import csv\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set start values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_firm = 5                                                                # the number of firms which we simulate\n",
    "time_horizon = 100                                                          \n",
    "num_market_events = time_horizon + 100                                      # the maximum number of market events\n",
    "mean_betw_customer_arrival = 1\n",
    "num_customer_events = int(time_horizon / mean_betw_customer_arrival + 100)  # the maximum number of customer events\n",
    "iter_diff_start = 100                                                       # how often to simulate with same start\n",
    "iter_same_start = 1                                                         # how often to simulate with \n",
    "                                                                            # different start\n",
    "watched_firm = 0\n",
    "current_line = 0\n",
    "setups = [\"S1\",\"S2\",\"S3\"]                                                   # set the names of different setups\n",
    "random.seed(3819469)                                                        # set a random seed to \n",
    "                                                                            # always get comparable results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set market start situation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "prices_firm = np.round(np.random.uniform(6,16, size = (iter_diff_start, num_firm)), decimals = 2)\n",
    "qual_firm = np.round(np.random.uniform(0.5,5.5, size = (iter_diff_start, num_firm)))\n",
    "rating_firm = np.round(np.random.uniform(90, 100, size = (iter_diff_start, num_firm)), decimals = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate all the necessary information for when firms adjusts their prices and how often they adjust their prices.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "time_betw_market_events_firm = np.random.uniform(0.8, 1.2, size = (iter_diff_start, iter_same_start, num_firm, num_market_events))\n",
    "\n",
    "time_adjusts_firm = time_betw_market_events_firm\n",
    "\n",
    "time_adjusts_firm[:,:,:,0] = np.random.uniform(0, 1, size = (iter_diff_start, iter_same_start, num_firm))\n",
    "\n",
    "for d in range(iter_diff_start):\n",
    "    for s in range(iter_same_start):\n",
    "        for c in range(num_firm):\n",
    "            for i in range(1, num_market_events):\n",
    "                time_adjusts_firm[d,s,c,i] = time_adjusts_firm[d,s,c,i-1] + time_betw_market_events_firm[d,s,c,i]\n",
    "                \n",
    "price_adjusts_before_horizon = np.array(\n",
    "    [[\n",
    "            set([\n",
    "                time_adjusts_firm[d,s].flatten()[i] \n",
    "                for i in range(len(time_adjusts_firm[d,s].flatten())) \n",
    "                if time_adjusts_firm[d,s].flatten()[i] < time_horizon])\n",
    "            for s in range(iter_same_start)]\n",
    "        for d in range(iter_diff_start)])\n",
    "        \n",
    "price_adjusts_before_horizon = np.array(\n",
    "    [[\n",
    "            [0] + sorted(price_adjusts_before_horizon[d,s]) + [time_horizon]\n",
    "            for s in range(iter_same_start)]\n",
    "        for d in range(iter_diff_start)])\n",
    "\n",
    "num_price_adjusts_before_horizon = np.array(\n",
    "    [[\n",
    "            len(price_adjusts_before_horizon[d,s])\n",
    "            for s in range(iter_same_start)]\n",
    "        for d in range(iter_diff_start)])\n",
    "\n",
    "firm_adjusting_at_time = np.array(\n",
    "    [[[\n",
    "                min(np.where(time_adjusts_firm[d,s]==price_adjusts_before_horizon[d,s][i])[0])\n",
    "                if price_adjusts_before_horizon[d,s][i]!=time_horizon and i!=0 \n",
    "                else -1  \n",
    "                for i in range(num_price_adjusts_before_horizon[d,s])]\n",
    "            for s in range(iter_same_start)]\n",
    "        for d in range(iter_diff_start)])    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate when customers will arrive and how they will score the quality and rating of all firms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "time_betw_customer_events = (mean_betw_customer_arrival * \n",
    "                             np.random.exponential(size = (iter_diff_start, iter_same_start, num_customer_events)))\n",
    "\n",
    "time_customer_events = time_betw_customer_events\n",
    "\n",
    "for d in range(iter_diff_start):\n",
    "    for s in range(iter_same_start):\n",
    "        for k in range(1, num_customer_events):\n",
    "            time_customer_events[d,s,k] = time_customer_events[d,s,k-1] + time_betw_customer_events[d,s,k]\n",
    "            \n",
    "rand_customer_score_qual = np.random.uniform(0, 1, size = (iter_diff_start, iter_same_start, num_customer_events))\n",
    "\n",
    "rand_customer_score_rating = np.random.uniform(0, 0.5, size = (iter_diff_start, iter_same_start, num_customer_events))\n",
    "\n",
    "customer_score_qual_rating = np.array(\n",
    "    [[[[\n",
    "                    (rand_customer_score_qual[d,s,c]*qual_firm[d,k] +\n",
    "                     rand_customer_score_rating[d,s,c]*(100-rating_firm[d,k]))\n",
    "                    for k in range(num_firm)] \n",
    "                for c in range(num_customer_events)]\n",
    "            for s in range(iter_same_start)]\n",
    "        for d in range(iter_diff_start)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create variables for watched firm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "count_adjusts_watched_firm = np.array(\n",
    "    [[\n",
    "            sum([\n",
    "                1 \n",
    "                for i in range(num_price_adjusts_before_horizon[d,s])\n",
    "                if firm_adjusting_at_time[d,s][i] == watched_firm]) + 1 \n",
    "            for s in range(iter_same_start)]\n",
    "        for d in range(iter_diff_start)])\n",
    "\n",
    "time_adjusts_watched_firm = np.array(\n",
    "    [[[\n",
    "                price_adjusts_before_horizon[d,s][i] \n",
    "                for i in range(num_price_adjusts_before_horizon[d,s]) \n",
    "                if firm_adjusting_at_time[d,s][i] == watched_firm]\n",
    "            for s in range(iter_same_start)] \n",
    "        for d in range(iter_diff_start)])\n",
    "\n",
    "for d in range(iter_diff_start):\n",
    "    for s in range(iter_same_start):\n",
    "        time_adjusts_watched_firm[d,s].insert(0,0)\n",
    "        time_adjusts_watched_firm[d,s].extend([time_horizon])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def price_if_smaller_else(price, price_threshold, min_price, prices_firm, c):\n",
    "    return(\n",
    "        price \n",
    "        if min([prices_firm[c] for firm in range(num_firm) if firm!=c]) < price_threshold \n",
    "        else round(\n",
    "            min([\n",
    "                    max(min_price, prices_firm[firm] - 0.5) for firm in range(num_firm) if firm != c])\n",
    "            ,2))\n",
    "s1 = {\n",
    "        '0' : lambda prices_firm, c: round(random.uniform(0,15),2),\n",
    "        '1' : lambda prices_firm, c: round(random.uniform(5,15),2),\n",
    "        '2' : lambda prices_firm, c: round(random.uniform(5,15),2),\n",
    "        '3' : lambda prices_firm, c: round(random.uniform(5,15),2),\n",
    "        '4' : lambda prices_firm, c: round(random.uniform(5,15),2)\n",
    "    }\n",
    "\n",
    "s2 = {\n",
    "        '0' : lambda prices_firm, c: price_if_smaller_else(10,5,4.9,prices_firm,c),\n",
    "        '1' : lambda prices_firm, c: price_if_smaller_else(10,5,4.9,prices_firm,c),\n",
    "        '2' : lambda prices_firm, c: price_if_smaller_else(10,5,4.9,prices_firm,c),\n",
    "        '3' : lambda prices_firm, c: price_if_smaller_else(10,5,4.9,prices_firm,c),\n",
    "        '4' : lambda prices_firm, c: price_if_smaller_else(10,5,4.9,prices_firm,c)\n",
    "    }\n",
    "\n",
    "s3 = {\n",
    "        '0' : lambda prices_firm, c: price_if_smaller_else(10,5,4.9,prices_firm,c),\n",
    "        '1' : lambda prices_firm, c: price_if_smaller_else(9,4,3.9,prices_firm,c),\n",
    "        '2' : lambda prices_firm, c: price_if_smaller_else(12,6,5.9,prices_firm,c),\n",
    "        '3' : lambda prices_firm, c: 11,\n",
    "        '4' : lambda prices_firm, c: 13\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def explanatory_vars(d, s, explanatory_var, count_adjust, setup):\n",
    "    return {\n",
    "        '0' : 1,\n",
    "        '1' : price_ranks_firm[setup,d,s][count_adjust],\n",
    "        '2' : 1 if price_ranks_firm[setup,d,s][count_adjust]==1 else 0,\n",
    "        '3' : qual_rank_firm[d],\n",
    "        '4' : rating_ranks[d,watched_firm],\n",
    "        '5' : prices_tracked[setup,d,s][count_adjust][watched_firm],\n",
    "        '6' : prices_tracked[setup,d,s][count_adjust][watched_firm] - min([prices_tracked[setup,d,s][count_adjust][j] for j in range(num_firm) if j!=watched_firm]) if num_firm > 1 else 0,\n",
    "        '7' : qual_firm[d,watched_firm],\n",
    "        '8' : rating_firm[d,watched_firm],\n",
    "        '9' : 15 - (prices_tracked[setup,d,s][count_adjust][watched_firm] + \n",
    "                    0.5*qual_firm[d,watched_firm] + \n",
    "                    0.25*(100-rating_firm[d,watched_firm])) \n",
    "        if((prices_tracked[setup,d,s][count_adjust][watched_firm] +\n",
    "            0.5*qual_firm[d,watched_firm] +\n",
    "            0.25*(100-rating_firm[d,watched_firm])) \n",
    "           < min([\n",
    "               prices_tracked[setup,d,s][count_adjust][j] +\n",
    "                  0.5*qual_firm[d,j] +\n",
    "                  0.25*(100-rating_firm[d,j]) \n",
    "               for j in range(num_firm) \n",
    "               if j!= watched_firm]))\n",
    "        else 0,\n",
    "    }[str(explanatory_var)]\n",
    "\n",
    "def calc_prices(prices_firm_adjust, adjusting_firm, setup):\n",
    "    prices_firm_adjust = np.array(prices_firm_adjust)\n",
    "    for comp in range(num_firm):\n",
    "        prices_firm_adjust[0:, comp] = setup[str(comp)](prices_firm_adjust[0],comp)\n",
    "    #weird indexing because there are actually no price adjustments at time 0 nor on the price horizon\n",
    "    for i in range(1,len(prices_firm_adjust)-1):        \n",
    "        prices_firm_adjust[i:,adjusting_firm[i]] = setup[str(adjusting_firm[i])](prices_firm_adjust[i-1],adjusting_firm[i])\n",
    "    return prices_firm_adjust"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we don't want to know what happens after or on the time horizon, we will always substract 1 when using the value num_price_adjusts_before_horizon in the following code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# get all the customers that appear between market events, so we don't always have to iterate through each customer\n",
    "gen = np.array(\n",
    "    [[[[\n",
    "                    x \n",
    "                    for x in range(num_customer_events) \n",
    "                    if price_adjusts_before_horizon[d,s][i]< time_customer_events[d,s][x]<= price_adjusts_before_horizon[d,s][i+1]] \n",
    "                for i in range(num_price_adjusts_before_horizon[d,s] - 1)] \n",
    "            for s in range(iter_same_start)] \n",
    "        for d in range(iter_diff_start)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the prices adjustments for each firm in the time horizon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prices_firm_adjust = np.array(\n",
    "    [[[[\n",
    "                    prices_firm[d,k]  \n",
    "                    for k in range(num_firm)]\n",
    "                for _ in range(num_price_adjusts_before_horizon[d,s] - 1)]\n",
    "            for s in range(iter_same_start)] \n",
    "        for d in range(iter_diff_start)])\n",
    "\n",
    "prices_firm_adjust = np.array(\n",
    "    [[[\n",
    "                np.array(calc_prices(prices_firm_adjust[d,s], firm_adjusting_at_time[d,s], globals()[\"s{}\".format(setup+1)])) \n",
    "                for s in range(iter_same_start)] \n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate the customer by completing the score using the current prices of each firm and then calculating a decision for each customer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "customer_score = np.array(\n",
    "    [[[[[[\n",
    "                            customer_score_qual_rating[d,s,c,k]+prices_firm_adjust[setup,d,s][i,k] \n",
    "                            for k in range(num_firm)] \n",
    "                        for c in gen[d,s][i]]\n",
    "                    for i in range(num_price_adjusts_before_horizon[d,s] - 1)]\n",
    "                for s in range(iter_same_start)]\n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])\n",
    "\n",
    "customer_decision = np.array(\n",
    "    [[[[[\n",
    "                        np.argmin(customer_score[setup,d,s][i][c]) \n",
    "                        if min(customer_score[setup,d,s][i][c])< random.uniform(5,15) \n",
    "                        else -1 \n",
    "                        for c in range(len(gen[d,s][i]))] \n",
    "                    for i in range(num_price_adjusts_before_horizon[d,s] - 1)]\n",
    "                for s in range(iter_same_start)]\n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Track data that is used to create further graphs and save them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "accumulated_profit = np.array(\n",
    "    [[[[\n",
    "                    sum([\n",
    "                        sum([\n",
    "                            prices_firm_adjust[setup,d,s][i,k] \n",
    "                            if customer_decision[setup,d,s][i][c]==k \n",
    "                            else 0 \n",
    "                            for c in range(len(gen[d,s][i]))]) \n",
    "                        for i in range(num_price_adjusts_before_horizon[d,s] - 1)])\n",
    "                    for k in range(num_firm)]\n",
    "                for s in range(iter_same_start)]\n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])\n",
    "\n",
    "prices_tracked = np.array(\n",
    "    [[[[[\n",
    "                        prices_firm_adjust[setup,d,s][i,k] \n",
    "                        for k in range(num_firm) \n",
    "                        if firm_adjusting_at_time[d,s][i]==watched_firm or price_adjusts_before_horizon[d,s][i]==0] \n",
    "                    for i in range(num_price_adjusts_before_horizon[d,s] - 1)]\n",
    "                for s in range(iter_same_start)]\n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])\n",
    "\n",
    "prices_tracked = np.array(\n",
    "    [[[[\n",
    "                    x \n",
    "                    for x in prices_tracked[setup,d,s] \n",
    "                    if x != []] \n",
    "                for s in range(iter_same_start)]\n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])\n",
    "\n",
    "index_time_adjusts = np.array(\n",
    "    [[[[[\n",
    "                        x \n",
    "                        for x in range(num_price_adjusts_before_horizon[d,s]) \n",
    "                        if price_adjusts_before_horizon[d,s][x]==time_adjusts_watched_firm[d,s][i]] \n",
    "                    for i in range(count_adjusts_watched_firm[d,s] + 1)]\n",
    "                for s in range(iter_same_start)]\n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])\n",
    "\n",
    "sales_betw_adjusts = np.array(\n",
    "    [[[[\n",
    "                    sum([\n",
    "                        sum([\n",
    "                            1 \n",
    "                            for c in range(len(gen[d,s][j])) \n",
    "                            if customer_decision[setup,d,s][j][c]==watched_firm])\n",
    "                        for j in range(num_price_adjusts_before_horizon[d,s]) \n",
    "                        if j>=index_time_adjusts[setup,d,s][i][0] and j<index_time_adjusts[setup,d,s][i+1][0]]) \n",
    "                    for i in range(count_adjusts_watched_firm[d,s])]\n",
    "                for s in range(iter_same_start)]\n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])\n",
    "\n",
    "price_ranks_firm = np.array(\n",
    "    [[[[\n",
    "                  1 +  \n",
    "                  sum([\n",
    "                      1 \n",
    "                      for i in range(num_firm) \n",
    "                      if prices_tracked[setup,d,s][j][watched_firm]> prices_tracked[setup,d,s][j][i] and i!=watched_firm]) \n",
    "                  +\n",
    "                  sum([\n",
    "                      1 \n",
    "                      for i in range(num_firm) \n",
    "                      if prices_tracked[setup,d,s][j][watched_firm]== prices_tracked[setup,d,s][j][i] and i!=watched_firm])/2\n",
    "                  for j in range(count_adjusts_watched_firm[d,s])]\n",
    "                for s in range(iter_same_start)] \n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])\n",
    "\n",
    "qual_rank_firm = np.array(\n",
    "    [\n",
    "          1 +  \n",
    "          sum([\n",
    "              1 \n",
    "              for i in range(num_firm) \n",
    "              if qual_firm[d,watched_firm]> qual_firm[d,i] and i!=watched_firm]) \n",
    "          +\n",
    "          sum([\n",
    "              1 \n",
    "              for i in range(num_firm) \n",
    "              if qual_firm[d,watched_firm]== qual_firm[d,i] and i!=watched_firm])/2\n",
    "          for d in range(iter_diff_start)])\n",
    "\n",
    "rating_ranks = np.array(\n",
    "    [[\n",
    "          1 +  \n",
    "          sum([\n",
    "              1 \n",
    "              for i in range(num_firm) \n",
    "              if rating_firm[d,j]> rating_firm[d,i] and i!=j]) \n",
    "          +\n",
    "          sum([\n",
    "              1 \n",
    "              for i in range(num_firm) \n",
    "              if rating_firm[d,j]== rating_firm[d,i] and i!=j])/2\n",
    "          for j in range(num_firm)] \n",
    "        for d in range(iter_diff_start)])\n",
    "\n",
    "for setup in range(len(setups)):\n",
    "    current_line = 0\n",
    "    with open('demand_learning_data_{}.csv'.format(setups[setup]), 'w') as file:\n",
    "        writer = csv.writer(file, delimiter='\\t')\n",
    "        for d in range(iter_diff_start):\n",
    "            for s in range(iter_same_start):\n",
    "                for w in range(count_adjusts_watched_firm[d,s]):\n",
    "                    current_line = current_line +1\n",
    "                    line = list([current_line])\n",
    "                    for k in range(1,10):\n",
    "                        line.append(explanatory_vars(d, s, k, w, setup))\n",
    "                    line.append(sales_betw_adjusts[setup,d,s][w])\n",
    "                    writer.writerow(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Track data for graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "arrival_out = np.array(\n",
    "    [[[\n",
    "                np.array([\n",
    "                    (time_customer_events[d,s,x],  x) \n",
    "                    for x in gen[d,s][i]]) \n",
    "                for i in range(num_price_adjusts_before_horizon[d,s] - 1)]\n",
    "            for s in range(iter_same_start)] \n",
    "        for d in range(iter_diff_start)])\n",
    "\n",
    "buying_out = np.array(\n",
    "    [[[[[\n",
    "                        [time_customer_events[d,s][gen[d,s][i][c]], prices_firm_adjust[setup,d,s][i,customer_decision[setup,d,s][i][c]], gen[d,s][i][c], customer_decision[setup,d,s][i][c]]\n",
    "                        for c in range(len(gen[d,s][i])) \n",
    "                        if customer_decision[setup,d,s][i][c]!=-1] \n",
    "                    for i in range(num_price_adjusts_before_horizon[d,s] - 1)] \n",
    "                for s in range(iter_same_start)]\n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])\n",
    "\n",
    "traject_out = np.array(\n",
    "    [[[[[\n",
    "                        [prices_firm_adjust[setup,d,s][i,k], price_adjusts_before_horizon[d,s][i] ,k]\n",
    "                        for k in range(num_firm)] \n",
    "                    for i in range(num_price_adjusts_before_horizon[d,s] - 1)] \n",
    "                for s in range(iter_same_start)]\n",
    "            for d in range(iter_diff_start)]\n",
    "        for setup in range(len(setups))])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Illusetupes examples of price trajectories over time as well as the arrival of interested (indicated by vertical bars on the x-axis) and buying customers (indicated by dots) in case of K=5.\n",
    "Note, due to different qualities and ratings of competitors customers do not always choose the cheapest competitor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdcAAAEWCAYAAAA0MN3QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3X98VPWd7/HXJ/yohFhAiKACSSkV\na7FOS9butl6dh9pe2sV2733UXgV3q6DpfbTbNuptt653F/Z2ubXW2rRb3TVSiHYTa9f+sGut1d0a\ntPZiC3YQVkBZCwgaiVVUGgViPvePcwaGkB8zk5M5ZzLv5+PBY3J+zDmfjMhnvt/z/X4/5u6IiIhI\ndKriDkBERGS0UXIVERGJmJKriIhIxJRcRUREIqbkKiIiEjElVxERkYgpuYqIiERMyVUqkpmdbWa/\nMrNXzOwlM3vUzP4oz/e6mc0dwdhazGybmfWa2WUjdR8RGTlKrlJxzOytwL3APwAnAKcAfwcciDOu\nHBuBTwOPxx2IiBRHyVUq0akA7n6nu7/p7q+7+wPu/kT2BDNbamZbzOxlM/u5mdWF+x8OT9loZvvN\n7H+Y2WVm9svcG+S2bs2s1cz+ycweNLPXzGxt9nr9cfeb3f3fgTei/sVFpDSUXKUSPQW8aWa3m9mH\nzWxK7kEz+zPgr4H/DtQCjwB3Arj7OeFpZ7p7jbvflec9lwBfBqYBGaBt+L+GiCSVkqtUHHd/FTgb\ncOA2oMvMfmJm08NTPgV8xd23uHsP8H+B1GCtzTz81N0fdvcDwHXAn5jZrGFcT0QSTMlVKlKYOC9z\n95nAfOBkoDk8XAd808z2mdk+4CXACJ7NFuvZnHvvD6958jCuJyIJpuQqFc/dtwKtBEkWgkT4KXef\nnPNngrv/aoBL/AGozm6Y2Yx+zpmVc7yGYCDVc1HELyLJo+QqFcfMTjOza8xsZrg9C7gEWBee8k/A\ntWb2rvD4JDO7KOcSLwBzcrY3Au8ys5SZHQes6Oe2Hwmn/4wnePb6mLs/2895mNn48DoGjDOz48xM\n/6+KlBH9DyuV6DXgfcBjZvYHgqS6GbgGwN1/BHwV+J6ZvRoe+3DO+1cAt4fdxp9w96eA/wP8G/A0\ncNTI4VA7sJygO3gBwQCngTwAvA68H2gJfz5nkPNFJGFMxdJFRpaZtQK73f1/xx2LiJSGWq4iIiIR\nU3IVERGJmLqFRUREIqaWq4iISMTGxh1AMaZNm+b19fVxhyEiUlY2bNjworvXDuP9J44dO3YVwZzw\nSm6c9QKbe3p6rliwYMHe/k4oy+RaX1/P+vXr4w5DRKSsmNnO4bx/7Nixq2bMmPHO2tral6uqqir2\nmWJvb691dXWd3tnZuQr4aH/nVPI3DxERKcz82traVys5sQJUVVV5bW3tKxxZ1e3Yc0oYj4iIlLeq\nSk+sWeHnMGAOVXIVERGJmJKriIiUjTFjxiw47bTTTs/+2bZt2/iHH364+rLLLoushOPPfvazmtNP\nP/2dY8eOXbBmzZopQ7/jWGU5oElERCrTW97ylt6tW7c+mbtv3rx5B88555zuvuceOnSIcePGFXyP\nOXPmHFyzZs2O66+/fvrQZ/dPyVVERMravffee/zXv/716Q899ND2q6+++uTnn39+3K5du8afcMIJ\nPR/84Adf/clPfjK5t7fXtm3bNuEzn/lM58GDB6vuuuuuqePHj+994IEHnp4+ffqbudebN2/eQYCq\nquI7d5VcRUSkYEuXMmvz5iN1jKMwfz7dq1fTbynGrAMHDlSddtpppwPMmjXrwIMPPviffc954okn\nqh977LGtNTU1/q1vfWvqU089NWHjxo1Pvv7661Xz5s2b/zd/8zd7tmzZ8uSyZctm3XrrrVP/9m//\ntt+5qsOh5CoikoeWFmhvjzeGxYuhsTHeGOLWX7dwXwsXLtxXU1NzeFTz+9///temTJnSO2XKlN6a\nmpo3L7roon0AZ5xxRvcTTzwR6ReELCVXEZE8tLfD2rVw7rnx3H/t2uA1Kcl1qBZmnCZOnNibuz1+\n/PjDibaqqorjjjvOsz/39PTYSMSg5Coikqdzz4WOjnjunU7Hc18pjqbiiIiI5Fi7dm319OnT333f\nffdNueqqq+rmzp37rkKvUZYl5xoaGlxrC4tIKWVbjnG3XIdzfzPb4O4Nxb5/48aNO84888wXi49g\ndNm4ceO0M888s76/Y2q5ioiIREzJVUREJGJKriIiIhFTchUREYlYyabimNlqYBGw193n9zn2v4Cv\nAbXuroflIgmRhIUTQIsnSPkpZcu1FVjYd6eZzQI+COwqYSwikof2dshk4o0hk0lGghcpRMlaru7+\nsJnV93PoG8AXgXtKFYuI5C+Vim/6CWjxBDnamDFjFrzjHe94Pbt9zz33bH/hhRfGrl69empra2sk\nq0atWLFi+ne/+91pY8aM8alTp/bcfvvtO0499dSDhVwj1hWazOyjwB5332g2+ApUZtYINALMnj27\nBNGJiEjSlKLk3IIFC7qvueaaLccff3zvV7/61dqrrrpq5k9/+tNnCrlGbMnVzKqB64AP5XO+u7cA\nLRAsIjGCoYmISBmJuuTchRde+Fr257PPPnv/XXfdNbXQmOJsub4deBuQbbXOBB43s7PcvTPGuERE\nZAhL71k6a/PezdGWnDtxfvfqj61OVMm5W2+9tfaCCy54pdDfJbbk6u6bgBOz22a2A2jQaGERERlI\nKUvO3XLLLSds3Lix+tZbb91WaJylnIpzJ5AGppnZbmC5u3+nVPcXkRhEMJenOTtaOT3saI6lOT5F\nG6qFGacoSs79+Mc/Pv7GG2886ZFHHtk2YcKEgh9FlnK08CVDHK8vUSgiidWyoYX2TcmZd5JJBa/p\n1qP3Lz5jMY0L8khK2bk8qVTksQ1bdo6Rkqv08eijj0747Gc/W3ffffc9fcopp/QUcw3VcxVJkPZN\n7WQ6M6RmJDAZhTKdQVLKK7nCsOfyNKWD18inA2mOjwzgC1/4wqzu7u4xF1100dsBTj755IO/+MUv\nthdyDSVXkYRJzUjRcVlH3GEAOWXOmnP2tabjCEUEgO7u7t/23bdo0aLXFi1a9BrATTfd9Fzusc99\n7nO/B36f3d6zZ8+mgY5l/epXv3pquHFqbWEREZGIKbmKiIhETMlVREQkYkquIiIiEVNyFRERiZiS\nq4hImfjLv2wCmuIOQ/KgqTgRStoCAFDAZH8RSby5c2MurpsApSg5d8MNN9SuWrWqtqqqiokTJ77Z\n0tKyc8GCBW8Ucg0l1wglbQGAgif7i4gkXClKzl1xxRW//+IXv9gF0NbWNqmpqWnWI4888nQh11By\njViiFgDQZH8RqQBRl5w74YQTDq9NvH///jFD1Rvvj5KrSEVpAfJ/dNHc3M++hYWspN/33MWAelJG\nh6WzINqSczC/G5JRcu4rX/lK7S233DL90KFDVQ8++GDBVXE0oEmkorRzJOGVWoZCErtIf7Ldwlu3\nbn2yv8QKA5ecO/nkk3v6lpzbsWPHW/q7xrXXXtv17LPPbl6xYsXu5cuXn1RonGq5ilScFNCR15lN\n4cDU3EXzm+5PB/vyevyRDl87GJmacRKf0V1yLuvKK6986Qtf+MLsQmNQy1VERCTHpk2bDrdm77rr\nrkl1dXUHCr2GWq4iIiI5brrpphMfeeSRt44dO9YnTZrU09ra+rtCr6HkKiIiZaMUJefWrFkz7C5v\nJVfJU2GjTOOnUakiEp+SPXM1s9VmttfMNufs+5qZbTWzJ8zsR2Y2uVTxSKHiHGVaKI1KFZF4lbLl\n2gp8G7gjZ9+DwLXu3mNmXwWuBf6qhDFJQfIfZRqvdNwBiEiFK1nL1d0fBl7qs+8Bd+8JN9cBM0sV\nj4iIyEhJ0jPXpcBdAx00s0bCh2izZxc85UgkWi0t0B5913NzuB40w1m6cvFiaNTzZpE4JWKeq5ld\nB/QAbQOd4+4t7t7g7g21tbWlC06kP+3tkEngM+hMZkSSvogUJvaWq5l9ElgEnO/uPtT5IomRSh29\ndFEEmsIWa9HFH9LpqEIRSaRSlJzLWrNmzZSlS5fOWbt27Zb+qu4MJtbkamYLCQYwnevuBQVeUYqc\nBdPcGa663hpFEP2s4D4UzYYRkYiVouQcwMsvv1x18803n/jud7/7D8W8v2TJ1czuJBjGOc3MdgPL\nCUYHvwV4MCzps87d/2epYiob2VkwySgTm59sj6mSq4iMsKhLzgFcc801p1xzzTWd3/jGN2YUE1PJ\nkqu7X9LP7u+U6v5lr4hZME2twarr0dSXDVdwzzeIdAS3FJHkWrp0FpsjLjk3f343q+MvOffoo49O\n2LNnz/hLLrnklcQnVxERkeHqr1u4r4FKzk2ZMqW3b8m5J5544qgvCG+++SZXXXXV7O9+97sFryec\nS8lVJCFaWiAT/pNR7Lik5rA7vmmA92eLn2dLyQ0lkwnGbYkcY4gWZpyGU3Ju3759Y55++unjzjvv\nvHkAL7744riPf/zjc+++++7thQxqUnIVKbn+R6iddRb8+JwMY8dCzf50UVeee0WQXZub+3//3LkZ\ntm/PP1umUsG0WZFKMXXq1Ddffvnljdnts846a96NN974bFmNFhapTAOPUBs7FmpqIDW3yEvXBC8D\ntzZTpFKLo55BJCJ9KLmKxOLYEWpNTZBJpYPps0UPQkuHr8W+XyTZSlFyLtevf/3rbcXEqeQqFa1l\nQwvtmwqfRJxdprCpiGUKmxeG773/6PdmUvDKlLXAuQVfU0SSRclVKlr7pnYynRlSM5IxamfSy+ey\neNFofsiZodB5Ws1FrF+Sn+xk7HTOPq18ItFQcpWKl5qRKrwbdljLFPb/3uwI4cYRSyZxS/qXBq18\nItFRchWREjlc2Kog2WlD0Q/CSoevHX22RYZPyVWkkhWxbnVBa1ZnwnPThd0DUA9trpYWmjPtwVSr\ns2uKWupbSisRJedEJCbZWUFJk6GoYhWjVns7qVfWxh2FFEAtV5FKV+C61QWtWZ0usk83XdjplSAz\n6VxYpRWzSlFy7lvf+tbU5cuXz5w+ffohgMbGxr1XX331i4VcQ8lVRETKRqlKzl144YUv33HHHbuK\nDFPJ9ShF1k3NirZ+ao5yKzcnIlJCI1FybriUXHMltW5qiuTPYhCRyrKUWWwm4pJzdLOa2EvOAfzs\nZz+bfOqpp9bMmTPnjW9/+9vPzp0791Ahv4qSa19F1E3NirZ+qoiI9DXSJecAPvGJT+y78sorX5ow\nYYLfcMMNtZdeeunb1q1b91QhcSq5iojkadFzxDPYKtPM3P3AFcDZ25MxFWeIFmachlNyDmDGjBmH\nu4mvvvrqri9/+cunFBpDyabimNlqM9trZptz9p1gZg+a2dPh65RSxSMiUqgL9pLMqUsSqZ07dx4e\nBdXe3j55zpw5bxR6jVK2XFuBbwN35Oz7EvDv7n69mX0p3P6rEsYkIlKYYTw6Klq6ie0ZcqbifLzE\nAVSWG2644cSf//znk8eMGeOTJ0/uaW1t3VHoNUqWXN39YTOr77P7YxzpZLmd4K+skquIxGTgwgLN\nzTD3imx/bFOpAgplmDsXalKvUOlVk0pRcu7mm2/eA+wZTpxxr9A03d2fBwhfT4w5HhGpWItJ3lSB\no2Uy56KpA+WhbAY0mdnhVb9nz55d1DWamiAzyPOS5vBYU/rIvsWLoVHrm4pUgMELCzQ1QfP2bLds\nR4liykqzfTs0NXWMQAEDGQlxt1xfMLOTAMLXY+YaZbl7i7s3uHtDbW1tSYLLZKBd65uKiGT19vb2\nHjO6thKFn0PvQMfjbrn+BPgkcH34es9I3mzIosvp4CX7zTBbXzPJWloG/gKQCXu4ovg9sp9d0xCP\nmtTSFxnVNnd1dZ1eW1v7SlVVlQ99+ujU29trXV1dk4DNA51TsuRqZncSpK9pZrYbWE6QVL9vZsuA\nXcBFpYpntGhvD1rYSVjMO9vlruQqMjr19PRc0dnZuaqzs3M+8fd8xqkX2NzT03PFQCeUcrTwJQMc\nOr9UMYxWqVT/RUfSrcFrR4QTzgd73lMOLX0RKd6CBQv2Ah+NO45yUMnfPEREREZE3M9cRfIz2MPl\nY2SHhKeHPLO5Mzy3dehzj75FQvriRSSR1HKV8pB9uJwUqVQwektEpB9quUr5GOjh8jHS4evQ5zaF\nLVZVMhKRKKnlKiIiEjElVxERkYhVWLdwE4PXizp6Ue4hF53o++6F+Q+kCdYH1YRQEZHRSC3XWGQA\nrasoIjJaVVjLNd+maAdwZKm/fBfKbro/HZw/5OCYdJ5xiIhIOaqo5Np0fxOZzoG7hZs7g+Tb1Bpk\n1cNr87bmd/1MZ4bUjFQeczIL6T4eXPMglyp4DqcWBhYRiURFJdf/9k9ruWzr7wY8PnfXfgCarw+S\n0v5gk5r787/HiROfh8c/FWycW0ZFjbUwsBQg05khnceXtuwXvKYCF+nI/aI72JfcxWcspnGB/s5K\n8lRUcj23/lzYN2ngEzprAILWJ0fyTWpuoTc6aYhWYDp87SjwwsfK1p7tr+u6oDmcWhhY8rT4jGQs\nnpHthVJylSSqqOQ65PDfdPgaZqrBEpdIpWpc0Jh/Qit2kY5WDr8v+72vbwGKfFrOInHRaGEREZGI\nKbmKiIhETMlVREQkYkquIiIiEausAU0iMqrkOyWoP4VOE8qkYP/2ZjKdR+bCD0bThCpbIpKrmV0F\nXAE4sAm43N3fiDcqEUmypEwJ6o+mCUnBydXMJgJvuPubUQRgZqcAnwNOd/fXzez7wMUcHowvInKs\ngqYE9afAaULpNNTsD+a9D/UeTROSIZOrmVURJLslwB8BB4C3mFkXcB/Q4u5PRxDHBDM7BFQDzw3z\nepIEGaJbRjkTTnLM63o5EyJVfEhEYpDPgKaHgLcD1wIz3H2Wu58I/BdgHXC9mV1abADuvge4EdgF\nPA+84u4P9D3PzBrNbL2Zre/q6ir2dlIqi4FUzDGo+JCIxCSfbuEL3P1Q353u/hLwA+AHZjau2ADM\nbArwMeBtwD7gX8zsUnf/5z73awFaABoaGrzY+0mJNBJtizFdSImi8Nx0PueOYm3AdQRfW2cDKwn6\nn0RkxA3Zcs0mVjNb0PeYmV2Ye06RLgB+5+5d4XV+CLx/GNcTkTaCLzc7CYYJ7gy32+IMSqRyFDKg\n6TYz+6S7bwIws0sImgj/OswYdgF/bGbVwOvA+cD6YV5TJF4tDNIl3f8a14OVD8xb9tn0MoLREbm6\nw/235Z5P/N33IqNQIYtIfBy43czeaWZXAp8GPjTcANz9MeBu4HGCaThVhN2/ImWrnSNle+PQN7EO\ntD9F8HxcRCKVd8vV3Z8xs4uBHwPPAh9y99ejCMLdlwPLo7iWDEMmM0jpuaGbVZG0vAaSyUCqzJpY\nKQaoKphdgODog5FUYco+m97REXQF91U3UEwiEqV8puJsInhqk3UCMAZ4zMxw93ePVHBSQosT3nxJ\npZIfY5KsJHjG2p2zrzrcLyIjLp+W66IRj0Li19g4SHF3yKfAu+rfJkh2VLBGC4vEIp/kusvdB536\nYmY21DkiUmJLUDIViUk+yfUhM/sBcI+778ruNLPxwNnAJwkWmmgdkQhFRMrQcIoK9NXcmWF/DWRS\naZruT9G8sP8R55Ic+STXhcBS4E4zyy70MIFgVO8DwDfcPc5xkSIiiZLkogJSGkMm17A6zS3ALeFK\nTNOA191930gHN7pFs/Bu8yBfYJsXRjl8VxMiRfI17KICfbWmgwHzmY5B/5+X5CioWLq7H3L35wla\nsgCY2bzIoxr1krDwbqE0IVJEJF8FlZwzs8nAN4DTzOwN4AmCNV8uH4HYRrHoFt5tGmTJ3ab708Gx\nPEtqiYhINApKrmFX8OVm9l+BF4F3E6wFLCIiIqG8u4XN7Dgzu9rMfkjQ7DoHuNPdh7u2sIhI8r3w\nArz6KqztgPp6aFMVBBlYIc9c7wDeBfwD8G3gneE+EZHRra0NntoG3hts79wZLLqiBCsDKKRbeJ67\nn5mz/ZCZbYw6IBGRkhl0Pe0c69bR3Hsac1l1ZF93NyxbBrfdNvD7opLRaP1yU0hy/a2Z/bG7rwMw\ns/cBj45MWCJStHwTRikkueBCIWtVHxigzNBA+6OWSvFvz2m0fjkpJLm+D/gLM8uu0jQb2JJd2F8L\n+IskQNKKGyS54MKQ62nnqK+naWczzcwllVtLsK6uZItp35suyW0kIoUk14UjFoWIRKOQhCH5W7kS\n/qIKenP2VVcH+0X6UUg91/6qQ4qIjH5LlsDfvwDbqoICnHV1QWJdosoI0r+C5rmKiMQhCY+RM89P\nh7cCqTR07Ig3GEm8ykquTcBgJQY0IE8kcZLyyDaVghOfjzsKKReJSK7hsoqrgPkEnS5L3f3/lTwQ\nLZ8rkjiJeoycjjsAKReJSK7AN4H73f3jYZ3Y6hG5i6pJiIhICcSeXM3srQRLKV4G4O4HgYNxxiQi\nIjIcBZWcGyFzgC5gjZn91sxWmdnEvieZWaOZrTez9V1dXaWPUkREJE9JSK5jgfcC/+ju7wH+AHyp\n70nu3uLuDe7eUFtbW+oYRURE8hZ7tzCwG9jt7o+F23fTT3KV8teyoYX2Te0jfp/mhRnmnrAfDv4S\ngO2diwY89+/P20/N+BqiH6mSfcDf1M8xDUsXGe1ib7m6eyfwrJnNC3edDzwZY0gyQto3tZPpHGwu\nVDT+7ZkT2f5STV7n1oyv4cSJ00c4or40LF1ktEtCyxXgs0BbOFL4GeDymOOREZKakaLjso7S3Gw8\nkIHUxftKc79chxunHaW/t4jELhHJ1d0zQEPcccgoE2fjUI1TkYqWiOQqMiIawz8iIiUW+zNXSaa2\nTW3UN9dT9XdV1DfX07apLe6QREqrDagn+FeyPtwWyZNarnKMtk1tNP5rI92HugHY+cpOGv81aAIu\nOUNVQKQCtBH0enSH2zvD7dlAqce/SVlSch3lMp0Z0q3pgt6zbvc6Drx54Kh93Ye6WXbPMm7bcNvA\n9wpnl6RbB44lNUNTUCRPGeJby3cdcKDPvm5gK0qukhcl11Fs8RnFjajpm1iH2p+v1IxU0TFJhYn7\nr8lgf9Xjjk3KgpLrKNa4oJHGBYWP6KlvrmfnKzuP2V83qW7QaTTZepsdKpAgwxX3YLR6gq7gvurQ\nIDnJiwY0yTFWnr+S6nFHFyaqHlfNyvNXxhSRSImt5NjaXNXhfpE8KLnKMZacsYSWC1uom1SHYdRN\nqqPlwhYNZpLKsQRoIWipWvjaEu4XyYO6haVfS85YomQqlW0JSqZSNLVcRUREIqbkKiIiEjElVxER\nkYgpuYqIiERMyVVERCRiSq4iIiIRU3IVERGJmJKriIhIxBKTXM1sjJn91szujTsWERGR4UhMcgU+\nD2yJOwgREZHhSkRyNbOZwJ8Cq+KORUREZLgSkVyBZuCLQG/cgYiIiAxX7MnVzBYBe919wxDnNZrZ\nejNb39XVVaLoREREChd7cgU+AHzUzHYA3wPOM7N/7nuSu7e4e4O7N9TW1pY6RhERkbzFnlzd/Vp3\nn+nu9cDFwC/c/dKYwxIRESla7Mm1UrW1QX09VFUFr21tcUckIiJRSVSxdHfvADpiDmPEtbVBYyN0\ndwfbO3cG2wBLVJxZRKTsJSq5JlEmA+l0tNdctw4OHDh6X3c3LFsGt91W2LUyGUilootNRESGT8l1\nEIsXj8x1+ybWofYPJpUauThFRKQ4Sq6DaGw80l0bpfr6oCu4r7o66OiI/n4iIlJaGtAUg5Urobr6\n6H3V1cF+EREpf0quMViyBFpagpaqWfDa0qLBTCIio4W6hWOyZImSqYjIaKWWq4iISMSUXEVERCKm\n5CoiIhIxJVcREZGIKbmKiIhETMlVREQkYkquIiIiEVNyFRERiZiSq4iISMSUXEVERCKm5CoiIhIx\nJVcRkTKxdi00NcUdheQj9oX7zWwWcAcwA+gFWtz9m/FGJSKSLIsXxx2BFCL25Ar0ANe4++Nmdjyw\nwcwedPcn4w5MRCQpGhuDP1IeYu8Wdvfn3f3x8OfXgC3AKfFGJSIiUrzYk2suM6sH3gM81s+xRjNb\nb2bru7q6Sh2aiIhI3hKTXM2sBvgB0OTur/Y97u4t7t7g7g21tbWlD1CkSG1tUF8PVVXBa1vb0cem\nTQOzYLDKo48efVxEylMSnrliZuMIEmubu/8w7nhEotLWFjwn6+4OtnfuPPq52eWXw6FDR7Z7emDp\n0uDnJUtKF6eIRCv25GpmBnwH2OLuN8Udj4w+LS3Q3h7PvdetgwMHjt7X3Q3LlgU/5ybWrIMH4brr\nlFxFylkSuoU/APw5cJ6ZZcI/H4k7KBk92tshk4nn3n0Ta+7+gY4B7No1MvGISGnE3nJ1918CFncc\nMrqlUtDRUfr71tcHXcF91dUFr/0dA5g9e8RCEpESSELLVWTUWrkSqquP3lddHexfuRLGjTv2PePH\nB8dEpHzF3nIVGc2yz02vuy7o6p09O0icuc9TP/95+P3vg5+nToVvflPPW0XKnZKryAhbsmTgZDnY\nMREpX+oWFhERiZiSq5TMYIspxHktEZGoqVtYSmKwxRQK7RaN8loiIiNByVUilclAOn3s/sEWU7jt\ntsLuUei1MplgKo6ISKkouUpkBqs3OdhiCoUq9FqplGphikhpmbvHHUPBGhoafP369XGHIQUYbDGF\nHTviu5ZIJTGzDe7eEHcclUADmqQkBltMIc5riYiMBCVXKYklS4IF9OvqgvJqdXXBdjEDkKK8lojI\nSFC3sIhIhVC3cOmo5SoiIhIxJVcREZGIKbmKiIhETMlVREQkYkquIiIiEVNyFRERiVgikquZLTSz\nbWa23cy+FHc8pbRixYq4QxjUSMQX5TWT/vnlGizWFStWFPy75Ht+OX1GSdT389PnKfmIfZ6rmY0B\nngI+COwGfgNc4u5PDvSe0TTP1cyI+7/BYEYiviivmfTPL9dgsZoZQEG/S76/ezl9RknU9/Mr589T\n81xLJwkt17OA7e7+jLsfBL4HfCzmmERERIqWhOR6CvBszvbucN9RzKzRzNab2fqurq6SBSciIlKo\nJCRX62ffMX0u7t7i7g3u3lBbW1uCsERERIqThOS6G5iVsz0TeC6mWERERIYtCcXSfwO8w8zeBuwB\nLgYqprT18uXL4w5hUCMRX5T/FSicAAAFRklEQVTXTPrnl2uwWIv5PfJ9Tzl9RknU9/PT5yn5iH20\nMICZfQRoBsYAq9190Mqco2m0sIhIqWi0cOkkoeWKu98H3Bd3HCIiIlFIwjNXERGRUUXJVUREJGJK\nriIiIhFTchUREYlYIkYLF8rMuoCdRb59GvBihOFEJalxQXJjU1yFUVyFSWpcUHxsde6uVXhKoCyT\n63CY2fokDkVPalyQ3NgUV2EUV2GSGhckOzYJqFtYREQkYkquIiIiEavE5NoSdwADSGpckNzYFFdh\nFFdhkhoXJDs2oQKfuYqIiIy0Smy5ioiIjCglVxERkYhVVHI1s4Vmts3MtpvZl+KOB8DMVpvZXjPb\nHHcsucxslpk9ZGZbzOw/zOzzcccEYGbHmdmvzWxjGNffxR1TLjMbY2a/NbN7444ll5ntMLNNZpYx\ns8SUlDKzyWZ2t5ltDf+u/UkCYpoXfk7ZP6+aWVPccQGY2VXh3/vNZnanmR0Xd0zSv4p55mpmY4Cn\ngA8SFGj/DXCJuz8Zc1znAPuBO9x9fpyx5DKzk4CT3P1xMzse2AD8WQI+LwMmuvt+MxsH/BL4vLuv\nizOuLDO7GmgA3urui+KOJ8vMdgAN7p6oRRHM7HbgEXdfZWbjgWp33xd3XFnhvxt7gPe5e7EL10QV\nyykEf99Pd/fXzez7wH3u3hpnXNK/Smq5ngVsd/dn3P0g8D3gYzHHhLs/DLwUdxx9ufvz7v54+PNr\nwBbglHijAg/sDzfHhX8S8Q3RzGYCfwqsijuWcmBmbwXOAb4D4O4Hk5RYQ+cD/xl3Ys0xFphgZmOB\nauC5mOORAVRScj0FeDZnezcJSBblwMzqgfcAj8UbSSDses0Ae4EH3T0RcQHNwBeB3rgD6YcDD5jZ\nBjNrjDuY0BygC1gTdqWvMrOJcQfVx8XAnXEHAeDue4AbgV3A88Ar7v5AvFHJQCopuVo/+xLR4kky\nM6sBfgA0ufurcccD4O5vunsKmAmcZWaxd6eb2SJgr7tviDuWAXzA3d8LfBj4TPg4Im5jgfcC/+ju\n7wH+ACRiLARA2E39UeBf4o4FwMymEPS2vQ04GZhoZpfGG5UMpJKS625gVs72TNSlMqjwmeYPgDZ3\n/2Hc8fQVdiF2AAtjDgXgA8BHw2eb3wPOM7N/jjekI9z9ufB1L/AjgsckcdsN7M7pebibINkmxYeB\nx939hbgDCV0A/M7du9z9EPBD4P0xxyQDqKTk+hvgHWb2tvAb6cXAT2KOKbHCgUPfAba4+01xx5Nl\nZrVmNjn8eQLBPzhb440K3P1ad5/p7vUEf7d+4e6JaFWY2cRwUBpht+uHgNhHp7t7J/Csmc0Ld50P\nxDpgro9LSEiXcGgX8MdmVh3+/3k+wVgISaCxcQdQKu7eY2Z/CfwcGAOsdvf/iDkszOxOIA1MM7Pd\nwHJ3/068UQFBS+zPgU3h802Av3b3+2KMCeAk4PZwFGcV8H13T9S0lwSaDvwo+PeYsUC7u98fb0iH\nfRZoC7/wPgNcHnM8AJhZNcHMgk/FHUuWuz9mZncDjwM9wG/RMoiJVTFTcUREREqlkrqFRURESkLJ\nVUREJGJKriIiIhFTchUREYmYkquIiEjElFxFIhJWePl03HGISPyUXEWiMxlQchURJVeRCF0PvD2s\nAfq1uIMRkfhoEQmRiITVg+5NUl1eEYmHWq4iIiIRU3IVERGJmJKrSHReA46POwgRiZ+Sq0hE3P33\nwKNmtlkDmkQqmwY0iYiIREwtVxERkYgpuYqIiERMyVVERCRiSq4iIiIRU3IVERGJmJKriIhIxJRc\nRUREIvb/AX4+GXDy+ByPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c50bb00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdcAAAEWCAYAAAA0MN3QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XuUXGWZ7/Hvk86FJA0kQCeYEBMR\nCSLRaDJeGSaLyxycCTrnHHGJOuogk1lnHDXA0qNzzghzXIyXYTDjiHNgkEQlIAPeEFHBIxHEEU0w\nkEguMJoAIZcWuTXB3Po5f1R1qHS6kq7u3b2r09/PWr2qa9eut556a1f/+t21a7+RmUiSpOKMKLsA\nSZIONYarJEkFM1wlSSqY4SpJUsEMV0mSCma4SpJUMMNVkqSCGa4aliLi1Ij4aUQ8HRG/i4h7IuIP\nennfjIgTBqiuEyPi2xHRXq3rBxExcyAeS9LAMVw17ETEEcCtwL8ARwFTgb8HdpRZV9UE4BZgJjAZ\n+Dnw7VIrktSw8AxNGm4iYi7ww8yccIB1zgc+AhxLJeAWZObGiLgL+ENgO5DA+4GxwAWZeWrN/RN4\nWWY+HBFLgN8DLwVeD9wHvCczN/ai1qOAJ4BjMvOJvjxfSYPPkauGo/XAnoj4ckS8OSIm1t4YEX8G\n/C3w34A24G7gBoDMPK262qsyszUzb+zlY74L+CRwDLASWNrL+50GbDFYpaHFcNWwk5nPAKdSGXn+\nG9AeEbdExOTqKn8FfCoz12TmbuAfgNkRMb0fD/vdzLwrM3cA/wt4Q0RMO9AdIuI44Ergon48rqQS\nGK4alqrB+b7MPA44BZgCLKrePB3454h4KiKeAn4HBJXPZvvq0ZrH7qi2OaXeyhHRBtwOfDEzb+jH\n40oqgeGqYS8z1wJLqIQsVILwrzJzQs3P2Mz8aZ0mngPGdV2JiGN7WGdaze2tVA6kerynxqq7qW8H\nbsnMyxp9PpLKZ7hq2ImIkyLi4upuV6q7Z88DflZd5f8CH4+IV1RvPzIizq1pYitwfM31+4FXRMTs\niDgMuLSHh/2T6td/RlP57PXezHy0+0rVI5l/ANyTmR/r1xOVVBrDVcPRs8DrgHsj4jkqoboauBgg\nM78JfAb4WkQ8U73tzTX3vxT4cnW38dszcz3wf4AfAg8BP+nhMa8HLqGyO3gOlQOcevJfgT8A/iIi\nOmp+XtyfJyxpcPlVHGmAVb+K81hm/u+ya5E0OBy5SpJUMMNVkqSCuVtYkqSCOXKVJKlgI8suoC+O\nOeaYnDFjRtllSNKQsmLFit9mZls/7j9p5MiR11D5TvhwHpx1Aqt37959wZw5c7b1tMKQDNcZM2aw\nfPnyssuQpCElIg46WcSBjBw58ppjjz325W1tbU+OGDFi2H6m2NnZGe3t7Sdv2bLlGuAtPa0znP/z\nkCQ15pS2trZnhnOwAowYMSLb2tqe5oWzuu2/ziDWI0ka2kYM92DtUu2HuhlquEqSVDDDVZI0ZLS0\ntMw56aSTTu76Wbdu3ei77rpr3Pve974DTuHYiO9973utJ5988stHjhw5Z/HixRMPfo/9DckDmiRJ\nw9OYMWM6165d+2DtspkzZ+487bTTtndfd9euXYwaNarhxzj++ON3Ll68eMOnP/3pyQdfu2eDFq4R\ncS0wH9iWmadUlx0F3AjMADYAb8/MJwerJknS0Hfrrbce/k//9E+T77zzzocvuuiiKZs3bx71yCOP\njD7qqKN2n3XWWc/ccsstEzo7O2PdunVjP/CBD2zZuXPniBtvvPHo0aNHd95+++0PTZ48eU9tezNn\nztwJMGJE33fuDubIdQnwBeArNcs+Bvy/zPx0RHysev1/DmJNkqQ+OP98pq1e/cI8xkU45RS2X3st\n+03FWGvHjh0jTjrppJMBpk2btuOOO+74z+7rPPDAA+Puvffeta2trfn5z3/+6PXr14+9//77H3z+\n+edHzJw585S/+7u/27RmzZoH3//+90+76qqrjv7EJz7R43dV+2PQwjUz74qIGd0WvxWYV/39y8Ay\nygzXq6lMDNZfmx+HrQW/VpN+CFNu3WfR45th29b6d9kxejO7Rh9ghRo3zZzEl2dNqXv7e9/wOOf2\n/F3pfYwaBWPGwNFffS8Tv3nu/re3jGJMy5he1bSPzTtg664G7jAKGAOTfsjjcWvdfmqkj7o7UJ+9\nd9V8zl135n7LJx+xg0mH138eETBmVB/7CHrZT9W+6dLDttWlp22s0T67aeYkWECvtp8uk7/2XiZ9\n54XtJwJ6GkT0eXuCOn3VrW9q7NgJu3bCk6N/yBOje+6vvev2Y7vqyU0zJ3HSJ7YzaXwb/33uLwtr\ndyjqabdwd2efffZTra2te49qfuMb3/jsxIkTOydOnNjZ2tq659xzz30KYNasWdsfeOCBQv9B6FL2\nZ66TM3MzQGZujohJ9VaMiAXAAoAXv3iApra8HlgJzO5nO1u3QUcHtLYWUBTQcULlstsfwG1bD/ww\nu0ZvZU9LBy17DlzHrPYOgAOG67lztjFragerNtVva091x8qYMTDxm+cy9lezeP4Vq164PffAHvr2\nx3DrLujYA60tvVi5WkjHKwDYxq11+6m3fdTdwfrs3HVnMqv9BFa1PbzP8kmH72L8mD08t2P/55EJ\nGbvZ1dc+gl70U9fer2r7dbatLj1tY430WVc/MYeDbj+1Jn3nXMavmcVzL19FvdOf92t7gh76qlvf\ndLNrJ7TumQ07OWi49nW76klXH+45Fkb3u7XiHGyEWabx48d31l4fPXr03q1oxIgRHHbYYdn1++7d\nu2Mgaig7XHstM6+mMrZk7ty5A/c9q9lUxs/9MW9h5XJZfxvqag9g9n7tLZx34IeZt6SywrL3HaSO\nefM4FXhq0YHWq7R16gn115lXW89o4NXQuuzUxus5aOMHXbl6sQyYzcLqC9rTXftc08H6bCVwApy6\nrPt/apXHO7KHjWzePFg5ex6zZ/exj7oagQP0U/X2rsefBz1tW1162sYa6rNqP1HN8ANtP/sYC7wG\njlx2at2n1K/tqVrbvg1Xr9f5AzB/HixaCbNnz2b2QbbDfte2T2OVPlz5fP+b0uAp+6s4WyPiRQDV\ny8L3e0uS1Igf//jH4yZPnvzK2267beKFF144/YQTTnhFo22UPXK9BXgv8Onq5bfLLUeS1My2b9++\n34fO8+fPf3b+/PnPAlxxxRWP1972oQ996Angia7rmzZtWlXvti5/9Ed/tH3r1q0P9KfOQRu5RsQN\nwH8AMyPisYh4P5VQPSsiHgLOql6XJGlIG8yjhc+rc9MZg1WDJEmDoezPXCVJOuQYrpIkFcxwlSSp\nYIarJEkFM1wlSUPGYEw5d+mll05+6Utf+ooTTzzx5De84Q0nrl+/vuGTY5X9PVdJknptMKacmzNn\nzvaLL754zeGHH975mc98pu3CCy887rvf/e6vG2nDcJUkDWlFTzl3zjnnPNv1+6mnntpx4403Ht1o\nTYarJKlh53/7/Gmrt60udsq5Sadsv/at1zbVlHNXXXVV25lnnvl0o8/FcJUkDRmDOeXcF7/4xaPu\nv//+cVddddW6Rus0XCVJDTvYCLNMRUw5961vfevwyy+//EV33333urFjxzY8E5vhKklSjXvuuWfs\nBz/4wem33XbbQ1OnTt3dlzYMV0mSanzkIx+Ztn379pZzzz33pQBTpkzZ+aMf/ejhRtowXCVJQ8Zg\nTDn305/+dH1/6/QkEpIkFcxwlSSpYIarJEkFM1wlSSqY4SpJUsEMV0mSCma4SpKGjMGYcu6zn/1s\n24knnnjySSeddPKcOXNmrlix4rBG2/B7rpKkIWMwppy74IILnvjoRz/aDrB06dIjFy5cOO3uu+9+\nqJE2DFdJ0pBW9JRzRx111N5zE3d0dLRE9Hj64QMyXCVJfXD+NCh2yjk4ZTs0x5Rzn/rUp9q++MUv\nTt61a9eIO+64o+FZcfzMVZI0ZHTtFl67du2DPQUr1J9ybsqUKbu7Tzm3YcOGMT218fGPf7z90Ucf\nXX3ppZc+dskll7yo0ToduUqS+uDQnnKuy1/+5V/+7iMf+ciLG63BkaskSTVWrVq1dzR74403Hjl9\n+vQdjbbhyFWSpBpXXHHFpLvvvvuIkSNH5pFHHrl7yZIlv2m0DcNVkjRkDMaUc4sXL+73Lm93C0uS\nVDDDVZKkghmukiQVrCnCNSIujIhfRcTqiLghIho+j6MkSc2i9HCNiKnAh4C5mXkK0AK8o9yqJEnq\nu9LDtWokMDYiRgLjgMcPsr4kSU2r9HDNzE3A5cAjwGbg6cy8vft6EbEgIpZHxPL29vbBLlOS1AQG\nY8q5LosXL54YEXPuuuuuhs+hXPr3XCNiIvBW4CXAU8BNEfHuzLyudr3MvBq4GmDu3Lm5X0OSpEPe\nYEw5B/Dkk0+OuPLKKye98pWvfK4v9y89XIEzgd9kZjtARHwDeCNw3QHvJUkSxU85B3DxxRdPvfji\ni7d87nOfO7YvNTVDuD4CvD4ixgHPA2cAy8stSZJ0QOefP43VBU85d8op27m2/Cnn7rnnnrGbNm0a\nfd555z09ZMM1M++NiJuB+4DdwC+p7v6VJKlWT7uFu6s35dzEiRM7u08598ADD+zzD8KePXu48MIL\nX/zVr3614fMJ1yo9XAEy8xLgkrLrkCT10kFGmGXqz5RzTz31VMtDDz102Omnnz4T4Le//e2ot73t\nbSfcfPPND/f0uW49TRGukiQ1g6OPPnrPk08+eX/X9de+9rUzL7/88kcbCVZogq/iSJJ0qHHkKkka\nMgZjyrlaP//5z9f1pU5HrpIkFcxwlSSpYIarJEkFM1wlSSqY4SpJUsEMV0mSCma4SpKGjMGYcu7z\nn//80RMnTnxV12NcccUVxzTaht9zlSQNGYM15dw555zz5Fe+8pVH+lim4SpJGtoGYsq5/jJcJUmN\nO59prKbgKefYzrWUPuUcwPe+970JJ554Yuvxxx//+y984QuPnnDCCbsaeSp+5ipJGjK6dguvXbv2\nwZ6CFepPOTdlypTd3aec27Bhw5ju93/729/+1COPPLJq/fr1D55++unPvvvd735Jo3U6cpUkNe4g\nI8wy9WfKOYBjjz12727iiy66qP2Tn/zk1EZrcOQqSVKNjRs37j0K6vrrr59w/PHH/77RNhy5SpJU\n47Of/eykH/zgBxNaWlpywoQJu5csWbKh0TYMV0nSkDEYU85deeWVm4BN/anT3cKSJBXMcJUkqWCG\nqySptzo7Ozv3O7p2OKr2Q2e92w1XSVJvrW5vbz9yuAdsZ2dntLe3HwmsrreOBzRJknpl9+7dF2zZ\nsuWaLVu2nMLwHpx1Aqt37959Qb0VDFdJUq/MmTNnG/CWsusYCobzfx6SJA0Iw1WSpIIZrpIkFWxY\nfea6cCGsXFn/9kXV2xbOe2HZO98JCxYMaFmSpEOMI9cDWLkSrr++7CokSUPNsBq5Llp0kBXmVS6W\nLatenTdwtUiSDl1NMXKNiAkRcXNErI2INRHxhrJrkiSpr5pl5PrPwPcz820RMRoYV3ZBkiT1Venh\nGhFHAKcB7wPIzJ3AzjJrkiSpP5pht/DxQDuwOCJ+GRHXRMT47itFxIKIWB4Ry9vb2we/SkmSeqkZ\nwnUk8BrgXzPz1cBzwMe6r5SZV2fm3Myc29bWNtg1SpLUa80Qro8Bj2XmvdXrN1MJW0mShqTSP3PN\nzC0R8WhEzMzMdcAZwIMD8VgLv7+QlVvqn0Vi0ZbKd3UWLlkIwMrZleXzlvT+Md456514zglJGt5K\nD9eqDwJLq0cK/xr4i5Lr6ZOu4DZcJWl4a4pwzcyVwNyBfpxFZx/kLBJLKhfL3rcMeOEkEssOdvKJ\nqnlL5vWhKknSoaYZPnOVJOmQYrhKklQww1WSpIIZrpIkFazhcI2I8RHRMhDFSJJ0KDhouEbEiIh4\nZ0R8NyK2AWuBzRHxq4j4x4h42cCXKUnS0NGbkeudwEuBjwPHZua0zJwE/CHwM+DTEfHuAaxRkqQh\npTffcz0zM3d1X5iZvwO+Dnw9IkYVXtmAWAjUP0MTLKpZrxeTq3e/99nVtpd0LZl3gLXfiaebkKRD\n00FHrl3BGhFzut8WEefUrqPeWglcX3YRkqQB0sgZmv4tIt6bmasAIuI8KkO87wxIZQOit0PRZQAs\nXFi9tqx391r4/Xk19973t33N62UdkqShqJFwfRtwc0S8CzgVeA/wxwNSlSRJQ1ivwzUzfx0R7wC+\nBTwK/HFmPj9glUmSNEQdNFwjYhWQNYuOAlqAeyOCzHzlQBUnSdJQ1Juv4swHzqn5eR2V3cFd19Un\nW4EZVF6CGcDSMovZ31L2LW9r/5vcuhV+9jMYMQJmzICl/XjKS1ctZcaiGYz4+xHMWDSDpauarP/2\nUduZP2NvZy5dWumICBg5En68DJ55ptJRRTzqkOqj5tP18nRtr08/BDwD/JimfMuqufRmt/AjmZkH\nWiEi4mDrqNZWYD3QWb2+kRe+lvOuUirax1Iq5WyvXt8IrKv+PrmPTS6F9euhsxNI2LgRFlSf8rsa\nfMpLVy1lwXcWsH1XpcCNT29kwXcqjb1rVhP03z66d+YOYD0s/WtY8GXYXl2+Z0/lMjsrHbV0aeMd\nU/uodfrozJe8mMmtfXwRh5GlSyvbZ9fL88aNlc/C9mqyt6yaTxwsEyNiGZXvs347Mx+pWT6ayoFN\n7wXuzMwlA1fmvubOnZvLly9v/I6v/jH85sj6t3ecAK0Pw+zKYcIrV0JHB7S29q75jtaV7GnpYPYW\nWNXWyvy3ze5xvVv/ZgWzpm5n1aax3W4JKnvc9zXrXbcyfs0snnv5qn2WZ1YGPS11Tka5J/fQEi20\njj7IE1jZAR17oLXa0B72/SAAoGUPdLZAS/22OjoqGdHSAmP3nMDzLQ/zcGulL595BjonrYTRHbCz\npo2Akb04meas9o69ffrMjmfI7NxvnYgRHDHmCG79m5XMmtoB85Yxfs0sOk5aVbefet1H3XXvs+46\nZkHrKph9Gvt3JpUB7I5uTfIqnpq8gRj9LC27xtV/YQ+gq59OO+f4Hvvozov+k1dPe55VmyrPt962\n1aWnbayhPtvbT1DZtnv7ZnrhvVjvfdj1fmvZ0+BrV1W7TQF7t5tVm1rZvXvfdVuo/Cs8+z230rp2\nFs+fVO2vnt+yfd+uelLtwxwPz7/oIcatvaDPTUXEiswc8Lmz1buR69nA+cANEfES4ClgLJV9XLcD\nn6tOdj70tT4Mk3649+qkBv/BH7VzMoyGVW1w08xJdde7acUEevyD2+My2HbOTfTUWkTlp56WaGFU\nSy/O7zG52zo9lbGnBXaNqrzydYwaDeys/P58y8M8OfqFvuxM4LkeOrSX+ztWtbXu7dOeQqN2+U0r\nKutNrvZbUL+fet1H3XXvs+5aV8Gkm6j7BHfsv2gS29jx3OHsYnf9+x1EVz/V66Prfz6BkSNeSIN6\n21aXnraxhvpsn35qoJ9r3ov13odd77e+qt2m4IXtpiddGfrE/Jv2zdI6L1Oft6ueVPuwsxN2/H4i\n44ppVQPsoCPXfVaunInpGOD5zHxqwKo6iD6PXJvGDCr7lbqbDmwY1Ep6NIPCy5sxo7IreL8mp8OG\nBtucsWgGG5/ev7HpR05nw8IGGxtwM+ixM2e0wMY99e/Wl46pbX5I9VHz6b69/obKK7mfJnnL9pYj\n18HT0Kw4mbkrMzdTGckCEBEzC6/qkHcZ7Pf/57jq8iYwAOVddhmM69bmuHGV5Q23dcZljBu1b2Pj\nRo3jsjOapP/2UaczL1uwf4fsvbmPHVP7qEOqj5pP9+31b4Hnuq/URG9ZNaHM7PUPMAFYDPwH8NdU\nPnNd3EgbRfzMmTMnh77rMnN6Zkb18royi9nfAJR33XWZ06dnRlQur+tHm9c9cF1O/9z0jEsjp39u\nel73QJP13z7qdGZXh0BmS0vlsr8dU/uoQ6qPmk/37fXu/5FN/ZbtDWB5DvLf6+H609Bu4S4R8V+A\n3wKvBH6bmYN6CsShv1tYkgafu4UHT6/P0BQRh/HCaDWBnwD/mpm/H6DaJEkakho5t/BXgGeBf6le\nP6+67O1FFyVJ0lDWSLjOzMxX1Vy/MyLuL7ogSZKGukaOFv5lRLy+60pEvA64p/iSJEka2hoZub4O\neE9EdJ2l6cXAmq4T+6cn8JckCWgsXM8esCokSTqENDKfa0/n7JEkSd00dIYmSZJ0cIarJEkFM1wl\nSSpY04RrRLRExC8j4taya5EkqT+aJlyBDwNryi5CkqT+aopwjYjjgD8Frim7FkmS+qspwhVYBHwU\n6Ky3QkQsiIjlEbG8vb198CqTJKlBpYdrRMwHtmXmigOtl5lXZ+bczJzb1tY2SNVJktS40sMVeBPw\nlojYAHwNOD0iriu3JEmS+q70cM3Mj2fmcZk5A3gH8KPMfHfJZUmS1Gelh6skSYeaRk7cP+Aycxmw\nrOQyJEnqF0eukiQVzHCVJKlghqskSQUzXCVJKpjhKklSwQxXSZIKZrhKklQww1WSpIIZrpIkFcxw\nlSSpYIarJEkFM1wlSSqY4SpJUsEMV0mSCma4SpJUMMNVkqSCGa6SJBXMcJUkqWCGqyRJBTNcJUkq\nmOEqSVLBDFdJkgpmuEqSVDDDVZKkghmukiQVzHCVJKlghqskSQUzXCVJKpjhKklSwQxXSZIKVnq4\nRsS0iLgzItZExK8i4sNl1yRJUn+MLLsAYDdwcWbeFxGHAysi4o7MfLDswiRJ6ovSR66ZuTkz76v+\n/iywBphablWSJPVd6eFaKyJmAK8G7u3htgURsTwilre3tw92aZIk9VrThGtEtAJfBxZm5jPdb8/M\nqzNzbmbObWtrG/wCJUnqpaYI14gYRSVYl2bmN8quR5Kk/ig9XCMigC8BazLzirLrkSSpv0oPV+BN\nwJ8Dp0fEyurPn5RdlCRJfVX6V3Ey8ydAlF2HJElFaYaRqyRJhxTDVZKkghmukiQVzHCVJKlghqsk\nSQUzXCVJKpjhKklSwQxXSZIKZrhKklQww1WSpIIZrpIkFcxwlSSpYIarJEkFM1wlSSqY4SpJUsEM\nV0mSCma4SpJUMMNVkqSCGa6SJBXMcJUkqWCGqyRJBTNcJUkqmOEqSVLBDFdJkgpmuEqSVDDDVZKk\nghmukiQVzHCVJKlghqskSQUzXCVJKpjhKklSwZoiXCPi7IhYFxEPR8THyq5nMF166aVll3BAA1Ff\nkW02e//Vqldr1/JGn0tv1x9KfdSMuvef/aneiMwst4CIFmA9cBbwGPAL4LzMfLDefebOnZvLly8f\npAoHVkRQ9mtwIANRX5FtNnv/1apXa9fyRp9Lb9cfSn3UjLr331Duz4hYkZlzy65jOGiGketrgYcz\n89eZuRP4GvDWkmuSJKnPmiFcpwKP1lx/rLpsHxGxICKWR8Ty9vb2QStOkqRGNUO4Rg/L9tvnkplX\nZ+bczJzb1tY2CGVJktQ3zRCujwHTaq4fBzxeUi2SJPVbM4TrL4CXRcRLImI08A7glpJrGjSXXHJJ\n2SUc0EDUV2Sbzd5/terV2rW80efS2/WHUh81o+79Z3+qN0o/WhggIv4EWAS0ANdm5mUHWv9QOlpY\nkgaLRwsPnpFlFwCQmbcBt5VdhyRJRWiG3cKSJB1SDFdJkgpmuEqSVDDDVZKkgjXF0cKNioh2YGMf\n734M8NsCyymKdTWuWWuzrsZYV2P6U9f0zPQsPINgSIZrf0TE8mY8FN26GtestVlXY6yrMc1al/bl\nbmFJkgpmuEqSVLDhGK5Xl11AHdbVuGatzboaY12Nada6VGPYfeYqSdJAG44jV0mSBpThKklSwYZV\nuEbE2RGxLiIejoiPlV0PQERcGxHbImJ12bXUiohpEXFnRKyJiF9FxIfLrgkgIg6LiJ9HxP3Vuv6+\n7JpqRURLRPwyIm4tu5ZaEbEhIlZFxMqIaJoppSJiQkTcHBFrq9vaG5qgppnVfur6eSYiFpZdF0BE\nXFjd7ldHxA0RcVjZNalnw+Yz14hoAdYDZ1GZoP0XwHmZ+WDJdZ0GdABfycxTyqylVkS8CHhRZt4X\nEYcDK4A/a4L+CmB8ZnZExCjgJ8CHM/NnZdbVJSIuAuYCR2Tm/LLr6RIRG4C5mdlUJ0WIiC8Dd2fm\nNdX5nMdl5lNl19Wl+ndjE/C6zOzriWuKqmUqle395Mx8PiL+HbgtM5eUWZd6NpxGrq8FHs7MX2fm\nTuBrwFtLronMvAv4Xdl1dJeZmzPzvurvzwJrgKnlVgVZ0VG9Oqr60xT/IUbEccCfAteUXctQEBFH\nAKcBXwLIzJ3NFKxVZwD/WXaw1hgJjI2IkcA44PGS61EdwylcpwKP1lx/jCYIi6EgImYArwbuLbeS\niuqu15XANuCOzGyKuoBFwEeBzrIL6UECt0fEiohYUHYxVccD7cDi6q70ayJifNlFdfMO4IayiwDI\nzE3A5cAjwGbg6cy8vdyqVM9wCtfoYVlTjHiaWUS0Al8HFmbmM2XXA5CZezJzNnAc8NqIKH13ekTM\nB7Zl5oqya6njTZn5GuDNwAeqH0eUbSTwGuBfM/PVwHNAUxwLAVDdTf0W4KayawGIiIlU9ra9BJgC\njI+Id5dbleoZTuH6GDCt5vpxuEvlgKqfaX4dWJqZ3yi7nu6quxCXAWeXXArAm4C3VD/b/BpwekRc\nV25JL8jMx6uX24BvUvmYpGyPAY/V7Hm4mUrYNos3A/dl5tayC6k6E/hNZrZn5i7gG8AbS65JdQyn\ncP0F8LKIeEn1P9J3ALeUXFPTqh449CVgTWZeUXY9XSKiLSImVH8fS+UPztpyq4LM/HhmHpeZM6hs\nWz/KzKYYVUTE+OpBaVR3u/4xUPrR6Zm5BXg0ImZWF50BlHrAXDfn0SS7hKseAV4fEeOq788zqBwL\noSY0suwCBktm7o6IvwF+ALQA12bmr0oui4i4AZgHHBMRjwGXZOaXyq0KqIzE/hxYVf18E+BvM/O2\nEmsCeBHw5epRnCOAf8/MpvraSxOaDHyz8veYkcD1mfn9ckva64PA0uo/vL8G/qLkegCIiHFUvlnw\nV2XX0iUz742Im4H7gN3AL/FUiE1r2HwVR5KkwTKcdgtLkjQoDFdJkgpmuEqSVDDDVZKkghmukiQV\nzHCVClKd4eWvy65DUvkMV6nEosodAAAAgElEQVQ4EwDDVZLhKhXo08BLq3OA/mPZxUgqjyeRkApS\nnT3o1maal1dSORy5SpJUMMNVkqSCGa5ScZ4FDi+7CEnlM1ylgmTmE8A9EbHaA5qk4c0DmiRJKpgj\nV0mSCma4SpJUMMNVkqSCGa6SJBXMcJUkqWCGqyRJBTNcJUkq2P8HXDQZ34nfQ5QAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x119a3b240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdcAAAEWCAYAAAA0MN3QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAH8xJREFUeJzt3X2UVPWd5/HPt2lQEAXUFgSBFg0w\nBqUTehOTOMLxIYckGOfsUVfA2RgwnbOTNQE5JjpMIrM5btQkhHHVia0BNNMoWc2DMZrITuQhGkjA\nNA+RB43hUZDWiKKtQNPf/aOqtGyqu+t236rfra736xxOd926detTP271t363fvf+zN0FAADiUxE6\nAAAAPQ3FFQCAmFFcAQCIGcUVAICYUVwBAIgZxRUAgJhRXAEAiBnFFWXJzC4ws2fN7A0z+5uZPWNm\n/yXPx7qZnV2gXKems7xmZgfM7Pdm9qlCPBeAwqkMHQAoNjM7SdLjkv6HpJ9I6iPp7yUdCpkr7S1J\nMyS9IMklXS7pl2Z2mru3BE0GIG/0XFGORkuSuz/k7kfd/R13f8rdN2RWMLMZZrbZzF43s9+Y2cj0\n8pXpVdab2Vtm9t/M7Foz+132E2T3bs1ssZn90MyWmdlBM1uR2V5b7v6uu29191ZJJumopEGSTo69\nFQAUDMUV5WibpKNm9oCZfcbMBmXfaWb/IOmfJf1XSVWSVkl6SJLc/cL0auPdvb+7L83zOadL+rak\nUyU1SmroaGUz2yDpXUmPSbrf3ffn+TwAEoDiirLj7m9KukCpw673SWoys8fMbHB6lS9L+o67b04f\niv3fkmra623m6VfuvtLdD0maK+kTZja8g4znSTpJ0jRJv2tvPQDJRHFFWUoXzmvd/QxJ4yQNlbQg\nffdISf+WHlB0QNLflDpEO6wbT7kr67nfSm9zaCcZ33X3hyTdZGbju/HcAIqM4oqy5+5bJC1WqshK\nqUL4ZXcfmPWvr7s/284m3pbUL3PDzIbkWGd41v39lfoO9eU8I/aWNCrPdQEkAMUVZcfMxprZHDM7\nI317uKSpklanV/mhpJvN7MPp+weY2ZVZm3hFHyx26yV92MxqzOx4SfNyPO1n06f/9FHqu9c17r6r\n7Upmdn5mPTPra2bfkDRY0ppuvWgARUVxRTk6KOnjktaY2dtKFdVNkuZIkrv/TNLtkh42szfT930m\n6/HzJD2QPmx8lbtvk/S/JP0/pU6hyfUd6RJJtyh1OHiCUgOccjlO0t2SXpO0R9JnJX3O3fPt5QJI\nAGOydKCwzGyxpN3u/i+hswAoDnquAADEjOIKAEDMOCwMAEDM6LkCABCzkrxw/6mnnurV1dWhYwBA\nSVm3bt2r7l7VjcefVllZeb9S54SXc+esVdKmlpaW6yZMmJDz0qQlWVyrq6u1du3a0DEAoKSY2Y7u\nPL6ysvL+IUOG/F1VVdXrFRUVZfudYmtrqzU1NZ2zb9+++yV9Ptc65fzJAwAQzbiqqqo3y7mwSlJF\nRYVXVVW9ofev6nbsOkXMAwAobRXlXlgz0u3Qbg2luAIAEDOKKwCgZPTq1WvC2LFjz8n827p1a5+V\nK1f2u/baa9udwjGqJ598sv8555zzd5WVlRMWLVo0qPNHHKskBzQBAMrTcccd17ply5bns5eNGTPm\n8IUXXtjcdt0jR46od+/ekZ9j1KhRhxctWrT9tttuG9z52rlRXAEAJe3xxx8/8fvf//7gp59++sUb\nbrhh6N69e3vv3Lmzz8knn9xy6aWXvvnYY48NbG1tta1bt/b9yle+su/w4cMVS5cuPaVPnz6tTz31\n1AuDBw8+mr29MWPGHJakioquH9yluAIAIpsxQ8M3bXp/HuM4jBun5oULdcxUjNkOHTpUMXbs2HMk\nafjw4YeWLVv2l7brbNiwod+aNWu29O/f3++8885Ttm3b1nf9+vXPv/POOxVjxowZ981vfnPP5s2b\nn585c+bwe++995RvfetbOc9V7Y7yKq6zJDWGDgEA3VAjaUHoEOHkOizc1uTJkw/079//vVHNn/zk\nJw8OGjSoddCgQa39+/c/euWVVx6QpHPPPbd5w4YNsX5AyCiv4goAiEVnPcyQTjjhhNbs23369Hmv\n0FZUVOj444/3zO8tLS1WiAzlVVzL+NMeAKB4OBUHAIAsK1as6Dd48ODznnjiiUGzZ88eefbZZ384\n6jZKcsq52tpa59rCABCNma1z99quPn79+vXbx48f/2qcmUrZ+vXrTx0/fnx1rvvouQIAEDOKKwAA\nMaO4AgAQM4orAAAxo7gCABAziisAADGjuAIASkYxppybN2/e4LPOOuvDo0ePPucTn/jE6G3btvWJ\nuo3yukITAKCkFWPKuQkTJjTPmTNn84knnth6++23V82ePfuMX/3qVy9F2UbRiquZLZQ0RdJ+dx+X\nXvZdSZdJOizpL5K+6O4HipUJAFD64p5y7rLLLjuY+f2CCy54a+nSpadEzVTMnutiSXdJejBr2TJJ\nN7t7i5ndLulmSd8oYiYAQBfM+MWM4Zv2b4p3yrnTxjUvvHxhoqacu/fee6suueSSN6K+lqIVV3df\naWbVbZY9lXVztaQripUnrHpJS0KH6IJpkupCh8hDKbRvqLYM3TZJ3odCt02+ktyGhVfMKefuueee\nk9evX9/v3nvv3Ro1Z5K+c50haWl7d5pZndJ71IgRI4qVqUCWKDWxbE3oIBFkJsIthTd10ts3ZFuG\nbJuk70NJ328kaUX6Z/g27KyHGVIcU879/Oc/P/F73/ve6atWrdrat2/fyBfhT0RxNbO5klokNbS3\njrvXK/XRUrW1taU328AxaiQtDx0igkmhA0SU5PadFPj5Q7XNpADPGVWS9xupNNqw9D3zzDN9r7/+\n+pFPPPHEC8OGDWvpyjaCF1cz+4JSA50u9lKcogcA0KPceOONw5ubm3tdeeWVZ0nS0KFDD//2t799\nMco2ghZXM5us1ACmie5+zDBqAACyNTc3/6ntsilTphycMmXKQUmaP3/+y9n3ffWrX31N0muZ23v2\n7NnY3n0Zzz777Lbu5izaRSTM7CFJv5c0xsx2m9lMpUYPnyhpmZk1mtkPi5UHAIBCKeZo4ak5Fv+o\nWM8PAECxcPlDAABiRnEFACBmFFcAAGJGcQUAIGYUVwBAySjGlHN33HFH1ejRo88ZO3bsORMmTBiz\nbt2646NuI/hFJAAAyFcxppy77rrrXvv617/eJEkNDQ0DZs2aNXzVqlUvRNkGxRUAUNLinnLu5JNP\nfu/axG+99VYvs5yXH+4QxRUA0AUzhkvxTjknjWuWkjHl3He+852qe+65Z/CRI0cqli1bFnlWHL5z\nBQCUjMxh4S1btjyfq7BK7U85N3To0Ja2U85t3779uFzbuPnmm5t27dq1ad68ebtvueWW06PmLK+e\n66xZUmNj5+sVXCbDpJAhIiqlzEnPGjJfuT53PpKeT5IapZr+0oLQOaTOepghxTHlXMaXvvSlv914\n442R5zml5woAQJaNGze+15tdunTpgJEjRx6Kuo3y6rkuSMTHPb3/yXh5wAxRTUr/XB4wQ74mpX8u\nD5ihI5PSP5fz3IkyKf1zecAMnZkUOkBZmD9//mmrVq06qbKy0gcMGNCyePHiv0bdRnkVVwBASSvG\nlHOLFi3q9iFvDgsDABAziisAADGjuAIAEDOKKwAAMaO4AgAQM4orAAAxo7gCAEpGMaacy1i0aNEg\nM5uwcuXKyNdQ5jxXAEDJKMaUc5L0+uuvV9x9992nnXfeeW935fEUVwBASYt7yjlJmjNnzrA5c+bs\n+8EPfjCkK5mKVlzNbKGkKZL2u/u49LKTJS2VVC1pu6Sr3P31YmUCAHTRjBnDtSnmKefGjWvWwvBT\nzj3zzDN99+zZ02fq1KlvJL64Slos6S5JD2Ytu0nSf7r7bWZ2U/r2N4qYCQBQQnIdFm6rvSnnBg0a\n1Np2yrkNGzZ84APC0aNHNXv27BE//vGPI19POFvRiqu7rzSz6jaLL9f7V6J+QKkrZlNcS1R9vbRk\nSegU78/PMGvW+8umTZPq6sLkyVcx2i9X23QmdNsVa7/Kp21CtkV9vfSxj0n9+0tnnx0mwwd00sMM\nqTtTzh04cKDXCy+8cPxFF100RpJeffXV3ldcccXZjzzyyIu5vtdtT+jRwoPdfa8kpX+e1t6KZlZn\nZmvNbG1TU1PRAiJ/S5YkZLrcNhobk1H0O7NkibRiRegUH5SEtkvKfhW6LZYskQ4cCPf85eKUU045\n+vrrr6/fs2fPxj179mwcP37821ELq1RCA5rcvV5SvSTV1tZ6J6sjkJoaafny0ClSMjkmTQqZIpqJ\nE4vTfvk+R1Larpj7VXvPk4S2GDgwIb1WdCp0cX3FzE53971mdrqk/Z0+AgBQtoox5Vy2P/zhD1u7\nkjP0YeHHJH0h/fsXJP0iYBYAAGJRtOJqZg9J+r2kMWa228xmSrpN0qVm9oKkS9O3AQAoacUcLTy1\nnbsuLlYGAACKIfRhYQAAehyKKwAAMaO4AgAQM4orAKBkFGPKuTvvvPOUQYMGjc88x/z580+Nuo3Q\n57kCAJC3Yk05d9lll73+4IMP7uxiTIorAKC0FWLKue6iuAIAopuh4dqkmKecU7MWKviUc5L05JNP\nDhw9enT/UaNGvXvXXXftOvvss49EeSl85woAKBmZw8Jbtmx5Pldhldqfcm7o0KEtbaec2759+3Ft\nH3/VVVcd2Llz58Zt27Y9f9FFFx285pprzoyak54rACC6TnqYIXVnyjlJGjJkyHuHiW+44Yamb3/7\n28OiZqDnCgBAlh07drw3CmrJkiUDR40a9W7UbdBzBQAgyx133HHab37zm4G9evXygQMHtixevHh7\n1G1QXAEAJaMYU87dfffdeyTt6U5ODgsDABAziisAADGjuAIA8tXa2tp6zOjacpRuh9b27qe4AgDy\ntampqWlAuRfY1tZWa2pqGiBpU3vrMKAJAJCXlpaW6/bt23f/vn37xqm8O2etkja1tLRc194KFFcA\nQF4mTJiwX9LnQ+coBeX8yQMAgIKguAIAEDOKKwAAMUtEcTWz2Wb2ZzPbZGYPmdnxoTMBANBVwYur\nmQ2T9FVJte4+TlIvSVeHTQUAQNcFL65plZL6mlmlpH6SXu5kfQAAEit4cXX3PZK+J2mnpL2S3nD3\np9quZ2Z1ZrbWzNY2NTUVOyYAAHkLXlzNbJCkyyWdKWmopBPM7Jq267l7vbvXunttVVVVsWMCAJC3\n4MVV0iWS/uruTe5+RNJPJX0ycCYAALosCcV1p6TzzayfmZmkiyVtDpwJAIAuC15c3X2NpEckPSdp\no1KZ6oOGAgCgGxJxbWF3v0XSLaFzAAAQh+A9VwAAehqKKwAAMaO4AgAQM4orAAAxo7gCABAziisA\nADGjuAIAELPIxdXMTjCzXoUIAwBAT9BpcTWzCjObZma/MrP9krZI2pue3Py7ZvahwscEAKB05NNz\nfVrSWZJuljTE3Ye7+2mS/l7Sakm35ZrFBgCAcpXP5Q8vSc9W8wHu/jdJj0p61Mx6x54MAIAS1WnP\nNVNYzWxC2/vM7LLsdQAAQLQBTfeZ2bmZG2Y2VdK/xB8JAIDSFmVWnCskPWJm0yVdIOm/S/p0QVIB\nAFDC8i6u7v6SmV0t6eeSdkn6tLu/U7BkAACUqE6Lq5ltlORZi06W1EvSGjOTu59XqHAAAJSifHqu\nUwqeAgCAHiSf4rrT3b2jFczMOlsHAIBykddFJMzsejMbkb3QzPqY2UVm9oCkLxQmHgAApSefnutk\nSTMkPWRmZ0o6IKmvUoX5KUk/cPfGwkWMz6xZUmPEpNOmSXV1hckDAOiZOi2u7v6upHsk3ZO+EtOp\nkt5x9wOFDhdaphBTXAEAUUQ5zzVzJaa9ZnaDpPmSZGZj3H1rIcLFbcGCaOtPmlSQGACAHi5ScTWz\ngZJ+IGmsmb0raYOkmZK+2J0Q6e3eL2mcUqf9zHD333dnmwAAhBJpPld3P+DuX5Q0T9IaSR+S9NMY\ncvybpF+7+1hJ4yVtjmGb6NEaJFUrtQtXp2+Xi3xfezm3ERBW3j1XMzte0j8pdelDl/Q7Sf+e/k62\ny8zsJEkXSrpWktz9sKTD3dkmeroGSXWSmtO3d6RvS9L0IImKJ9/X3t56IyQNLnxMoMxFOSz8oKSD\nkv5P+vbU9LKruplhlKQmSYvMbLykdZK+5u5vZ69kZnVK/xUZMWLEMRvJzyxJ+Q8Xjvodbf4aJdUU\nauMF1ChpUrv3Fq692lot6VCbZc1KfUNxn0qjfY9ty/zar7PX3tl6WxW2uHa8D+VSvP2qFPYblIoo\nh4XHuPtMd386/a9O0pgYMlRK+qhSveCPSHpb0k1tV3L3enevdffaqqqqGJ42pBpJ00KHiGiakvOH\np23RaLs86e3bnbbs7LV3tp4rXNskaR/KJen7jfTiizVKdhsiI0rP9U9mdr67r5YkM/u4pGdiyLBb\n0m53X5O+/YhyFNc4zPq11Lgv//Uzp+LURNiXp507TXUTeuK5O+8dOGjXrFmpn8uXFzpLtVKHOdsa\nKangTx6D3G2ZX/tVK7/X3tF6ofbPzvehXIq3XyXfXXct0BVXhE6BfETpuX5c0rNmtt3Mtkv6vaSJ\nZrbRzDZ0NYC775O0y8wyveCLJT3f1e2F1LivUUs2LgkdowzcKqlfm2X90st7unxfezm3ERBelJ7r\n5IKlkK6X1GBmfSS9pG6e2tOeBZOjfXmTOc91eZ4Pm7R4UqTto6syA3fmStqp1CCdW9XzBzNJ+b/2\ncm4jILwo87nmOsYUi/TlE2sLtX30RNNVvoUi39dezm0EhBXpPFcAANA5iisAADGLdPlDoDvq19UH\nGfDVmB7tPWlx7vuTNML75dPrNWlxcgbFtdd2SWizYu9Pne1HbSWhjRAOxbUTjY35X8D/vTdfnuvH\npVSmxVuycYka9zWqZkhyztNr3Jc63yopfwT3D16i/Qlro7aS0mZJ3J8yktJGCIfi2oFpyT6fXFLp\nTYtXM6RGy69dXtTn7GjUdxJHeIdoo/bkarsktVkx2yrK2QNJaiOEQXHtQF1dtKKVOVyU76k7cWBa\nPABIHgY0AQAQM4orcmpokKqrpYqK1M8GZisrmExbr1ghvfmm9MoroRMhl1dekVav5j2B/FBccYyG\nhtTh8B07JPfUz7o6/pgUQnZbS6n23raNtk6ahobU/8uhQ7wnkB++cw2kYWOD5v7nXO18Y6dGDBih\nWy++VdPP7drVdKKMaM7H6tWpPyLZmpulmTOl++7L/ZhMjiiTHBRTe20UaoR3Rq62bm2V5s6Vpifk\n4kpt267YbZaE/Wru3NT/S7aO3hOFaKMVK6SJE+PbHgqL4hpAw8YG1f2yTs1HUhNZ73hjh+p+mRo5\nFbXAFmJEc9s/9p0tz6ipSeYI6yRmymivTXfuLG6O9iSh7ZKwX7X3/9HZeyJOEyeGbwfkj+Ias8Z9\njZ0Ow1+9e7UOHf3gu7L5SLNm/mKm7lvXQdcwlz6Srm3/7q6cyF5d/f5hymwjR5bmtF8djfoOMcI7\nW3ttPWJE0aPklKvtQrdZCCNGRHtPlGMb4YP4zjVG086dltcJ7W0La2fLu6qrU+DdeqvUr81sZf36\npZYjXrnauqKCtk4a3hOIip5rjOom1OXVS6xeUK0dbxz7MXjkgJGxnhDf1RPZM9/1zZ2bOhw2YkTq\nj0hSvgPsSdq2dZ/jpDPPpK2ThvcEoqLnGsCtF9+qfr0/+DG4X+9+uvXi5HwMnj5d2r49NYhj+3b+\niBRSdluff740eHDoRMiF9wSioLgGMP3c6aq/rF4jB4yUyTRywEjVX1bf5dHCSLaGjQ2qXlCtin+t\nUPWCajVs5PwNoKfjsHAg08+dTjEtA3GODAdQOiiuPVw+o5eLJckzmBSqjaKODF+xY4Umjkz+yYxJ\n2K+Suj9lFKqNaobUaMFkhiEnHcW1B5t2brJOiqsZUpO4TIXOE3Vk+MSRExPXRm0lJV8S96eMpOZC\n8Zi7h84QWW1tra9duzZ0DKBTHY0M3z5re/EDoayZ2Tp3rw2doxwwoAkooFIYGQ4gfokprmbWy8z+\nZGaPh84CxIWR4UB5StJ3rl+TtFnSSaGDAHFiZDhQfhLRczWzMyR9TtL9obMAANBdiSiukhZI+rqk\n1vZWMLM6M1trZmubmpqKlwwAgIiCF1czmyJpv7uv62g9d69391p3r62qqipSOgAAogteXCV9StLn\nzWy7pIclXWRm/xE2EgAAXRe8uLr7ze5+hrtXS7pa0m/d/ZrAsQAA6LLgxRUAgJ4mSafiyN2XS1oe\nOAYAAN1CzxUAgJhRXAEAiBnFFQCAmFFcAQCIGcUVAICYUVwBAIgZxRUAgJhRXAEAiBnFFQCAmFFc\nAQCIGcUVAICYUVwBAIgZxRUAgJhRXAEAiBnFFQCAmFFcAQCIGcUVAICYUVwBAIgZxRUAgJhRXAEA\niBnFFQCAmAUvrmY23MyeNrPNZvZnM/ta6EwAAHRHZegAklokzXH358zsREnrzGyZuz8fOhgAAF0R\nvOfq7nvd/bn07wclbZY0LGwqAAC6LnhxzWZm1ZI+ImlN2CQAAHRdYoqrmfWX9KikWe7+Zo7768xs\nrZmtbWpqKn5AAADylIjiama9lSqsDe7+01zruHu9u9e6e21VVVVxAwIAEEHw4mpmJulHkja7+/zQ\neQAA6K7gxVXSpyT9o6SLzKwx/e+zoUMBANBVwU/FcfffSbLQOQAAiEsSeq4AAPQoFFcAAGJGcQUA\nIGYUVwAAYkZxBQAgZhRXAABiRnEFACBmFFcAAGJGcQUAIGYUVwAAYkZxBQAgZhRXAABiRnEFACBm\nFFcAAGJGcQUAIGYUVwAAYkZxBQAgZhRXAABiRnEFACBmFFcAAGJGcQUAIGYUVwAAYpaI4mpmk81s\nq5m9aGY3hc5TTPPmzQsdoUOFyBfnNpPeftk6yjpv3rzIryXf9UupjZKobfvRnsiHuXvYAGa9JG2T\ndKmk3ZL+KGmquz/f3mNqa2t97dq1RUpYWGam0P8HHSlEvji3mfT2y9ZRVjOTpEivJd/XXkptlERt\n26+U29PM1rl7begc5SAJPdePSXrR3V9y98OSHpZ0eeBMAAB0WRKK6zBJu7Ju704v+wAzqzOztWa2\ntqmpqWjhAACIKgnF1XIsO+aYi7vXu3utu9dWVVUVIRYAAF2ThOK6W9LwrNtnSHo5UBYAALqtMnQA\npQYwfcjMzpS0R9LVkqaFjVQ8t9xyS+gIHSpEvji3mfT2y9ZR1q68jnwfU0ptlERt24/2RD6CjxaW\nJDP7rKQFknpJWujut3a0fk8aLQwAxcJo4eJJQs9V7v6EpCdC5wAAIA5J+M4VAIAeheIKAEDMKK4A\nAMSM4goAQMwSMVo4KjNrkrSjiw8/VdKrMcaJS1JzScnNRq5oyBVNUnNJXc820t25Ck8RlGRx7Q4z\nW5vEoehJzSUlNxu5oiFXNEnNJSU7G1I4LAwAQMworgAAxKwci2t96ADtSGouKbnZyBUNuaJJai4p\n2dmgMvzOFQCAQivHnisAAAVFcQUAIGZlVVzNbLKZbTWzF83sptB5JMnMFprZfjPbFDpLNjMbbmZP\nm9lmM/uzmX0tdCZJMrPjzewPZrY+netfQ2fKZma9zOxPZvZ46CzZzGy7mW00s0YzS8yUUmY20Mwe\nMbMt6X3tEwnINCbdTpl/b5rZrNC5JMnMZqf3+01m9pCZHR86E3Irm+9czayXpG2SLlVqgvY/Sprq\n7s8HznWhpLckPeju40JmyWZmp0s63d2fM7MTJa2T9A8JaC+TdIK7v2VmvSX9TtLX3H11yFwZZnaD\npFpJJ7n7lNB5Msxsu6Rad0/URRHM7AFJq9z9fjPrI6mfux8InSsj/Xdjj6SPu3tXL1wTV5ZhSu3v\n57j7O2b2E0lPuPvikLmQWzn1XD8m6UV3f8ndD0t6WNLlgTPJ3VdK+lvoHG25+153fy79+0FJmyUN\nC5tK8pS30jd7p/8l4hOimZ0h6XOS7g+dpRSY2UmSLpT0I0ly98NJKqxpF0v6S+jCmqVSUl8zq5TU\nT9LLgfOgHeVUXIdJ2pV1e7cSUCxKgZlVS/qIpDVhk6SkD702StovaZm7JyKXpAWSvi6pNXSQHFzS\nU2a2zszqQodJGyWpSdKi9KH0+83shNCh2rha0kOhQ0iSu++R9D1JOyXtlfSGuz8VNhXaU07F1XIs\nS0SPJ8nMrL+kRyXNcvc3Q+eRJHc/6u41ks6Q9DEzC3443cymSNrv7utCZ2nHp9z9o5I+I+kr6a8j\nQquU9FFJ/+7uH5H0tqREjIWQpPRh6s9L+r+hs0iSmQ1S6mjbmZKGSjrBzK4JmwrtKafiulvS8Kzb\nZ4hDKh1Kf6f5qKQGd/9p6DxtpQ8hLpc0OXAUSfqUpM+nv9t8WNJFZvYfYSO9z91fTv/cL+lnSn1N\nEtpuSbuzjjw8olSxTYrPSHrO3V8JHSTtEkl/dfcmdz8i6aeSPhk4E9pRTsX1j5I+ZGZnpj+RXi3p\nscCZEis9cOhHkja7+/zQeTLMrMrMBqZ/76vUH5wtYVNJ7n6zu5/h7tVK7Vu/dfdE9CrM7IT0oDSl\nD7t+WlLw0enuvk/SLjMbk150saSgA+bamKqEHBJO2ynpfDPrl35/XqzUWAgkUGXoAMXi7i1m9j8l\n/UZSL0kL3f3PgWPJzB6SNEnSqWa2W9It7v6jsKkkpXpi/yhpY/r7TUn6Z3d/ImAmSTpd0gPpUZwV\nkn7i7ok67SWBBkv6WervsSolLXH3X4eN9J7rJTWkP/C+JOmLgfNIksysn1JnFnw5dJYMd19jZo9I\nek5Si6Q/icsgJlbZnIoDAECxlNNhYQAAioLiCgBAzCiuAADEjOIKAEDMKK4AAMSM4grEJD3Dyz+F\nzgEgPIorEJ+BkiiuACiuQIxuk3RWeg7Q74YOAyAcLiIBxCQ9e9DjSZqXF0AY9FwBAIgZxRUAgJhR\nXIH4HJR0YugQAMKjuAIxcffXJD1jZpsY0ASUNwY0AQAQM3quAADEjOIKAEDMKK4AAMSM4goAQMwo\nrgAAxIziCgBAzCiuAADE7P8DyJmSBD0/bw0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11db1a400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "last_time = 8\n",
    "adjustments = [t for t in range(num_price_adjusts_before_horizon[0,0] - 1) if traject_out[0,0,0][t][0][1]<last_time]\n",
    "colors = [\"blue\",\"green\",\"yellow\",\"red\",\"magenta\"]\n",
    "firms = range(num_firm)\n",
    "arrivals = np.concatenate(np.array([arrival_out[0,0][t] for t in adjustments if arrival_out[0,0][t].size]).squeeze())\n",
    "if arrivals.ndim == 1:\n",
    "    arrivals = (np.array([arrival_out[0,0][t] for t in adjustments if arrival_out[0,0][t].size]).squeeze())\n",
    "arrivals = np.array([arrival for arrival in arrivals if arrival[0]<last_time])\n",
    "for setup in range(len(setups)):\n",
    "    for color, firm in zip(colors, firms):\n",
    "        plt.step([traject_out[setup,0,0][t][0][1] for t in adjustments], [traject_out[setup,0,0][t][firm][0] for t in adjustments], c = color, label='Firm {}'.format(firm+1), where=\"post\")\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "    buying_data = np.concatenate([x for x in [buying_out[setup,0,0][t] for t in adjustments] if x != []])\n",
    "    buying_data = np.array([buying for buying in buying_data if buying[0]<last_time])\n",
    "    for color, firm in zip(colors, firms):\n",
    "        buying_data_firm = np.array([buying_data[x] for x in range(len(buying_data)) if buying_data[x,3]==firm])\n",
    "        if len(buying_data_firm>0):\n",
    "            plt.scatter(buying_data_firm[:,0] ,buying_data_firm[:,1], c = color)\n",
    "    plt.plot(arrivals[:,0],[0 for _ in arrivals[:,0]],'k|')\n",
    "    plt.xlabel('t')\n",
    "    plt.ylabel('p$_t$(k)')\n",
    "    plt.title('Setup {}'.format(setup+1))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}