How to read the posterior file.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to read the posterior file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior samples are stored as datasets in the `samples` group of `posterior.hdf`. Below we demonstrate how to read the samples and make a simple histogram of the 1D marginal distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "from matplotlib import pyplot\n",
    "\n",
    "import numpy\n",
    "import h5py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp = h5py.File('posterior.hdf', 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The available parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dec\n",
      "distance\n",
      "inclination\n",
      "lambda1\n",
      "lambda2\n",
      "loglikelihood\n",
      "mass1\n",
      "mass2\n",
      "polarization\n",
      "ra\n",
      "spin1z\n",
      "spin2z\n",
      "tc\n"
     ]
    }
   ],
   "source": [
    "for param in sorted(fp['samples'].keys()):\n",
    "    print param"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: In the posterior file, `mass1` is defined to be $\\geq$ `mass2`. Likewise, `spin1z` and `lambda1` are the spin and tidal parameter of the more massive compact object. All quantities are in the detector frame."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the 1D marginal distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we plot the 1D marginal distribution of the distance, with vertical lines at the 5, 50, and 95th percentiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "/* Put everything inside the global mpl namespace */\n",
       "window.mpl = {};\n",
       "\n",
       "\n",
       "mpl.get_websocket_type = function() {\n",
       "    if (typeof(WebSocket) !== 'undefined') {\n",
       "        return WebSocket;\n",
       "    } else if (typeof(MozWebSocket) !== 'undefined') {\n",
       "        return MozWebSocket;\n",
       "    } else {\n",
       "        alert('Your browser does not have WebSocket support.' +\n",
       "              'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
       "              'Firefox 4 and 5 are also supported but you ' +\n",
       "              'have to enable WebSockets in about:config.');\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
       "    this.id = figure_id;\n",
       "\n",
       "    this.ws = websocket;\n",
       "\n",
       "    this.supports_binary = (this.ws.binaryType != undefined);\n",
       "\n",
       "    if (!this.supports_binary) {\n",
       "        var warnings = document.getElementById(\"mpl-warnings\");\n",
       "        if (warnings) {\n",
       "            warnings.style.display = 'block';\n",
       "            warnings.textContent = (\n",
       "                \"This browser does not support binary websocket messages. \" +\n",
       "                    \"Performance may be slow.\");\n",
       "        }\n",
       "    }\n",
       "\n",
       "    this.imageObj = new Image();\n",
       "\n",
       "    this.context = undefined;\n",
       "    this.message = undefined;\n",
       "    this.canvas = undefined;\n",
       "    this.rubberband_canvas = undefined;\n",
       "    this.rubberband_context = undefined;\n",
       "    this.format_dropdown = undefined;\n",
       "\n",
       "    this.image_mode = 'full';\n",
       "\n",
       "    this.root = $('<div/>');\n",
       "    this._root_extra_style(this.root)\n",
       "    this.root.attr('style', 'display: inline-block');\n",
       "\n",
       "    $(parent_element).append(this.root);\n",
       "\n",
       "    this._init_header(this);\n",
       "    this._init_canvas(this);\n",
       "    this._init_toolbar(this);\n",
       "\n",
       "    var fig = this;\n",
       "\n",
       "    this.waiting = false;\n",
       "\n",
       "    this.ws.onopen =  function () {\n",
       "            fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n",
       "            fig.send_message(\"send_image_mode\", {});\n",
       "            if (mpl.ratio != 1) {\n",
       "                fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n",
       "            }\n",
       "            fig.send_message(\"refresh\", {});\n",
       "        }\n",
       "\n",
       "    this.imageObj.onload = function() {\n",
       "            if (fig.image_mode == 'full') {\n",
       "                // Full images could contain transparency (where diff images\n",
       "                // almost always do), so we need to clear the canvas so that\n",
       "                // there is no ghosting.\n",
       "                fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n",
       "            }\n",
       "            fig.context.drawImage(fig.imageObj, 0, 0);\n",
       "        };\n",
       "\n",
       "    this.imageObj.onunload = function() {\n",
       "        fig.ws.close();\n",
       "    }\n",
       "\n",
       "    this.ws.onmessage = this._make_on_message_function(this);\n",
       "\n",
       "    this.ondownload = ondownload;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_header = function() {\n",
       "    var titlebar = $(\n",
       "        '<div class=\"ui-dialog-titlebar ui-widget-header ui-corner-all ' +\n",
       "        'ui-helper-clearfix\"/>');\n",
       "    var titletext = $(\n",
       "        '<div class=\"ui-dialog-title\" style=\"width: 100%; ' +\n",
       "        'text-align: center; padding: 3px;\"/>');\n",
       "    titlebar.append(titletext)\n",
       "    this.root.append(titlebar);\n",
       "    this.header = titletext[0];\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_canvas = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var canvas_div = $('<div/>');\n",
       "\n",
       "    canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n",
       "\n",
       "    function canvas_keyboard_event(event) {\n",
       "        return fig.key_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    canvas_div.keydown('key_press', canvas_keyboard_event);\n",
       "    canvas_div.keyup('key_release', canvas_keyboard_event);\n",
       "    this.canvas_div = canvas_div\n",
       "    this._canvas_extra_style(canvas_div)\n",
       "    this.root.append(canvas_div);\n",
       "\n",
       "    var canvas = $('<canvas/>');\n",
       "    canvas.addClass('mpl-canvas');\n",
       "    canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n",
       "\n",
       "    this.canvas = canvas[0];\n",
       "    this.context = canvas[0].getContext(\"2d\");\n",
       "\n",
       "    var backingStore = this.context.backingStorePixelRatio ||\n",
       "\tthis.context.webkitBackingStorePixelRatio ||\n",
       "\tthis.context.mozBackingStorePixelRatio ||\n",
       "\tthis.context.msBackingStorePixelRatio ||\n",
       "\tthis.context.oBackingStorePixelRatio ||\n",
       "\tthis.context.backingStorePixelRatio || 1;\n",
       "\n",
       "    mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n",
       "\n",
       "    var rubberband = $('<canvas/>');\n",
       "    rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n",
       "\n",
       "    var pass_mouse_events = true;\n",
       "\n",
       "    canvas_div.resizable({\n",
       "        start: function(event, ui) {\n",
       "            pass_mouse_events = false;\n",
       "        },\n",
       "        resize: function(event, ui) {\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "        stop: function(event, ui) {\n",
       "            pass_mouse_events = true;\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "    });\n",
       "\n",
       "    function mouse_event_fn(event) {\n",
       "        if (pass_mouse_events)\n",
       "            return fig.mouse_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    rubberband.mousedown('button_press', mouse_event_fn);\n",
       "    rubberband.mouseup('button_release', mouse_event_fn);\n",
       "    // Throttle sequential mouse events to 1 every 20ms.\n",
       "    rubberband.mousemove('motion_notify', mouse_event_fn);\n",
       "\n",
       "    rubberband.mouseenter('figure_enter', mouse_event_fn);\n",
       "    rubberband.mouseleave('figure_leave', mouse_event_fn);\n",
       "\n",
       "    canvas_div.on(\"wheel\", function (event) {\n",
       "        event = event.originalEvent;\n",
       "        event['data'] = 'scroll'\n",
       "        if (event.deltaY < 0) {\n",
       "            event.step = 1;\n",
       "        } else {\n",
       "            event.step = -1;\n",
       "        }\n",
       "        mouse_event_fn(event);\n",
       "    });\n",
       "\n",
       "    canvas_div.append(canvas);\n",
       "    canvas_div.append(rubberband);\n",
       "\n",
       "    this.rubberband = rubberband;\n",
       "    this.rubberband_canvas = rubberband[0];\n",
       "    this.rubberband_context = rubberband[0].getContext(\"2d\");\n",
       "    this.rubberband_context.strokeStyle = \"#000000\";\n",
       "\n",
       "    this._resize_canvas = function(width, height) {\n",
       "        // Keep the size of the canvas, canvas container, and rubber band\n",
       "        // canvas in synch.\n",
       "        canvas_div.css('width', width)\n",
       "        canvas_div.css('height', height)\n",
       "\n",
       "        canvas.attr('width', width * mpl.ratio);\n",
       "        canvas.attr('height', height * mpl.ratio);\n",
       "        canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n",
       "\n",
       "        rubberband.attr('width', width);\n",
       "        rubberband.attr('height', height);\n",
       "    }\n",
       "\n",
       "    // Set the figure to an initial 600x600px, this will subsequently be updated\n",
       "    // upon first draw.\n",
       "    this._resize_canvas(600, 600);\n",
       "\n",
       "    // Disable right mouse context menu.\n",
       "    $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n",
       "        return false;\n",
       "    });\n",
       "\n",
       "    function set_focus () {\n",
       "        canvas.focus();\n",
       "        canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    window.setTimeout(set_focus, 100);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>')\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items) {\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) {\n",
       "            // put a spacer in here.\n",
       "            continue;\n",
       "        }\n",
       "        var button = $('<button/>');\n",
       "        button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n",
       "                        'ui-button-icon-only');\n",
       "        button.attr('role', 'button');\n",
       "        button.attr('aria-disabled', 'false');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "\n",
       "        var icon_img = $('<span/>');\n",
       "        icon_img.addClass('ui-button-icon-primary ui-icon');\n",
       "        icon_img.addClass(image);\n",
       "        icon_img.addClass('ui-corner-all');\n",
       "\n",
       "        var tooltip_span = $('<span/>');\n",
       "        tooltip_span.addClass('ui-button-text');\n",
       "        tooltip_span.html(tooltip);\n",
       "\n",
       "        button.append(icon_img);\n",
       "        button.append(tooltip_span);\n",
       "\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    var fmt_picker_span = $('<span/>');\n",
       "\n",
       "    var fmt_picker = $('<select/>');\n",
       "    fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n",
       "    fmt_picker_span.append(fmt_picker);\n",
       "    nav_element.append(fmt_picker_span);\n",
       "    this.format_dropdown = fmt_picker[0];\n",
       "\n",
       "    for (var ind in mpl.extensions) {\n",
       "        var fmt = mpl.extensions[ind];\n",
       "        var option = $(\n",
       "            '<option/>', {selected: fmt === mpl.default_extension}).html(fmt);\n",
       "        fmt_picker.append(option)\n",
       "    }\n",
       "\n",
       "    // Add hover states to the ui-buttons\n",
       "    $( \".ui-button\" ).hover(\n",
       "        function() { $(this).addClass(\"ui-state-hover\");},\n",
       "        function() { $(this).removeClass(\"ui-state-hover\");}\n",
       "    );\n",
       "\n",
       "    var status_bar = $('<span class=\"mpl-message\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {\n",
       "    // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n",
       "    // which will in turn request a refresh of the image.\n",
       "    this.send_message('resize', {'width': x_pixels, 'height': y_pixels});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_message = function(type, properties) {\n",
       "    properties['type'] = type;\n",
       "    properties['figure_id'] = this.id;\n",
       "    this.ws.send(JSON.stringify(properties));\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_draw_message = function() {\n",
       "    if (!this.waiting) {\n",
       "        this.waiting = true;\n",
       "        this.ws.send(JSON.stringify({type: \"draw\", figure_id: this.id}));\n",
       "    }\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    var format_dropdown = fig.format_dropdown;\n",
       "    var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n",
       "    fig.ondownload(fig, format);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_resize = function(fig, msg) {\n",
       "    var size = msg['size'];\n",
       "    if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {\n",
       "        fig._resize_canvas(size[0], size[1]);\n",
       "        fig.send_message(\"refresh\", {});\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_rubberband = function(fig, msg) {\n",
       "    var x0 = msg['x0'] / mpl.ratio;\n",
       "    var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;\n",
       "    var x1 = msg['x1'] / mpl.ratio;\n",
       "    var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;\n",
       "    x0 = Math.floor(x0) + 0.5;\n",
       "    y0 = Math.floor(y0) + 0.5;\n",
       "    x1 = Math.floor(x1) + 0.5;\n",
       "    y1 = Math.floor(y1) + 0.5;\n",
       "    var min_x = Math.min(x0, x1);\n",
       "    var min_y = Math.min(y0, y1);\n",
       "    var width = Math.abs(x1 - x0);\n",
       "    var height = Math.abs(y1 - y0);\n",
       "\n",
       "    fig.rubberband_context.clearRect(\n",
       "        0, 0, fig.canvas.width, fig.canvas.height);\n",
       "\n",
       "    fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_figure_label = function(fig, msg) {\n",
       "    // Updates the figure title.\n",
       "    fig.header.textContent = msg['label'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_cursor = function(fig, msg) {\n",
       "    var cursor = msg['cursor'];\n",
       "    switch(cursor)\n",
       "    {\n",
       "    case 0:\n",
       "        cursor = 'pointer';\n",
       "        break;\n",
       "    case 1:\n",
       "        cursor = 'default';\n",
       "        break;\n",
       "    case 2:\n",
       "        cursor = 'crosshair';\n",
       "        break;\n",
       "    case 3:\n",
       "        cursor = 'move';\n",
       "        break;\n",
       "    }\n",
       "    fig.rubberband_canvas.style.cursor = cursor;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_message = function(fig, msg) {\n",
       "    fig.message.textContent = msg['message'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_draw = function(fig, msg) {\n",
       "    // Request the server to send over a new figure.\n",
       "    fig.send_draw_message();\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_image_mode = function(fig, msg) {\n",
       "    fig.image_mode = msg['mode'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Called whenever the canvas gets updated.\n",
       "    this.send_message(\"ack\", {});\n",
       "}\n",
       "\n",
       "// A function to construct a web socket function for onmessage handling.\n",
       "// Called in the figure constructor.\n",
       "mpl.figure.prototype._make_on_message_function = function(fig) {\n",
       "    return function socket_on_message(evt) {\n",
       "        if (evt.data instanceof Blob) {\n",
       "            /* FIXME: We get \"Resource interpreted as Image but\n",
       "             * transferred with MIME type text/plain:\" errors on\n",
       "             * Chrome.  But how to set the MIME type?  It doesn't seem\n",
       "             * to be part of the websocket stream */\n",
       "            evt.data.type = \"image/png\";\n",
       "\n",
       "            /* Free the memory for the previous frames */\n",
       "            if (fig.imageObj.src) {\n",
       "                (window.URL || window.webkitURL).revokeObjectURL(\n",
       "                    fig.imageObj.src);\n",
       "            }\n",
       "\n",
       "            fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n",
       "                evt.data);\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "        else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == \"data:image/png;base64\") {\n",
       "            fig.imageObj.src = evt.data;\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        var msg = JSON.parse(evt.data);\n",
       "        var msg_type = msg['type'];\n",
       "\n",
       "        // Call the  \"handle_{type}\" callback, which takes\n",
       "        // the figure and JSON message as its only arguments.\n",
       "        try {\n",
       "            var callback = fig[\"handle_\" + msg_type];\n",
       "        } catch (e) {\n",
       "            console.log(\"No handler for the '\" + msg_type + \"' message type: \", msg);\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        if (callback) {\n",
       "            try {\n",
       "                // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n",
       "                callback(fig, msg);\n",
       "            } catch (e) {\n",
       "                console.log(\"Exception inside the 'handler_\" + msg_type + \"' callback:\", e, e.stack, msg);\n",
       "            }\n",
       "        }\n",
       "    };\n",
       "}\n",
       "\n",
       "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n",
       "mpl.findpos = function(e) {\n",
       "    //this section is from http://www.quirksmode.org/js/events_properties.html\n",
       "    var targ;\n",
       "    if (!e)\n",
       "        e = window.event;\n",
       "    if (e.target)\n",
       "        targ = e.target;\n",
       "    else if (e.srcElement)\n",
       "        targ = e.srcElement;\n",
       "    if (targ.nodeType == 3) // defeat Safari bug\n",
       "        targ = targ.parentNode;\n",
       "\n",
       "    // jQuery normalizes the pageX and pageY\n",
       "    // pageX,Y are the mouse positions relative to the document\n",
       "    // offset() returns the position of the element relative to the document\n",
       "    var x = e.pageX - $(targ).offset().left;\n",
       "    var y = e.pageY - $(targ).offset().top;\n",
       "\n",
       "    return {\"x\": x, \"y\": y};\n",
       "};\n",
       "\n",
       "/*\n",
       " * return a copy of an object with only non-object keys\n",
       " * we need this to avoid circular references\n",
       " * http://stackoverflow.com/a/24161582/3208463\n",
       " */\n",
       "function simpleKeys (original) {\n",
       "  return Object.keys(original).reduce(function (obj, key) {\n",
       "    if (typeof original[key] !== 'object')\n",
       "        obj[key] = original[key]\n",
       "    return obj;\n",
       "  }, {});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.mouse_event = function(event, name) {\n",
       "    var canvas_pos = mpl.findpos(event)\n",
       "\n",
       "    if (name === 'button_press')\n",
       "    {\n",
       "        this.canvas.focus();\n",
       "        this.canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    var x = canvas_pos.x * mpl.ratio;\n",
       "    var y = canvas_pos.y * mpl.ratio;\n",
       "\n",
       "    this.send_message(name, {x: x, y: y, button: event.button,\n",
       "                             step: event.step,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "\n",
       "    /* This prevents the web browser from automatically changing to\n",
       "     * the text insertion cursor when the button is pressed.  We want\n",
       "     * to control all of the cursor setting manually through the\n",
       "     * 'cursor' event from matplotlib */\n",
       "    event.preventDefault();\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    // Handle any extra behaviour associated with a key event\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.key_event = function(event, name) {\n",
       "\n",
       "    // Prevent repeat events\n",
       "    if (name == 'key_press')\n",
       "    {\n",
       "        if (event.which === this._key)\n",
       "            return;\n",
       "        else\n",
       "            this._key = event.which;\n",
       "    }\n",
       "    if (name == 'key_release')\n",
       "        this._key = null;\n",
       "\n",
       "    var value = '';\n",
       "    if (event.ctrlKey && event.which != 17)\n",
       "        value += \"ctrl+\";\n",
       "    if (event.altKey && event.which != 18)\n",
       "        value += \"alt+\";\n",
       "    if (event.shiftKey && event.which != 16)\n",
       "        value += \"shift+\";\n",
       "\n",
       "    value += 'k';\n",
       "    value += event.which.toString();\n",
       "\n",
       "    this._key_event_extra(event, name);\n",
       "\n",
       "    this.send_message(name, {key: value,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onclick = function(name) {\n",
       "    if (name == 'download') {\n",
       "        this.handle_save(this, null);\n",
       "    } else {\n",
       "        this.send_message(\"toolbar_button\", {name: name});\n",
       "    }\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {\n",
       "    this.message.textContent = tooltip;\n",
       "};\n",
       "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to  previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Pan axes with left mouse, zoom with right\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n",
       "\n",
       "mpl.extensions = [\"eps\", \"jpeg\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n",
       "\n",
       "mpl.default_extension = \"png\";var comm_websocket_adapter = function(comm) {\n",
       "    // Create a \"websocket\"-like object which calls the given IPython comm\n",
       "    // object with the appropriate methods. Currently this is a non binary\n",
       "    // socket, so there is still some room for performance tuning.\n",
       "    var ws = {};\n",
       "\n",
       "    ws.close = function() {\n",
       "        comm.close()\n",
       "    };\n",
       "    ws.send = function(m) {\n",
       "        //console.log('sending', m);\n",
       "        comm.send(m);\n",
       "    };\n",
       "    // Register the callback with on_msg.\n",
       "    comm.on_msg(function(msg) {\n",
       "        //console.log('receiving', msg['content']['data'], msg);\n",
       "        // Pass the mpl event to the overridden (by mpl) onmessage function.\n",
       "        ws.onmessage(msg['content']['data'])\n",
       "    });\n",
       "    return ws;\n",
       "}\n",
       "\n",
       "mpl.mpl_figure_comm = function(comm, msg) {\n",
       "    // This is the function which gets called when the mpl process\n",
       "    // starts-up an IPython Comm through the \"matplotlib\" channel.\n",
       "\n",
       "    var id = msg.content.data.id;\n",
       "    // Get hold of the div created by the display call when the Comm\n",
       "    // socket was opened in Python.\n",
       "    var element = $(\"#\" + id);\n",
       "    var ws_proxy = comm_websocket_adapter(comm)\n",
       "\n",
       "    function ondownload(figure, format) {\n",
       "        window.open(figure.imageObj.src);\n",
       "    }\n",
       "\n",
       "    var fig = new mpl.figure(id, ws_proxy,\n",
       "                           ondownload,\n",
       "                           element.get(0));\n",
       "\n",
       "    // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n",
       "    // web socket which is closed, not our websocket->open comm proxy.\n",
       "    ws_proxy.onopen();\n",
       "\n",
       "    fig.parent_element = element.get(0);\n",
       "    fig.cell_info = mpl.find_output_cell(\"<div id='\" + id + \"'></div>\");\n",
       "    if (!fig.cell_info) {\n",
       "        console.error(\"Failed to find cell for figure\", id, fig);\n",
       "        return;\n",
       "    }\n",
       "\n",
       "    var output_index = fig.cell_info[2]\n",
       "    var cell = fig.cell_info[0];\n",
       "\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.handle_close = function(fig, msg) {\n",
       "    var width = fig.canvas.width/mpl.ratio\n",
       "    fig.root.unbind('remove')\n",
       "\n",
       "    // Update the output cell to use the data from the current canvas.\n",
       "    fig.push_to_output();\n",
       "    var dataURL = fig.canvas.toDataURL();\n",
       "    // Re-enable the keyboard manager in IPython - without this line, in FF,\n",
       "    // the notebook keyboard shortcuts fail.\n",
       "    IPython.keyboard_manager.enable()\n",
       "    $(fig.parent_element).html('<img src=\"' + dataURL + '\" width=\"' + width + '\">');\n",
       "    fig.close_ws(fig, msg);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.close_ws = function(fig, msg){\n",
       "    fig.send_message('closing', msg);\n",
       "    // fig.ws.close()\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.push_to_output = function(remove_interactive) {\n",
       "    // Turn the data on the canvas into data in the output cell.\n",
       "    var width = this.canvas.width/mpl.ratio\n",
       "    var dataURL = this.canvas.toDataURL();\n",
       "    this.cell_info[1]['text/html'] = '<img src=\"' + dataURL + '\" width=\"' + width + '\">';\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Tell IPython that the notebook contents must change.\n",
       "    IPython.notebook.set_dirty(true);\n",
       "    this.send_message(\"ack\", {});\n",
       "    var fig = this;\n",
       "    // Wait a second, then push the new image to the DOM so\n",
       "    // that it is saved nicely (might be nice to debounce this).\n",
       "    setTimeout(function () { fig.push_to_output() }, 1000);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>')\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items){\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) { continue; };\n",
       "\n",
       "        var button = $('<button class=\"btn btn-default\" href=\"#\" title=\"' + name + '\"><i class=\"fa ' + image + ' fa-lg\"></i></button>');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    // Add the status bar.\n",
       "    var status_bar = $('<span class=\"mpl-message\" style=\"text-align:right; float: right;\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "\n",
       "    // Add the close button to the window.\n",
       "    var buttongrp = $('<div class=\"btn-group inline pull-right\"></div>');\n",
       "    var button = $('<button class=\"btn btn-mini btn-primary\" href=\"#\" title=\"Stop Interaction\"><i class=\"fa fa-power-off icon-remove icon-large\"></i></button>');\n",
       "    button.click(function (evt) { fig.handle_close(fig, {}); } );\n",
       "    button.mouseover('Stop Interaction', toolbar_mouse_event);\n",
       "    buttongrp.append(button);\n",
       "    var titlebar = this.root.find($('.ui-dialog-titlebar'));\n",
       "    titlebar.prepend(buttongrp);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(el){\n",
       "    var fig = this\n",
       "    el.on(\"remove\", function(){\n",
       "\tfig.close_ws(fig, {});\n",
       "    });\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(el){\n",
       "    // this is important to make the div 'focusable\n",
       "    el.attr('tabindex', 0)\n",
       "    // reach out to IPython and tell the keyboard manager to turn it's self\n",
       "    // off when our div gets focus\n",
       "\n",
       "    // location in version 3\n",
       "    if (IPython.notebook.keyboard_manager) {\n",
       "        IPython.notebook.keyboard_manager.register_events(el);\n",
       "    }\n",
       "    else {\n",
       "        // location in version 2\n",
       "        IPython.keyboard_manager.register_events(el);\n",
       "    }\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    var manager = IPython.notebook.keyboard_manager;\n",
       "    if (!manager)\n",
       "        manager = IPython.keyboard_manager;\n",
       "\n",
       "    // Check for shift+enter\n",
       "    if (event.shiftKey && event.which == 13) {\n",
       "        this.canvas_div.blur();\n",
       "        event.shiftKey = false;\n",
       "        // Send a \"J\" for go to next cell\n",
       "        event.which = 74;\n",
       "        event.keyCode = 74;\n",
       "        manager.command_mode();\n",
       "        manager.handle_keydown(event);\n",
       "    }\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    fig.ondownload(fig, null);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.find_output_cell = function(html_output) {\n",
       "    // Return the cell and output element which can be found *uniquely* in the notebook.\n",
       "    // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n",
       "    // IPython event is triggered only after the cells have been serialised, which for\n",
       "    // our purposes (turning an active figure into a static one), is too late.\n",
       "    var cells = IPython.notebook.get_cells();\n",
       "    var ncells = cells.length;\n",
       "    for (var i=0; i<ncells; i++) {\n",
       "        var cell = cells[i];\n",
       "        if (cell.cell_type === 'code'){\n",
       "            for (var j=0; j<cell.output_area.outputs.length; j++) {\n",
       "                var data = cell.output_area.outputs[j];\n",
       "                if (data.data) {\n",
       "                    // IPython >= 3 moved mimebundle to data attribute of output\n",
       "                    data = data.data;\n",
       "                }\n",
       "                if (data['text/html'] == html_output) {\n",
       "                    return [cell, data, j];\n",
       "                }\n",
       "            }\n",
       "        }\n",
       "    }\n",
       "}\n",
       "\n",
       "// Register the function which deals with the matplotlib target/channel.\n",
       "// The kernel may be null if the page has been refreshed.\n",
       "if (IPython.notebook.kernel != null) {\n",
       "    IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAgAElEQVR4nO3df3QV9Z3/8ZGQewNqclUIKZgbsQVqyWJ7ll5K3SXbTdtVYJOzXUiXnkJWW34ssFj5obQF0p/ChlTOHkD7Q4H9rkI2RSEWdO22ZXtsEKzblF9BaCnYCIKtbVKBJMDN6/sHx1suucEkM3Pnc2eej3Pef2TuJIwz70/eL5PMHUsAAAAIFMvrAwAAAEB6EQABAAAChgAIAAAQMARAAACAgCEAAgAABAwBEAAAIGAIgAAAAAFDAAQAAAgYAiAAAEDAEAABAAAChgAIAAAQMARAAACAgCEAAgAABAwBEAAAIGAIgAAAAAFDAAQAAAgYAiAAAEDAEAABAAAChgAIAAAQMARAAACAgCEAAgAABAwBEAAAIGAIgAAAAAFDAAQAAAgYAiAAAEDAEAABAAAChgAIAAAQMARAAACAgCEAAgAABAwBEAAAIGAIgAAAAAFDAAQAAAgYAiAAAEDAEAABAAAChgAIAAAQMARAAACAgCEAAgAABAwBEAAAIGAIgAAAAAFDAAQAAAgYAiAAAEDAEAABAAAChgAIAAAQMARAAACAgCEAAgAABAwBEAAAIGAIgAAAAAFDAAQAAAgYAiAAAEDAEAABAAAChgAIAAAQMARAAACAgCEAAgAABAwBEAAAIGAIgAAAAAFDAAQAAAgYAqAN8Xhczc3NamlpUWtrK0VRFEVRGVAtLS1qbm5WPB73Okp4hgBoQ3NzsyzLoiiKoigqA6u5udnrKOEZAqANLS0tiQby+v9m/FinTp1KLNJTp055fjx+Ls61ucW1oSjn650f4LS0tHgdJTxDALShtbVVlmWptbXV60PxpbNnzyYG39mzZ70+HF/jXJuLawM4j/lNALSFBnIXgy99ONfm4toAzmN+EwBtoYHc1d7ersrKSlVWVqq9vd3rw/E1zrW5uDaA85jfBEBbaCAAADIP85sAaAsNBABA5mF+EwBtoYHc1dnZqbNnz+rs2bPq7Oz0+nB8jXNtLq4N4DzmNwHQFhrIXfzxe/pwrs3FtQGcx/wmANpCA7mLwZc+nGtzcW0A5zG/CYC20EDuYvClD+faXFwbwHnMbwKgLTSQuxh86cO5NhfXBnAe85sAaAsN5C4GX/pwrs3FtQGcx/wmANpCA7mLwZc+nGtzcW0A5zG/CYC20EDuYvClD+faXFwbwHnMb4MC4Lp161RUVKRwOKxYLKa9e/dec/+6ujqNGjVK4XBYxcXF2rlzZ9LrnZ2dWr58uQoKCpSTk6PS0lIdPXo08fquXbsS31SvrpdffrlHx0wDuautrU1TpkzRlClT1NbW5vXh+Brn2lxcG8B5zG9DAmBtba1CoZA2bNigQ4cOaebMmYpEIjpz5kzK/RsaGpSVlaXq6mo1NTVp2bJlys7O1oEDBxL7rFq1Snl5edq+fbv27dunsrIyDR8+PPENtKOjQ2+88UZSff7zn9fw4cN7/GarNBAAAJmH+W1IAIzFYpo3b17i43g8rqFDh2rlypUp96+oqNCkSZOSto0bN06zZ8+WdPmnfwUFBVq9enXi9ZaWFoXDYW3ZsiXl17xw4YIGDx6sr33taz0+bhoIJit6aMe7FgAEEfPbgADY0dGhrKwsbdu2LWn7jBkzVFZWlvJzCgsLtWbNmqRtK1as0JgxYyRJx44dk2VZamxsTNpnwoQJWrBgQcqvuXXrVvXr10/Nzc09PnYaCCYjAAJAasxvAwLgyZMnZVmWdu/enbR9yZIlisViKT8nOztbmzdvTtq2fv165efnS7r8K2LLsnTq1KmkfaZOnaqKioqUX/Oee+7RPffcc81jbW9vV2tra6Kam5sD30Bu4o/f7elNAORcm4trAziPAEgAlCQ1NzerX79+2rp16zWPtaqqKuVNI0FuIDcx+OwhAPoD1wZwHgHQgABowq+Av/a1r2nw4MG6cOHCNY+VnwCmF4PPnp4EwHeq8IGtiXNd+MBWfk1sENYB4DwCoAEBULp8E8j8+fMTH8fjcQ0bNuyaN4FMnjw5adv48eO73ARSU1OTeL21tTXlTSCdnZ0aPny4Fi1a1OvjpoHcxeCzhwDoD6wDwHnMb0MCYG1trcLhsDZt2qSmpibNmjVLkUhEp0+fliRNnz5dS5cuTezf0NCg/v37q6amRocPH1ZVVVXKt4GJRCKqr6/X/v37VV5envQ2MO/40Y9+JMuydPjw4V4fNw3kLgafPQRAf2AdAM5jfhsSACVp7dq1ikajCoVCisVi2rNnT+K1kpISVVZWJu1fV1enkSNHKhQKafTo0d2+EfSQIUMUDodVWlqqI0eOdPl3p02bpo9+9KN9OmYayF0MPnsIgP7AOgCcx/w2KABmIhrIXQw+ewiA/sA6AJzH/CYA2kIDuautrU0TJ07UxIkTeQRWH/QmAEYXPaMBt4/VgNvHKrroGQKgQVgHgPOY3wRAW2ggmKw3AZA3iwYQJMxvAqAtNBBMRgAEgNSY3wRAW2ggmIwACACpMb8JgLbQQO46e/asBg4cqIEDB/LH733Q25tArssO67rsMDeBGIZ1ADiP+U0AtIUGchd3P9rDXcD+wDoAnMf8JgDaQgO5i8HXPad+vUsANB/rAHAe85sAaAsN5C4GX/cIgMHBOgCcx/wmANpCA7mLwdc9AmBwsA4A5zG/CYC20EDuYvB1jwAYHKwDwHnMbwKgLTSQuxh83SMABgfrAHAe85sAaAsN5K7z58+rpKREJSUlOn/+vNeHYxTHA+DCpxUuLFa4sFiFC58mABqEdQA4j/lNALSFBoJXnA6A1yoA8BvmNwHQFhoIXiEAAkDfMb8JgLbQQPCKaQGQIAkgkzC/CYC20EDuOnv2rAYNGqRBgwbxx+9XcfxvAB/Yqn4DctVvQG6fbgIhALqHdQA4j/lNALSFBnIXdz92z40AaOcuYAKge1gHgPOY3wRAW2ggdzH4ukcADA7WAeA85jcB0BYayF0Mvu4RAIODdQA4j/lNALSFBnIXg697BMDgYB0AzmN+EwBtoYHcxeDrXjoDoFOFvmEdAM5jfhMAbaGB3MXg6x4BMDhYB4DzmN8EQFtoIHedP39eY8eO1dixY3kE1lUcD4ALn1aoYIRCBSO6PAqOAOgt1gHgPOY3AdAWGghecSOguV0AYArmNwHQFhoIXvE6zLkVAAmSANKB+U0AtIUGgle8DnMEQACZjPlNALSFBnLXuXPnVFRUpKKiIp07d87rwzGK0+GscOFWZeXmKys3X4ULvbsJhADYFesAcB7zmwBoCw3kLu5+7J7jAdCQu4AJgF2xDgDnMb8JgLbQQO5i8HWPABgcrAPAecxvAqAtNJC7GHzdIwAGB+sAcB7zmwBoCw3kLgZf9wiAwcE6AJzH/CYA2kIDuYvB1z0CYHCwDgDnMb8JgLbQQO5i8HWPABgcrAPAecxvAqAtNJC7zp07pw984AP6wAc+wNtfXMXxALhwq7JviSr7lihvA2MY1gHgPOY3AdAWGghecSOguV1O/XcBgF3MbwKgLTQQvOJ1mCMAAshkzG8CoC00ELzidZgjAALIZMxvAqAtNJC7+Nun7jkdzvgbQHOxDgDnMb8JgLbQQO7i7sfuOR4AuQvYWKwDwHnMb4MC4Lp161RUVKRwOKxYLKa9e/dec/+6ujqNGjVK4XBYxcXF2rlzZ9LrnZ2dWr58uQoKCpSTk6PS0lIdPXq0y9fZsWOHYrGYcnJyFIlEVF5e3uNjpoHcxeDrXiYGwHQGST9hHQDOY34bEgBra2sVCoW0YcMGHTp0SDNnzlQkEtGZM2dS7t/Q0KCsrCxVV1erqalJy5YtU3Z2tg4cOJDYZ9WqVcrLy9P27du1b98+lZWVafjw4Wpra0vss3XrVt1000167LHHdOTIER06dEj/9V//1ePjpoHcxeDrHgEwOFgHgPOY34YEwFgspnnz5iU+jsfjGjp0qFauXJly/4qKCk2aNClp27hx4zR79mxJl3/6V1BQoNWrVydeb2lpUTgc1pYtWyRJFy9e1LBhw/T444/3+bhpIHcx+LpHAAwO1gHgPOa3AQGwo6NDWVlZ2rZtW9L2GTNmqKysLOXnFBYWas2aNUnbVqxYoTFjxkiSjh07Jsuy1NjYmLTPhAkTtGDBAknS3r17ZVmWNmzYoA9+8IMqKCjQ3XffnfRTxHdDA7mLwdc9AmBwsA4A5zG/DQiAJ0+elGVZ2r17d9L2JUuWKBaLpfyc7Oxsbd68OWnb+vXrlZ+fL+nyr4gty9KpU6eS9pk6daoqKiokSVu2bJFlWYpGo9q6dateeeUVTZs2TbfccoveeuutlP9ue3u7WltbE9Xc3Bz4BnITg697BMDgYB0AziMABjgAPvXUU7IsS9/5zncSr7e3t2vQoEH69re/nfLfraqqSnwjvrKC3EBuOnfunIqKilRUVMTbX1zF8QC4cKuycvOVlZvv2tvAEAD7hnUAOI8AaEAA9OpXwD/5yU9kWZZefPHFpH1isZi+9KUvpfx3+QkgTOF1CCMAAshkBEADAqB0OXTNnz8/8XE8HtewYcOueRPI5MmTk7aNHz++y00gNTU1iddbW1uTbgJ55+MrbwK5cOGC8vPzk34qeC00EHrLqYDjdQgjAALIZMxvQwJgbW2twuGwNm3apKamJs2aNUuRSESnT5+WJE2fPl1Lly5N7N/Q0KD+/furpqZGhw8fVlVVVcq3gYlEIqqvr9f+/ftVXl7e5W1g7r//fg0bNkwvvPCCXn31VX3uc59Tfn6+/vCHP/TouGkg9BYBkAAIwHvMb0MCoCStXbtW0WhUoVBIsVhMe/bsSbxWUlKiysrKpP3r6uo0cuRIhUIhjR49uts3gh4yZIjC4bBKS0t15MiRpH0uXLigRYsWKT8/XzfeeKM+/vGP6+DBgz0+ZhrIXefPn9fYsWM1duxYnT9/3uvDcYTX4am7Klz4tEIFIxQqGKHChU97fjwEwD/z4zoAvMb8NigAZiIayF1+vPvR6/DUbQDkLmBj+XEdAF5jfhMAbaGB3OXHwed1eCIAZh4/rgPAa8xvAqAtNJC7/Dj4vA5PBMDM48d1AHiN+U0AtIUGcpcfB5/X4YkAmHn8uA4ArzG/CYC20EDu8uPg8zo8EQAzjx/XAeA15jcB0BYayF1+HHxehycCYObx4zoAvMb8JgDaQgO56+zZsxo0aJAGDRrkm8HndXi6VgDsNyBX/QbkEgAN48d1AHiN+U0AtIUGQm95HZ78UABgF/ObAGgLDYTe8jo8+aEAwC7mNwHQFhoIveV1ePJDAYBdzG8CoC00kLvOnz+vkpISlZSU+OYRWF6Hp+6qcOHTChcWK1xYzKPgDOPHdQB4jflNALSFBnKXH+9+9Do8dRsAuQvYWH5cB4DXmN8EQFtoIHf5cfB5HZ4IgJnHj+sA8BrzmwBoCw3kLj8OPq/DEwEw8/hxHQBeY34TAG2hgdzlx8HndXgiAGYeP64DwGvMbwKgLTSQu/w4+LwOTwTAzOPHdQB4jflNALSFBnKXHwef1+GJAJh5/LgOAK8xvwmAttBA7jp79qwGDhyogQMH+mbweR2erhUAr8sO67rsMAHQMH5cB4DXmN8EQFtoIPSW1+HJDwUAdjG/CYC20EDoLa/Dkx8KAOxifhMAbaGB0Ftehyc/FADYxfwmANpCA7mrra1NEydO1MSJE9XW1ub14TjC6/DUXUUXPaMBt4/VgNvHKrroGc+PhwD4Z35cB4DXmN8EQFtoIHf58e5Hr8NTd8VdwOby4zoAvMb8JgDaQgO5y4+Dz+vwRADMPH5cB4DXmN8EQFtoIHdl2uDzOhgRAP0p09YBkAmY3wRAW2ggd2Xa4PM6GBEA/SnT1gGQCZjfBEBbaCB3Zdrg8zoYEQD9KdPWAZAJmN8EQFtoIHdl2uDzOhgRAP0p09YBkAmY3wRAW2ggd2Xa4PM6GBEA/SnT1gGQCZjfBEBbaCBcyetgFJQCALuY3wRAW2ggXMnrYEQREgH0DPObAGgLDYQreR16KAIggJ5hfhMAbaGB3NXW1qYpU6ZoypQpGfEILK9Dj52KLnpGA0fdpYGj7jL+UXBBC4CZtg6ATMD8JgDaQgO5K9P++N3r0GOnMukmkKAFwExbB0AmYH4TAG2hgdyVaYPP69BDACQAAugZ5jcB0BYayF2ZNvi8Dj0EQAIggJ5hfhMAbaGB3JVpg8/r0EMAJAAC6BnmNwHQFhrIXZk2+LwOPQRAAiCAnmF+EwBtoYHclWmDz+vQQwAkAALoGeY3AdAWGshdmTb4vA49BEACIICeYX4bFADXrVunoqIihcNhxWIx7d2795r719XVadSoUQqHwyouLtbOnTuTXu/s7NTy5ctVUFCgnJwclZaW6ujRo0n7FBUVJb6xvlMrV67s8THTQO7q7OzU2bNndfbsWXV2dnp9OO/K69Bjp6IP/kCFD2xV4QNbFX3wB54fDwHwzzJtHQCZgPltSACsra1VKBTShg0bdOjQIc2cOVORSERnzpxJuX9DQ4OysrJUXV2tpqYmLVu2TNnZ2Tpw4EBin1WrVikvL0/bt2/Xvn37VFZWpuHDhye9kWpRUZG+9rWv6Y033khUb/4PmwbClbwOPZQ/AyAA5zG/DQmAsVhM8+bNS3wcj8c1dOjQbn8aV1FRoUmTJiVtGzdunGbPni3p8v8xFxQUaPXq1YnXW1paFA6HtWXLlsS2oqIirVmzps/HTQPhSl6HHooACKBnmN8GBMCOjg5lZWVp27ZtSdtnzJihsrKylJ9TWFjYJbitWLFCY8aMkSQdO3ZMlmWpsbExaZ8JEyZowYIFiY+Lioo0ZMgQ3XzzzfrgBz+o6upqXbx4scfHTgO5q729XZWVlaqsrFR7e7vXh/OuvA49diq6aJuuLy7V9cWlii7a5vnxEAD/LNPWAZAJmN8GBMCTJ0/Ksizt3r07afuSJUsUi8VSfk52drY2b96ctG39+vXKz8+XdPlXxJZl6dSpU0n7TJ06VRUVFYmPv/Wtb2nXrl3at2+fHnvsMUUiET3wwAPdHmt7e7taW1sT1dzcHPgGclOm/fG716HHTnETiLkybR0AmYAAGPAAeLUnnnhC/fv37/b/squqqrrcNBL0BnJTpg0+r0MPAZAACKBnCIAGBEAvfwV8tYMHD8qyLL366qspX+cngOmVaYPP69BDACQAAugZAqABAVC6fBPI/PnzEx/H43ENGzbsmjeBTJ48OWnb+PHju9wEUlNTk3i9tbW1y00gV3vyySfVr18//eEPf+jRcdNA7sq0wed16CEAEgAB9Azz25AAWFtbq3A4rE2bNqmpqUmzZs1SJBLR6dOnJUnTp0/X0qVLE/s3NDSof//+qqmp0eHDh1VVVZXybWAikYjq6+u1f/9+lZeXJ70NzO7du7VmzRr98pe/1LFjx/Tkk09q8ODBmjFjRo+PmwZyl0mDz+tAQwAkAJqwDgC/YH4bEgAlae3atYpGowqFQorFYtqzZ0/itZKSElVWVibtX1dXp5EjRyoUCmn06NHdvhH0kCFDFA6HVVpaqiNHjiRe/7//+z+NGzdOeXl5ysnJ0R133KGHH364V3fZ0UDuMmnweR1oCIAEQBPWAeAXzG+DAmAmooHcZdLg8zrQEAAJgCasA8AvmN8EQFtoIHd1dnbqzTff1Jtvvun5I7C8DjRuV/TBH+jWf31Kt/7rUzwKzjAmrQPAL5jfBEBbaKDg8DrQUMENgACcx/wmANpCAwWH14GGIgACcA7zmwBoCw3krvb2ds2dO1dz5871/BFYXgcatyu6aJtu+NAk3fChSTwKzjAmrQPAL5jfBEBbaCB3mfTH714HGreLm0DMZdI6APyC+U0AtIUGcpdJg8/rQEMAJACasA4Av2B+EwBtoYHcZdLg8zrQEAAJgCasA8AvmN8EQFtoIHeZNPi8DjQEQAKgCesA8AvmNwHQFhrIXSYNPq8DDQGQAGjCOgD8gvlNALSFBnKXSYPP60BDACQAmrAOAL9gfhMAbaGB3GXS4PM60BAACYAmrAPAL5jfBEBbaCB3xeNxHT9+XMePH1c8Hvf0WLwONG5X9MFnNWzOExo25wlFH3zW8+MhAP6ZSesA8AvmNwHQFhooOLwONFRwAyAA5zG/CYC20EDB4XWgoQiAAJzD/CYA2kIDuaujo0OLFy/W4sWL1dHR4emxeB1o3K7o4m3KjX1KubFPKbqYR8GZxKR1APgF85sAaAsN5C6T/vjd60DjdnETiLlMWgeAXzC/CYC20EDuMmnweR1oCIAEQBPWAeAXzG8CoC00kLtMGnxeBxoCIAHQhHUA+AXzmwBoCw3kLpMGn9eBhgBIADRhHQB+wfwmANpCA7nLpMHndaAhABIATVgHgF8wvwmAttBA7jJp8HkdaAiABEAT1gHgF8xvAqAtNJC7TBp8XgcaAmBwQ6JJ6wDwC+Y3AdAWGshd8XhcBw8e1MGDBz1/BJbXYcXtij74rN5z33q95771vngUnJ8CoEnrAPAL5jcB0BYaKDi8DitUcAMgAOcxvwmAttBAweF1WKEIgACcw/wmANpCA7mro6NDVVVVqqqq8vwRWF6HFbcrunib8u6apry7pvniUXB+CoAmrQPAL5jfBEBbaCB3mfTH716HFbeLm0DMZdI6APyC+U0AtIUGcpdJg8/rsEIAJACasA4Av2B+EwBtoYHcZdLg8zqsEAAJgCasA8AvmN8EQFtoIHela/B5HURMKAKguQiAgPOY3wRAW2ggdxEACYAEQAIg4AbmNwHQFhrIXQRAAiABkAAIuIH5TQC0hQZyFwGQAEgAJAACbmB+EwBtoYHcdenSJb388st6+eWXdenSJdf+Ha+DiAkVXVKvghmPqGDGI4ouqff8eAiAf5audQAECfObAGgLDeQPXgcRigAIIL2Y3wRAW2ggf/A6iFAEQADpxfwmANpCA7mro6ND1dXVqq6udvURWF4HERMqunibIn9zryJ/cy+PgjNMutYBECTMbwKgLTSQu7gJJH3FTSDm4iYQwHnMbwKgLTSQuwiABEACIAEQcAPz26AAuG7dOhUVFSkcDisWi2nv3r3X3L+urk6jRo1SOBxWcXGxdu7cmfR6Z2enli9froKCAuXk5Ki0tFRHjx5N+bXa29t15513yrIsNTY29viYaSB3EQAJgARAAiDgBua3IQGwtrZWoVBIGzZs0KFDhzRz5kxFIhGdOXMm5f4NDQ3KyspSdXW1mpqatGzZMmVnZ+vAgQOJfVatWqW8vDxt375d+/btU1lZmYYPH662trYuX2/BggW65557CICGIQASAAmABEDADcxvQwJgLBbTvHnzEh/H43ENHTpUK1euTLl/RUWFJk2alLRt3Lhxmj17tqTLP/0rKCjQ6tWrE6+3tLQoHA5ry5YtSZ/33HPP6f3vf78OHTpEADQMAZAASAAkAAJuYH4bEAA7OjqUlZWlbdu2JW2fMWOGysrKUn5OYWGh1qxZk7RtxYoVGjNmjCTp2LFjKcPchAkTtGDBgsTHp0+f1rBhw/Tzn/9cx48ff9cA2N7ertbW1kQ1NzcHvoHcRAAkABIACYCAGwiABgTAkydPyrIs7d69O2n7kiVLFIvFUn5Odna2Nm/enLRt/fr1ys/Pl3T5V8SWZenUqVNJ+0ydOlUVFRWSLv+U8O6779bXv/51SepRAKyqqkp8I76ygtxAbiIAEgAJgARAwA0EwAAHwH//93/XXXfdlXi0Ej8BNM+lS5e0a9cu7dq1i0fBuVzRJfUaMu1hDZn2MI+CM0y61gEQJARAAwKgV78CLi8vV79+/ZSVlZUoy7KUlZWlGTNm9OjYaSB/8DqIUARAAOnF/DYgAEqXbwKZP39+4uN4PK5hw4Zd8yaQyZMnJ20bP358l5tAampqEq+3trYm3QTy2muv6cCBA4l64YUXZFmWtm7dqubm5h4dNw3kD14HEYoACCC9mN+GBMDa2lqFw2Ft2rRJTU1NmjVrliKRiE6fPi1Jmj59upYuXZrYv6GhQf3791dNTY0OHz6sqqqqlG8DE4lEVF9fr/3796u8vLzbt4GRevYr4KvRQO66cOGC1q1bp3Xr1unChQuu/TteBxETKrp4u27+xBzd/Ik5ii7e7vnxEAD/LF3rAAgS5rchAVCS1q5dq2g0qlAopFgspj179iReKykpUWVlZdL+dXV1GjlypEKhkEaPHt3tG0EPGTJE4XBYpaWlOnLkSLf/PgHQPNwEkr4K4k0gPSkTcBMI4Dzmt0EBMBPRQO4iAKavCICpywQEQMB5zG8CoC00kLsIgOkrAmDqMgEBEHAe85sAaAsN5C4CYPqKAJi6TEAABJzH/CYA2kIDuYsAmL4iAKYuExAAAecxvwmAttBA7iIApq8IgKnLBARAwHnMbwKgLTSQuwiA6SsCYOoyAQEQcB7zmwBoCw3krosXL2rHjh3asWOHLl686Nq/43XIMKGiS+o1eEqVBk+pCsyj4HpSJkjXOgCChPlNALSFBvIHr0MGZW4B8CfmNwHQFhrIH7wOGZS5BcCfmN8EQFtoIHdduHBBGzdu1MaNG3kUnMsVXbxdt0z8gm6Z+IXAPAquJ2WCdK0DIEiY3wRAW2ggdznxx+9eB4hMKW4CSV0m4CYQwHnMbwKgLTSQuwiA6SsCYOoyAQEQcB7zmwBoCw3kLgJg+ooAmLpMQAAEnMf8JgDaQgO5iwCYviIApi4TEAAB5zG/CYC20EDuIgCmrwiAqcsEBEDAecxvAqAtNJC7CIDpKwJg6jIBARBwHvObAGgLDeQuAmD6igCYukxAAAScx/wmANpCA7nr4sWLqqurU11dXZ8fgeV1gMiUii6p16DypRpUvpRHwV1RJnBiHQBIxvwmANpCA5nP6wBBZXYB8CfmNwHQFhrIfF4HCCqzC4A/Mb8JgLbQQO7iV8DpK34FnLpMwK+AAecxvwmAttBA7uImkPQVN4GkLhNwEwjgPOY3AdAWGshdBMD0FQEwdZmAAAg4j5/GSQcAACAASURBVPlNALSFBnIXATB9RQBMXSYgAALOY34TAG2hgdxFAExfEQBTlwkIgIDzmN8EQFtoIHe92+DzOhz4qQiAqcsEBEDAecxvAqAtNJC7CIDpKwJg6jIBARBwHvObAGgLDeQuAmD6igCYukxAAAScx/wmANpCA7nrwoUL2rhxozZu3KgLFy50ed3rcOCnii7erlsmfkG3TPyCoou3e348ppQJ3m0dAOg95jcB0BYayFtehwPK/wXAn5jfBEBbaCBveR0OKP8XAH9ifhMAbaGB3HXx4kXt2LFDO3bsSPkILK/DgZ8quqReg6dUafCUKh4Fd0WZ4N3WAYDeY34TAG2hgdzFTSDpK24CSV0m4CYQwHnMbwKgLTSQuwiA6SsCYOoyAQEQcB7zmwBoCw3kLgJg+ooAmLpMQAAEnMf8JgDaQgO5iwCYviIApi4TEAAB5zG/CYC20EDuIgCmrwiAqcsEBEDAecxvAqAtNJC7CIDpKwJg6jIBARBwHvObAGgLDeQuAmD6igCYukxAAAScx/wmANpCA7nrwoULWrdundatW8ej4Fyu6OLtuvkTc3TzJ+bwKLgrygTvtg4A9B7z26AAuG7dOhUVFSkcDisWi2nv3r3X3L+urk6jRo1SOBxWcXGxdu7cmfR6Z2enli9froKCAuXk5Ki0tFRHjx5N2ufv//7vVVhYqHA4rIKCAn32s5/VyZMne3zMNJC3vA4HlP8LgD8xvw0JgLW1tQqFQtqwYYMOHTqkmTNnKhKJ6MyZMyn3b2hoUFZWlqqrq9XU1KRly5YpOztbBw4cSOyzatUq5eXlafv27dq3b5/Kyso0fPhwtbW1JfZ55JFH9NJLL+nEiRNqaGjQ+PHjNX78+B4fNw3kLa/DAeX/AuBPzG9DAmAsFtO8efMSH8fjcQ0dOlQrV65MuX9FRYUmTZqUtG3cuHGaPXu2pMs//SsoKNDq1asTr7e0tCgcDmvLli3dHkd9fb2uu+66Hv+ahQZy16VLl7Rr1y7t2rVLly5d6vK61+HATxVdUq8h0x7WkGkP8yi4K8oE77YOAPQe89uAANjR0aGsrCxt27YtafuMGTNUVlaW8nMKCwu1Zs2apG0rVqzQmDFjJEnHjh2TZVlqbGxM2mfChAlasGBByq/51ltvqaKiQnfddVe3x9re3q7W1tZENTc3B76B3MRNIOkrbgJJXSbgJhDAeQRAAwLgyZMnZVmWdu/enbR9yZIlisViKT8nOztbmzdvTtq2fv165efnS7r8K2LLsnTq1KmkfaZOnaqKioqkbQ8++KAGDhwoy7L0kY98RL///e+7PdaqqqrEN+IrK8gN5CYCYPqKAJi6TEAABJxHACQA6ne/+52OHDmiH/7wh7rrrrs0ceJEdXZ2pvx3+QlgehEA01cEwNRlAgIg4DwCoAEB0JRfAUtKBLqrw2h3aCB3EQDTVwTA1GUCAiDgPOa3AQFQunwTyPz58xMfx+NxDRs27Jo3gUyePDlp2/jx47vcBFJTU5N4vbW19V1vAnnttddkWZZ27drVo+OmgdxFAExfEQBTlwkIgIDzmN+GBMDa2lqFw2Ft2rRJTU1NmjVrliKRiE6fPi1Jmj59upYuXZrYv6GhQf3791dNTY0OHz6sqqqqlG8DE4lEVF9fr/3796u8vDzpbWD27NmjtWvXqrGxUSdOnNCPf/xjffSjH9V73/tetbe39+i4aSB3EQDTVwTAvpfX6wBA7zG/DQmAkrR27VpFo1GFQiHFYjHt2bMn8VpJSYkqKyuT9q+rq9PIkSMVCoU0evTobt8IesiQIQqHwyotLdWRI0cSr+/fv18f+9jHdPPNNyscDuu2227TnDlz9Prrr/f4mGkgdxEA01cEwL6X1+sAQO8xvw0KgJmIBnJXR0eHqqurVV1drY6Oji6vez34/VTRxdsU+Zt7FfmbexVdvM3z48mk8nodAOg95jcB0BYayFteD36KKnrIjL8TBNA7zG8CoC00kLe8HvwUVfQQARDIRMxvAqAtNJC7Ll26pJdfflkvv/wyj4JzuaJL6lUw4xEVzHiER8H1srxeBwB6j/lNALSFBnIXN4Gkr7gJpO/l9ToA0HvMbwKgLTSQuwiA6SsCYN/L63UAoPeY3wRAW2ggdxEA01cEwL6X1+sAQO8xvwmAttBA7iIApq8IgH0vr9cBgN5jfhMAbaGB3EUATF8RAPteXq8DAL3H/CYA2kIDuYsAmL4iAPa9vF4HAHqP+U0AtIUGchcBMH1FAOx7eb0OAPQe85sAaAsN5K6Ojg5VVVWpqqqKR8G5XNHF25R31zTl3TWNR8H1srxeBwB6j/lNALSFBvKW14Ofoooe4kkgQCZifhMAbaGBvOX14KeooocIgEAmYn4TAG2hgdwVj8d18OBBHTx4UPF4vMvrXg9+P1X0wWf1nvvW6z33rVf0wWc9P55MKq/XAYDeY34TAG2hgdxT9BA3JqSzONd9L7dxEwjgPOY3AdAWGsg9hJL0Fue67+U2AiDgPOY3AdAWGsg9hJL0Fue67+U2AiDgPOY3AdAWGsg9hJL0Fue67+U2AiDgPOY3AdAWGsg9hJL0Fue67+U2AiDgPOY3AdAWGsg9hJL0Fue67+U2AiDgPOY3AdAWGsg9hJL0Fue67+U2AiDgPOY3AdAWGsg9RQ9dfjxZbuxTyo19iseTuVyca3fLjo6ODi1evFiLFy/mUXCAQ5jfBEBbaCD3eD2wKcrJAmAW5jcB0BYayD1eD2yKcrIAmIX5TQC0hQZyT9FDlx9PNmzOExo25wkeT+Zyca7dLTvi8biOHz+u48eP8yg4wCHMbwKgLTSQe4oe4saEdBbn2t2yg5tAAOcxvwmAttBA7iGUpLc41+6WHQRAwHnMbwKgLTSQewgl6S3OtbtlBwEQcB7zmwBoCw3kHkJJeotz7W7ZQQAEnMf8JgDaQgO5h1CS3uJcu1t2EAAB5zG/CYC20EDuIZSktzjX7pYdBEDAecxvAqAtNJB7CCXpLc61u2UHARBwHvObAGgLDeSeood2KLpom2740CTd8KFJii7i8WRuFufa3bKjvb1dc+fO1dy5c9Xe3u7QCgOCjflNALSFBnKP1wObopwsAGZhfhMAbaGB3OP1wKYoJwuAWZjfBEBbaCD3FD20Q9EHf6Bb//Up3fqvTyn64A88H+J+Ls61u2VHZ2en3nzzTb355pvq7Ox0aIUBwcb8JgDaQgO5p+ghbkxIZ3Gu3S07uAkEcB7zmwBoCw3kHkJJeotz7W7ZQQAEnMf8JgDaQgO5h1CS3uJcu1t2EAAB5zG/DQqA69atU1FRkcLhsGKxmPbu3XvN/evq6jRq1CiFw2EVFxdr586dSa93dnZq+fLlKigoUE5OjkpLS3X06NHE68ePH9d9992n2267TTk5Obr99tu1YsUKdXR09PiYaSD3EErSW5xrd8sOAiDgPOa3IQGwtrZWoVBIGzZs0KFDhzRz5kxFIhGdOXMm5f4NDQ3KyspSdXW1mpqatGzZMmVnZ+vAgQOJfVatWqW8vDxt375d+/btU1lZmYYPH662tjZJ0vPPP69//ud/1gsvvKBjx46pvr5e+fn5WrRoUY+PmwZyD6EkvcW5drfsIAACzmN+GxIAY7GY5s2bl/g4Ho9r6NChWrlyZcr9KyoqNGnSpKRt48aN0+zZsyVd/ulfQUGBVq9enXi9paVF4XBYW7Zs6fY4qqurNXz48B4fNw3kHkJJeotz7W7ZQQAEnMf8NiAAdnR0KCsrS9u2bUvaPmPGDJWVlaX8nMLCQq1ZsyZp24oVKzRmzBhJ0rFjx2RZlhobG5P2mTBhghYsWNDtsXz5y1/WX/7lX3b7ent7u1pbWxPV3Nwc+AZyC6EkvcW5drfsIAACziMAGhAAT548KcuytHv37qTtS5YsUSwWS/k52dnZ2rx5c9K29evXKz8/X9LlXxFblqVTp04l7TN16lRVVFSk/Jq/+tWvlJubq+9+97vdHmtVVVXiG/GVFeQGckvRQ5cfT3Z9camuLy7l8WQuF+fa3bKjvb1dlZWVqqys5FFwgEMIgARASdLrr7+u9773vfrc5z53zWPlJ4Dp4/XApignC4BZCIAGBECvfwV88uRJjRgxQtOnT1c8Hu/VsdNA7vF6YFOUkwXALMxvAwKgdPkmkPnz5yc+jsfjGjZs2DVvApk8eXLStvHjx3e5CaSmpibxemtra5ebQF5//XWNGDFC//RP/6RLly71+rhpIPcUPXT58WSFD2xV4QNbeTyZy8W5drfs6Ozs1NmzZ3X27FkeBQc4hPltSACsra1VOBzWpk2b1NTUpFmzZikSiej06dOSpOnTp2vp0qWJ/RsaGtS/f3/V1NTo8OHDqqqqSvk2MJFIRPX19dq/f7/Ky8uT3gbm9ddf1/ve9z6Vlpbq9ddf1xtvvJGonqKB3FP0EDcmpLM41+6WHdwEAjiP+W1IAJSktWvXKhqNKhQKKRaLac+ePYnXSkpKVFlZmbR/XV2dRo4cqVAopNGjR3f7RtBDhgxROBxWaWmpjhw5knh948aNKW/osKyenxIayD2EkvQW59r76g4BEHAe89ugAJiJaCD3EErSW5xr76s7BEDAecxvAqAtNJB7CCXpLc6199UdAiDgPOY3AdAWGsg9hJL0Fufa++oOARBwHvObAGgLDeQeQkl6i3PtfXWHAAg4j/lNALSFBnIPoSS9xbn2vrpDAAScx/wmANpCA7mn6KEdii56RgNH3aWBo+5SdNEzng9oPxfn2vvqTltbm6ZMmaIpU6Yk3sYKgD3MbwKgLTSQe7wexhRlYgFwBvObAGgLDeQerwctRZlYAJzB/CYA2kIDucfrQUtRJhYAZzC/CYC20EDuKXqIGxPSWZxrc+vKa8NNIIAzmN8EQFtoIPcQSrwLGZxrs4oACDiP+U0AtIUGcg+hxLuQwbk2qwiAgPOY3wRAW2gg9xBKvAsZnGuzigAIOI/5TQC0hQZyD6HEu5DBuTarCICA85jfBEBbaKC+IZSYV5xrc4sACDiP+U0AtIUG6htCiXnFuTa3CICA85jfBEBbaKC+6engiy56RgNuH6sBt4/l8WQuF+fa3Lry2vAoOMAZzG8CoC00UN94PVApKlMLgDOY3wRAW2igvvF6iFJUphYAZzC/CYC20EB94/UQpahMLQDOYH4TAG2hgfqmp8Ou8IGtui47rOuyw9yY4HJxrs2tK68NN4EAzmB+EwBtoYH6pjeDjztT0xcyONdmFncBA85jfhMAbaGB+oZQYl5xrs0tAiDgPOY3AdAWGqhvCCXmFefa3CIAAs5jfhMAbaGB+oZQYl5xrs0tAiDgPOY3AdAWGqhvCCXmFefa3CIAAs5jfhMAbaGB+oZQYl5xrs0tAiDgPOY3AdAWGqhvejz4Fj6tcGGxwoXFKlz4tOeD2M/FuTa3rrw258+f93r5Ar7A/CYA2kID9Y3XA5WiMrUAOIP5TQC0hQbqG6+HKEVlagFwBvObAGgLDdQ3Xg9RisrUAuAM5jcB0BYaqCsnh13hA1vVb0Cu+g3I5cYEl4tzbW5deW24CQRwBvObAGgLDdSV04OPO1PTFzI412YWdwEDzmN+EwBtoYG6IpRkZnGuzS0CIOA85jcB0BYaqCtCSWYW59rcIgACzmN+EwBtoYG6IpRkZnGuza3eXhsA7475TQC0hQbqilCSmcW5NrcIgIDzmN8EQFtooK4IJZlZnGtziwAIOI/5TQC0hQbqytHBt/BphQpGKFQwgseTuR0yONfGVm+vDYB3x/wmANpCA3Xl9bCkKOrdCwg65rdBAXDdunUqKipSOBxWLBbT3r17r7l/XV2dRo0apXA4rOLiYu3cuTPp9c7OTi1fvlwFBQXKyclRaWmpjh49mrTPN77xDY0fP14DBgxQXl5er4+ZBurK68FGUdS7FxB0zG9DAmBtba1CoZA2bNigQ4cOaebMmYpEIjpz5kzK/RsaGpSVlaXq6mo1NTVp2bJlys7O1oEDBxL7rFq1Snl5edq+fbv27dunsrIyDR8+XG1tbYl9VqxYoUceeUQLFy4kADrE68FGUdS7FxB0zG9DAmAsFtO8efMSH8fjcQ0dOlQrV65MuX9FRYUmTZqUtG3cuHGaPXu2pMs//SsoKNDq1asTr7e0tCgcDmvLli1dvt7GjRsJgA5xckgVLtyqrNx8ZeXmq3AhNya4WZxrc8uNawMEHfPbgADY0dGhrKwsbdu2LWn7jBkzVFZWlvJzCgsLtWbNmqRtK1as0JgxYyRJx44dk2VZamxsTNpnwoQJWrBgQZev19MA2N7ertbW1kQ1NzcHvoGu5ujg487U9IUMzrWx5ca1AYKOAGhAADx58qQsy9Lu3buTti9ZskSxWCzl52RnZ2vz5s1J29avX6/8/HxJl39FbFmWTp06lbTP1KlTVVFR0eXr9TQAVlVVJb4RX1lBbqCrmT74KM51phUBEHAeAZAAKImfADrJ9MFHca4zrQiAgPMIgAYEwEz6FfDVaKCuTB98FOc604oACDiP+W1AAJQu3wQyf/78xMfxeFzDhg275k0gkydPTto2fvz4LjeB1NTUJF5vbW3lJpA0MH3wUZzrTCsCIOA85rchAbC2tlbhcFibNm1SU1OTZs2apUgkotOnT0uSpk+frqVLlyb2b2hoUP/+/VVTU6PDhw+rqqoq5dvARCIR1dfXa//+/SovL+/yNjCvvfaaGhsb9dWvflU33HCDGhsb1djYqLfffrtHx00DdWX64KM415lWBEDAecxvQwKgJK1du1bRaFShUEixWEx79uxJvFZSUqLKysqk/evq6jRy5EiFQiGNHj262zeCHjJkiMLhsEpLS3XkyJGkfSorK1Pe1LFr164eHTMN1JWjg2/hVmXfElX2LVHemsTtkMG5NrbcuDZA0DG/DQqAmYgG6srrYUlR1LsXEHTMbwKgLTRQV14PNoqi3r2AoGN+EwBtoYG68nqwURT17gUEHfObAGhL0Boo3UOKv0vjXFP8DSDghqDN71QIgDYErYHSPvi4M5VzTXEXMOCCoM3vVAiANgStgfww+CjOdaYVARBwXtDmdyoEQBuC1kB+GHwU5zrTigAIOC9o8zsVAqANQWsgPww+inOdaUUABJwXtPmdCgHQhqA1kB8GH8W5zrTy6toAfha0+Z0KAdCGoDVQUAZfEItzbW4RAAHnBW1+p0IAtMFPDeT1kEs5+BZuVVZuvrJy83lrEs51YMurawP4mZ/md18RAG3wUwN5PeQoijKrAD/z0/zuKwKgDX5qIK+HDUVRZhXgZ36a331FALTBTw3k9bChKMqsAvzMT/O7rwiANvipgbweNqmqcOHTChWMUKhghAoXPu358fi5ONfmllfXBvAzP83vviIA2uCnBvJ6yKUcfNyZyrmmuAsYcIGf5ndfEQBt8FMDeT3kTBp8QSzOtblFAASc56f53VcEQBv81EBeDzmTBl8Qi3NtbhEAAef5aX73FQHQBj81kNdDzqTBF8TiXJtbJl8bIFP5aX73FQHQBj81kNeDJNMGn9+Kc21umXxtnPreAqSbn+Z3XxEAbfBTA3k9SDJt8PmtONfmlsnXxqnvLUC6+Wl+9xUB0AY/NZDXg6S7wddvQK76Dcg1bvD5rTjX5pbJ18ap7y1AuvlpfvcVAdAGPzWQ14OEoqjMK6e+twDp5qf53VcEQBv81EBeDxKKojKvnPreAqSbn+Z3XxEAbfBTA3k9SCiKCm4B6ean+d1XBEAb/NRAXg+AVFW48GmFC4sVLizm8WSc68BWEK4NkG5+mt99RQC0wU8N5PUASDn4DL770W/FuTa3gnBtgHTz0/zuKwKgDX5qIK8HQFAHnynFuTa3gnBtgHTz0/zuKwKgDX5qIK8HQFAHnynFuTa3gnBtgHTz0/zuKwKgDX5qIK8HQFAHnynFuTa3gnBtgHTz0/zuKwKgDX5qIK8HQFAHnynFuTa3gnBtgHTz0/zuKwKgDZnSQF5/c2fwmV+ca3MrCNcGSLdMmd9uIgDakCkN5PU3dzuD77rssK7LDvt28JlSnGtzKwjXBki3TJnfbiIA2pApDeT1N3eKoii7BTgpU+a3mwiANmRKA3n9jZuiKMpuAU7KlPntJgKgDZnSQF5/46YoirJbgJMyZX67iQBoQ6Y0kNffuPta0UXPaMDtYzXg9rGKLnrG8+Pxc3GuzS2uTc/Lqe+H8L9Mmd9uIgDakCkN5PU35b5WEO5+NKU41+YW1yb9Bf/LlPntJgKgDSY0kNffKN0sBh/nmuLaeFFOfe+FuUyY314jANpgQgN5/Y3SzWLwca4pro0X5dT3XpjLhPntNWMC4Lp161RUVKRwOKxYLKa9e/dec/+6ujqNGjVK4XBYxcXF2rlzZ9LrnZ2dWr58uQoKCpSTk6PS0lIdPXo0aZ+33npLn/nMZ3TjjTcqLy9P9913n95+++0eH7MJDeT1N0o3i8HHuaa4Nl5UOr/3whsmzG+vGREAa2trFQqFtGHDBh06dEgzZ85UJBLRmTNnUu7f0NCgrKwsVVdXq6mpScuWLVN2drYOHDiQ2GfVqlXKy8vT9u3btW/fPpWVlWn48OFqa2tL7HP33Xfrzjvv1J49e/Tiiy/qfe97n6ZNm9bj43a7gbz+Juh1Mfg41xTXxu9lmkw85r4gABoSAGOxmObNm5f4OB6Pa+jQoVq5cmXK/SsqKjRp0qSkbePGjdPs2bMlXf7pX0FBgVavXp14vaWlReFwWFu2bJEkNTU1ybIs/fznP0/s8/zzz+u6667TyZMne3TcBEB3i8HHuaa4Nn6vdM6LdH4d0xEADQiAHR0dysrK0rZt25K2z5gxQ2VlZSk/p7CwUGvWrEnatmLFCo0ZM0aSdOzYMVmWpcbGxqR9JkyYoAULFkiSnnjiCUUikaTXL168qKysLD3zzDMp/9329na1trYm6re//a0sy1Jzc3PSdqeq8At1ga5hc/9fYvANm/v/PD8ePxfn2tzi2lCZWE7NODdma2trq5qbm2VZllpaWnoeWHzG8wB48uRJWZal3bt3J21fsmSJYrFYys/Jzs7W5s2bk7atX79e+fn5ki7/itiyLJ06dSppn6lTp6qiokKS9M1vflMjR47s8rUHDx6sRx99NOW/W1VVlfhGTFEURVFUZldzc3PPwooPEQCvcq0AePVPAP/4xz/q2LFjamlpce3/Uqg//5+aWz9ppbg+fi6uj9nF9fGmWlpa1NzcrHg83rOw4kOeB8BM+hUwvNHayt9qmIzrYzauj9m4PvCK5wFQunwTyPz58xMfx+NxDRs27Jo3gUyePDlp2/jx47vcBFJTU5N4vbW1NeVNIK+88kpinxdeeKFXN4EgPfgGaTauj9m4Pmbj+sArRgTA2tpahcNhbdq0SU1NTZo1a5YikYhOnz4tSZo+fbqWLl2a2L+hoUH9+/dXTU2NDh8+rKqqqpRvAxOJRFRfX6/9+/ervLw85dvAfOhDH9LevXv1s5/9TCNGjOjV28AgPfgGaTauj9m4Pmbj+sArRgRASVq7dq2i0ahCoZBisZj27NmTeK2kpESVlZVJ+9fV1WnkyJEKhUIaPXp0t28EPWTIEIXDYZWWlurIkSNJ+7z11luaNm2abrjhBuXm5uree+/t1RtBIz3a29tVVVWl9vZ2rw8FKXB9zMb1MRvXB14xJgACAAAgPQiAAAAAAUMABAAACBgCIAAAQMAQAAEAAAKGAAhP/PSnP9XkyZP1nve8R5ZldXkj8Hfu4i4oKFBOTo5KS0t19OjRpH3eeustfeYzn9GNN96ovLw83XfffdzF7ZCHH35YY8eO1Q033KDBgwervLxcr776atI+bW1tmjt3rm6++WZdf/31+tSnPpV466Z3vPbaa5o4caIGDBigwYMHa/Hixbp48WI6/1N86dFHH9Vf/MVf6MYbb9SNN96oj3zkI3ruuecSr3NtzLJy5UpZlqX7778/sY1rBK8RAOGJ5557Tl/+8pf1zDPPpAyAq1atUl5enrZv3659+/aprKws5fs43nnnndqzZ49efPFFve997+N9HB3yd3/3d9q4caMOHjyoX/7yl5o4caKi0ajOnj2b2GfOnDkqLCzUj3/8Y73yyiv6yEc+oo9+9KOJ1y9duqTi4mJ9/OMfV2Njo5577jkNGjRIX/ziF734T/KVZ599Vjt37tTRo0d15MgRfelLX1J2drYOHjwoiWtjkpdfflm33XabxowZkxQAuUbwGgEQnrs6AL7zJJfVq1cntrW0tKR8ksvPf/7zxD7PP/88T3JxyZtvvinLsvTTn/5U0uXrkZ2dre9///uJfQ4fPizLsvTSSy9Juhzy+/Xrl/RTjccee0y5ubnq6OhI739AANx00016/PHHuTYGefvttzVixAj9z//8j0pKShIBkGsEExAA4bmrAyDPcjbPr371K1mWlXjazo9//GNZlqU//vGPSftFo1E98sgjkqTly5frzjvvTHr9N7/5jSzL0i9+8Yv0HHgAXLp0SVu2bFEoFNKhQ4e4NgaZMWOGvvCFL0hSUgDkGsEEBEB47uoA2NDQIMuydOrUqaT9pk6dqoqKCknSN7/5TY0cObLL1xo8eLAeffRRdw84YOLxuCZNmqS77rorse2pp55SKBTqsu+HP/xhPfjgg5KkmTNn6pOf/GTS6+fOnZNlWUl/r4a+2b9/v66//nplZWUpLy8v8TQkro0ZtmzZouLi4sSfrVwZALlGMAEBEJ4jAJptzpw5KioqUnNzc2IbA8x7HR0d+tWvfqVXXnlFS5cu1aBBg3To0CGujQF++9vfKj8/X/v27UtsIwDCNARAeI5fAZtr3rx5uvXWW/Wb3/wmaTu/wjJPaWmpZs2axbUxwLZt22RZlrKyshJlWZauu+46ZWVl6Uc/+hHXCJ4jAMJz3d0EUlNTk9jW2tqa8iaQV155JbHPCy+8wE0gDuns7NS8efM0dOjQLm+/I/35BhpTUgAABtJJREFUj9i3bt2a2Pbqq6+m/CP2M2fOJPb5zne+o9zcXB5874KPfexjqqys5NoY4E9/+pMOHDiQVGPHjtVnP/tZHThwgGsEIxAA4Ym3335bjY2NamxslGVZeuSRR9TY2KjXXntN0uW3gYlEIqqvr9f+/ftVXl6e8m1gPvShD2nv3r362c9+phEjRvA2MA75l3/5F+Xl5el///d/9cYbbyTq/PnziX3mzJmjaDSqn/zkJ3rllVc0fvx4jR8/PvH6O29j8clPflK//OUv9d///d8aPHgwb2PhgKVLl+qnP/2pjh8/rv3792vp0qW67rrr9MMf/lAS18ZEV/4KWOIawXsEQHhi165dsiyrS1VWVkr68xtBDxkyROFwWKWlpTpy5EjS13jrrbc0bdo03XDDDcrNzdW9997LG0E7JNW1sSxLGzduTOzzzhvZ3nTTTRo4cKD+4R/+QW+88UbS1zlx4oTuueceDRgwQIMGDdKiRYt4I1sH3HfffSoqKlIoFNLgwYNVWlqaCH8S18ZEVwdArhG8RgAEAAAIGAIgAABAwBAAAQAAAoYACAAAEDAEQAAAgIAhAAIAAAQMARAAACBgCIAAAAABQwAEAAAIGAIgAABAwBAAAWS8r3/96xo3blyvP+/3v/+9Bg8erOPHjzt/UD3w6U9/WjU1NZ782wCCjQAIIONNmTJFs2fP7vXnPfDAA/r85z+f+LiyslKWZaX8WnPnzk16XrUTDhw4oJtuukktLS2OfU0A6AkCIICMN2LECD322GO9+pxz584pNzdXL730UmJbZWWlCgsLlZeXp/Pnzye2t7W1KRKJKBqNOhoAJWns2LFat26do18TAN4NARBARnnllVf013/918rJydEHP/hB7dmzR/369UsKcj3x/e9/X4MHD07aVllZqfLychUXF+vJJ59MbH/qqac0ZswYlZeXJwXAkpISzZs3T/PmzVNubq5uueUWLVu2TJ2dnYl94vG4/u3f/k3vfe97FQqFVFhYqG984xuJ17/61a/qr/7qr3p5FgDAHgIggIxx+PBh3XjjjVq2bJl+/etfa+vWrSooKFC/fv107ty5Xn2tBQsW6O67707a9k4AfOSRR1RaWprYXlpaqjVr1qQMgDfccIPuv/9+vfrqq3ryySc1cOBAffe7303s8+CDD+qmm27Spk2b9Otf/1ovvviivve97yVef/755xUKhdTe3t7LswEAfUcABJAx/vZv/1bTp09P2jZlyhSNGjVKkvStb31LX/ziF3v0tcrLy3XfffclbXsnAL755psKh8M6ceKETpw4oZycHP3ud79LGQDvuOOOpJ/4PfTQQ7rjjjskSX/6058UDoeTAt/V9u3bJ8uydOLEiR4dNwA4gQAIICOcOHFClmXpF7/4RdL2adOm6dOf/rQk6d5779VTTz3Vo6/3yU9+UnPnzk3a9k4AlKRPfepT+spXvqKqqir94z/+oySlDID33ntv0tfYvn27+vfvr0uXLmnv3r2yLEu/+c1vuj2Oo0ePyrIsNTU19ei4AcAJBEAAGeGdYHXx4sWk7e9///u1cuVKSdKHP/xh7du3r0df7zOf+YymTZuWtO3KALhjxw7ddtttuu2227Rz505JvQ+A+/fvf9cAuGfPHlmWpd/97nc9Om4AcAIBEEBGePbZZ7v8rd/OnTtlWZaef/55dXZ2KhKJqKOjo0dfb/Xq1brzzjuTtl0ZAC9duqShQ4dq2LBhunTpkqTUAfADH/hA0tdYunRp4lfAbW1tGjBgwDV/Bfz444/r1ltv7dExA4BTCIAAMsLJkycVDoc1d+5cHTt2TPX19Ro6dKgsy9Ibb7yhY8eOJYJXT+zfv1/9+/fXH/7wh8S2KwOgJLW2tqq1tTXxcXc3gTzwwAN69dVXtXnzZl1//fX69re/ndjnK1/5im666Sb9x3/8h37961/rpZde0uOPP570b179t4gA4DYCIICM8Z//+Z+69dZblZOTowkTJmjp0qUaMmSIJKm+vl5Tpkzp1deLxWJJYe3qAHi1VAFw7ty5mjNnjnJzc3XTTTfpS1/6Upe3gfnGN76hoqIiZWdnKxqN6uGHH5Z0+SeEeXl5vX4LGwCwiwAIwBe++c1v6itf+UqvPmfHjh264447FI/H+/RvlpSU6P777+/T50rSo48+qk984hN9/nwA6CsCIABfmDZtmm655RYVFRWpqKioxz8NXLNmjX7729/26d+0GwC/973v6dVXX+3z5wNAXxEAAaCP7AZAAPAKARAAACBgCIAAAAABQwAEAAAIGAIgAABAwBAAAQAAAoYACAAAEDAEQAAAgIAhAAIAAATM/wdzvLKh4jF5jQAAAABJRU5ErkJggg==\" width=\"640\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "distance = fp['samples/distance'][:]\n",
    "percentiles = numpy.percentile(distance, [5, 50, 95])\n",
    "fig, ax = pyplot.subplots()\n",
    "ax.hist(distance, density=True, bins=50)\n",
    "for dl in percentiles:\n",
    "    ax.axvline(dl, linestyle='--', color='k')\n",
    "ax.set_xlabel('$d_L$ (Mpc)')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metadata about the run, such as the parameters that were kept fixed through the run, are stored in the file's `.attrs`. In addition, the log likelihood of the noise is stored in the `samples` group's `.attrs` as `lognl`. Subtracting this from the `loglikelihood` gives the log likelihood ratio, from which we can get the coherent matched filter SNR. For example, the maximum likelihood SNR is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.955152900190067\n"
     ]
    }
   ],
   "source": [
    "loglikelihood = fp['samples/loglikelihood'][:]\n",
    "loglr = loglikelihood - fp['samples'].attrs['lognl']\n",
    "# the coherent SNR = sqrt(2 * loglr), Note: This is marginalized over phase, which is less than the peak SNR before marginalization\n",
    "print((2*loglr.max())**0.5)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}