{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Modeling of the thickness of the sensor\n",
"\n",
"In this notebook we will re-use the experiment done at ID28 and previously calibrated and model in 3D the detector.\n",
"\n",
"This detector is a Pilatus 1M with a 450µm thick silicon sensor. Let's first have a look at the absorption coefficients of this sensor material: https://physics.nist.gov/PhysRefData/XrayMassCoef/ElemTab/z14.html\n",
"\n",
"First we retieve the results of the previous step, then calculate the absorption efficiency:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"%pylab nbagg"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/users/kieffer/VirtualEnvs/py3/lib/python3.5/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
" from ._conv import register_converters as _register_converters\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"wavelength: 6.968e-11m,\t dist: 2.845e-01m,\t poni1: 8.865e-02m,\t poni2: 8.931e-02m,\t energy: 17.793keV\n"
]
}
],
"source": [
"import numpy\n",
"import fabio, pyFAI, pyFAI.units, pyFAI.detectors, pyFAI.azimuthalIntegrator\n",
"import json\n",
"with open(\"id28.json\") as f:\n",
" calib = json.load(f)\n",
"\n",
"thickness = 450e-6\n",
"wavelength = calib[\"wavelength\"]\n",
"dist = calib[\"param\"][calib['param_names'].index(\"dist\")]\n",
"poni1 = calib[\"param\"][calib['param_names'].index(\"poni1\")]\n",
"poni2 = calib[\"param\"][calib['param_names'].index(\"poni2\")]\n",
"energy = pyFAI.units.hc/(wavelength*1e10)\n",
"print(\"wavelength: %.3em,\\t dist: %.3em,\\t poni1: %.3em,\\t poni2: %.3em,\\t energy: %.3fkeV\" % \n",
" (wavelength, dist, poni1, poni2, energy))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Absorption coeficient at 17.8 keV"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"µ = 1537.024214 m^-1 hence absorption efficiency for 450µm: 49.9 %\n"
]
}
],
"source": [
"# density from https://en.wikipedia.org/wiki/Silicon\n",
"rho = 2.3290 # g/cm^3\n",
"\n",
"#Absorption from https://physics.nist.gov/PhysRefData/XrayMassCoef/ElemTab/z14.html\n",
"# Nota: enegies are in MeV !\n",
"Si_abs = \"\"\"\n",
" 2.00000E-03 2.777E+03 2.669E+03 \n",
" 3.00000E-03 9.784E+02 9.516E+02 \n",
" 4.00000E-03 4.529E+02 4.427E+02 \n",
" 5.00000E-03 2.450E+02 2.400E+02 \n",
" 6.00000E-03 1.470E+02 1.439E+02 \n",
" 8.00000E-03 6.468E+01 6.313E+01 \n",
" 1.00000E-02 3.389E+01 3.289E+01 \n",
" 1.50000E-02 1.034E+01 9.794E+00 \n",
" 2.00000E-02 4.464E+00 4.076E+00 \n",
" 3.00000E-02 1.436E+00 1.164E+00 \n",
" 4.00000E-02 7.012E-01 4.782E-01 \n",
" 5.00000E-02 4.385E-01 2.430E-01 \n",
" 6.00000E-02 3.207E-01 1.434E-01 \n",
" 8.00000E-02 2.228E-01 6.896E-02 \n",
" 1.00000E-01 1.835E-01 4.513E-02 \n",
" 1.50000E-01 1.448E-01 3.086E-02 \n",
" 2.00000E-01 1.275E-01 2.905E-02 \n",
" 3.00000E-01 1.082E-01 2.932E-02 \n",
" 4.00000E-01 9.614E-02 2.968E-02 \n",
" 5.00000E-01 8.748E-02 2.971E-02 \n",
" 6.00000E-01 8.077E-02 2.951E-02 \n",
" 8.00000E-01 7.082E-02 2.875E-02 \n",
" 1.00000E+00 6.361E-02 2.778E-02 \n",
" 1.25000E+00 5.688E-02 2.652E-02 \n",
" 1.50000E+00 5.183E-02 2.535E-02 \n",
" 2.00000E+00 4.480E-02 2.345E-02 \n",
" 3.00000E+00 3.678E-02 2.101E-02 \n",
" 4.00000E+00 3.240E-02 1.963E-02 \n",
" 5.00000E+00 2.967E-02 1.878E-02 \n",
" 6.00000E+00 2.788E-02 1.827E-02 \n",
" 8.00000E+00 2.574E-02 1.773E-02 \n",
" 1.00000E+01 2.462E-02 1.753E-02 \n",
" 1.50000E+01 2.352E-02 1.746E-02 \n",
" 2.00000E+01 2.338E-02 1.757E-02 \"\"\"\n",
"data = numpy.array([[float(i) for i in line.split()] for line in Si_abs.split(\"\\n\") if line])\n",
"energy_tab, mu_over_rho, mu_en_over_rho = data.T\n",
"abs_18 = numpy.interp(energy, energy_tab*1e3, mu_en_over_rho) \n",
"mu = abs_18*rho*1e+2\n",
"eff = 1.0-numpy.exp(-mu*thickness)\n",
"\n",
"print(\"µ = %f m^-1 hence absorption efficiency for 450µm: %.1f %%\"%(mu, eff*100))\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"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,iVBORw0KGgoAAAANSUhEUgAAAroAAAILCAYAAAAHaz/JAAAgAElEQVR4nOzdd3RUdf7/8ZuEGqq4wBKFQRFEsa3AWkM0wGJXUHTX8nXXFRv+rOCAoKFo1KWoiC5IVRcF2QXLAKEZegsg0gLSIbRQkgAhbWZevz+QKyHcNHJzZybPxzn3nOVm7twP8P7+fk/GO/caAgAAAEKQ4fQCAAAAADsQugAAAAhJhC4AAABCEqELAACAkEToAgAAICQRugAAAAhJhC4AAABCEqELAACAkEToAgAAICQRugAAAAhJhC4AAABCEqELAACAkEToAgAAICQRugAAAAhJhC4AAABCEqELAACAkEToAgAAICQRugAAAAhJhC4AAABCEqELAACAkEToAgAAICQRugAAAAhJhC4AAABCEqELAACAkEToAgAAICQRugAAAAhJhC4AAABCEqELAACAkEToAgAAICQRugAAAAhJhC4AAABCEqELAACAkEToAgAAICQRugAAAAhJhC4AAABCEqELAACAkEToAgAAICQRugAAAAhJhC6AchcTEyOXy5Vv35NPPinDyP//JMXFxckwDO3YsaP8Focyxd8hACcRugDKxP79+/X666+rVatWqlmzpmrVqqXLLrtMjzzyiP73v//le21FC12fz6fjx4+f13scP35c/fr10z333KNGjRrJMAx16tTpnK/dsWOHDMModPvPf/5TrPOuXbtWDz/8sJo0aaJq1arJ5XLpscce07p164p1fFn/HW7YsEGGYeiuu+4q9HU//PCDDMPQq6++WibnBRCcCF0A523nzp1q0KCBqlWrpn/+858aPny4hg8frldeeUXXXnut7rnnnnyvz8nJUXZ2dr595wrdvLw8ZWVlye/32/57KGsbN27UCy+8oGbNmik8PFyGYSgyMlI33XSThgwZUuLwPR2vjRo10r333lto6J44cUJfffXVObc6deqoSpUqOnToUJHnXLVqlapVq6aoqCj1799fo0ePltvtVp06dVS9enWtX7++yPew4x8rN954oyIiIrR3717L19x///0yDEMbNmwos/MCCD6ELoDz9uKLL8owDH333Xfn/Pn+/fuLfI9zhW4wys3N1csvv6zw8HC1atVKffv21eTJkzV9+nSNHz9ezzzzjP7whz/ooosuUkJCQrHfNzs7WykpKeavCwtdK0uWLJFhGHrooYeK9fqnnnpKhmFozZo1+fZPnTpVhmHI7XYX+R52hO6oUaNkGIbi4+PP+fODBw+qUqVKuvHGG8vsnACCU/D//yoAHNepUycZhlHsTynP99KFjIwMvfnmm2rZsqWqVq2qevXq6ZZbbtE333yT73Xr16/Xgw8+qAsvvFBVqlRRixYt1L9//wKfJp8+z6ZNm9S7d29ddNFFqlKliq655hpNmzateH8Ikrxer+68807Vrl27wFrOXn/37t1VpUqVEr3/mUoTuqfDtbiB/dBDD8kwDKWmpubbv3z5chmGoQEDBhT5Huf6O/T5fHrppZdkGIbi4uLyvX727Nnq2LGj6tSpo6pVq+rqq6/Wv//973yvOX78uGrWrKnmzZuf85yDBg2SYRgaNWpUsX6fAEIXoQvgvL3wwgsyDENDhw4t1mUG5xO6aWlpatWqlQzD0IMPPqiPP/5YQ4cO1RNPPKHHHnvMfN2qVavMa4V79eql4cOHm//J/4477pDP5ytwnhtuuEG33nqrPvzwQ33wwQdq1KiRKleuXOxPI/v06aPatWsX+M/lJ06cUF5enqRTn8wePXpUkvT222/rggsuKBCSxVHS0D0dhy6XK9/vvTCnPznt1KmTli9frpSUFCUmJur6669X48aNi/VJ/dl/h1lZWXrwwQcVERFRIERHjhypsLAw3XDDDfrXv/6lzz77TJ07d5ZhGOrRo0e+156O9gULFhQ455VXXqkaNWro2LFjxfp9AghdhC6A87Zt2zbVrl1bhmGocePGevTRR/Xhhx9q1apV53z9+YTu888/L8Mw9OmnnxZ43zMD7pZbblF4eHiBNXTr1k2GYWjChAkFznP33XfnC/UVK1bIMAz16tWryD+DvXv3qlq1avriiy/MfUuWLNG1114rwzBUuXJlPfroo/roo4/M37vP51OrVq00cODAIt//bCUN3dGjR8swDPXr16/Yx/h8PvXt21c1a9bM90W2W2+9VQcOHCjWe5z5d3j06FHdeuutqlGjhjweT77X7du3T1WrVlXXrl0LvMdLL72k8PBwbdu2zdy3ePFiGYahv//97/leu3TpUhmGoX/84x/F/n0CCF2ELoAysWPHDnXv3l2NGzfOF0VXX321Vq5cme+1pQ1dn8+nCy64QC1atCj0k+ODBw/KMAzde++9BX62e/du89Pgs88za9asAq+vWbNmvtdaGTZsmC6++GJ5vV5J0qFDh1SvXj1de+21+uKLLzR58mR16dJFNWrUyPd779evn2699dYi3/9sJQ3dG2+8UeHh4dq1a1eJzjN8+HB16tRJn376qb777jv169dPtWrVUtu2bZWenl7k8af/bOfNm6crrrhC9evX14oVKwq8btiwYeZlFYcOHcq3zZ49W4ZhaOTIkfmOueKKK1SjRo18l8yc/ofMokWLSvT7BBCaCF0AZW7fvn369ttvzUsF/vjHP+rIkSPmz0sbuqcD9lyf+p1p2bJlhX4SW6dOHf3pT38qcJ6tW7cWeK3L5dJtt91W6Pkk6bHHHtMTTzxh/nrIkCGqVauWDh8+nO91N910U77f+4gRIyyvNS1MSUL39C25SnpN7+k7LJx9icKMGTNkGIZ69+5d5Huc/rOtXbu2qlWrpk2bNp3zdac/qS9sO/ua4MGDB8swDI0ePVqSlJmZqdq1a+vyyy8v0e8TQOgidAHY6tFHH5VhGPrqq6/MfYEauue6FtflcikmJqbQ80lSx44d853vhRdeUNu2bQu8rkePHvl+771799bNN99c5PufrSTh+tprr8kwDE2ePLnY75+bm6sqVaoUuDXcabVq1SrWXQ1O/9k+++yzMgxDjz32mPmp95mee+45GYahcePGafbs2efczrx0QZJSU1NVuXJl88/viy++kGEY+te//lXs3yeA0EboArDV6f8k/f7775v77L50ITU11fLShT179hS4xVZZhG7Xrl313HPPmb9+5513VL9+feXk5OR73b333mv+3k+cOKGLLrrI1mt0c3JyVL9+fdWvX1+5ubnFfv99+/bJMAzdeeedBX7m9/sVGRmpNm3aFPk+Z/7Znv4EtmvXruaX804bMmSIDMPQDz/8UOw1StKDDz5o3jEjJiZGlSpVKvb1wwBCH6EL4LwlJibq5MmTBfb7fD517NixQMCcz5fRTt/hYcSIEQXOd2b83nLLLYqIiChwD9jTnxx+/fXXhZ7ntOKGbv/+/XXNNdeYv960aZOqVKmirl276pdfftHWrVs1cOBAhYeH6+KLL9bs2bPVpk0bXX755crIyCjy/c9W3NCdPHmyDMPQa6+9ZvmazMxMJScna9++feY+n8+nCy+8UDVq1ND27dvzvX7ixIkyDEMvvPBCkec/+8/2k08+UVhYmB544IF8/wjYs2ePqlatquuvv16ZmZkF3ic9Pb3AbeEkafr06eY/XE6/LwCcRugCOG9333236tatqyeeeEIffvihxo4dq/j4eLVu3VqGYej222/Pd0eE8wndo0ePqmXLluYng8OGDdPHH3+sf/zjH3r88cfN161atUo1atRQ7dq11bt3bw0fPlz33XefGYjnur3Y+YTuunXrFBYWpvnz55v7xo4dqxo1apjXmF5xxRXq06ePeReGxx57rES3Fvvkk080cOBADRw4UIZhqHnz5uavv/zyy3Mec8cdd8gwDG3cuNHyfRMTE2UYhp588skC5zMMQ/Xr19dbb72lkSNH6vnnn1flypV14YUXaufOnUWu+Vx/tp9//rnCwsJ011135YvXsWPHKjw8XI0bN1bfvn01atQoxcfH629/+5uqV69+zr8fn8+X7wuQP/74Y5FrAlBxELoAztvSpUv12muvqU2bNmrQoIEqVaqkOnXq6MYbb9SQIUMKfBJ3vg+MSEtLU8+ePdWsWTNVrlxZ9erV06233qpJkyble926devUpUsX1atXT5UrV1bz5s3Vr18/ywdGnE/oSqcuX2jevHm+L6BlZGRo4cKFWrNmjXw+n/bv36+kpKRzfmpZFJfLZflFrXOtcffu3QoPDy/yGmCr0JWkWbNm6S9/+YsuuugiVa5cWY0aNdITTzxR4FNeK1Z/tuPHj1d4eLg6duyY778GLFq0SA888IDq169vnu+2227T4MGDlZWVdc5zvPXWWzIMQ1FRUee8/hdAxUXoAkAZOXTokJo1a6arrrrK8u4C0qlrc0vyxTAAQOkQugBQhvbs2aM///nPqly5sv7+97/ru+++06ZNm7Rr1y4tWrRIcXFxatSokRo3blyqJ6IBAIqP0AWAMpaXl6fPP//cfCramVvjxo3Vv3//Yj1sAQBwfghdALBRamqqkpKStHDhwgL3gQUA2IvQBQAAQEgidAEAABCSCF0AAACEJEIXAAAAIYnQtUlOTo6SkpK0a9cupaSksLGxsbGxsbHZsu3atUtJSUn5HquNUwhdmyQlJVk+wYiNjY2NjY2Nray3pKQkp/Mn4BC6Ntm1a5c5dE7/S4+NjY2NjY0tdLfTH67t2rXL6fwJOISuTVJSUmQYhlJSUpxeCgAACGE0hzVC1yYMHQAAKA80hzVC1yYMHQAAKA80hzVC1yYMHQAAKA80hzVC1yYMHQAAKA80hzVC1yYMHQAAKA80hzVC1yYMHQAAKA80hzVC1yYMHQAAKA80hzVC1yYMHQAAKA80hzVC1yYMHQAAKA80hzVC1yYMHQAAKA80h7WQDd0tW7bomWee0dVXX63w8HDFxMQU67icnBz16NFDDRs2VGRkpDp06KBNmzaV+PwMHQAAKA80h7WQDd3vvvtOjRs3VteuXdWiRYtih+6zzz6rOnXqaMyYMUpISFB0dLQuuugipaenl+j8DB0AACgPNIe1kA1dn89n/u/777+/WKG7Z88eRURE6PPPPzf3HTlyRDVq1NAHH3xQovMzdAAAoDzQHNZCNnTPVNzQHTNmjMLCwpSWlpZvf+fOnYv9ifBp5TV0fr/f1vcHAACBjdC1RuieoWfPnrrooosK7H/zzTfVsGHDEp3T7qE7np0n939/Ufy0jba8PwAACA6ErjVC9wxPP/20WrVqVWD/oEGDVLly5UKPzcjIUEpKirklJSXZOnTDf9oil9ujpr08WrbtsC3nAAAAgY/QtUbonuF8QjcuLk6GYRTY7Bq67DyvOn04Xy63R7e8P1fHsnJtOQ8AAAhshK41QvcM53PpQnl/oitJG/Zm6LI3p8nl9qjn5DW2nQcAAAQuQtcaoXuG019GO/tWYl26dAnYL6N9lrhVLrdHLrdHszYcsPVcAAAg8BC61gjdM5y+vdjo0aPNfUePHlXNmjUD9vZiXp9fD362WC63R60HztLh49m2ng8AAAQWQtdayIZuZmamJk+erMmTJ6tt27a68sorzV+npqZKkmJjYxUbG5vvuGeffVZ169bV2LFjNXPmTMXExAT8AyN2Hj6hK96aIZfbo25fJHHLMQAAKhBC11rIhu6OHTvO+eUwwzCUmJgoSYqJiSnwSW92drZef/11NWjQQNWrV1eHDh2UnJxc4vOX99B9vXyXeQnDt0m7y+WcAADAeYSutZANXaeV99D5/X49NW6FXG6PWr2doN1HMsvlvAAAwFmErjVC1yZODN3BY1n604BZcrk9eujfi+X1cQkDAAChjtC1RujaxKmhm7l+v3kJw/CftpTruQEAQPkjdK0RujZxcujc//1FLrdHzXpP09o9JfsSHQAACC6ErjVC1yZODt2J7Dy1+9dPcrk9ih2cqJM53nJfAwAAKB+ErjVC1yZOD92qXUd1ae9TT01767t1jqwBAADYz+nmCGSErk0CYeiGztpsXq/706aDjq0DAADYJxCaI1ARujYJhKHL9fp0//BFvz01bTZPTQMAIAQFQnMEKkLXJoEydNsP/f7UtH+OX8FT0wAACDGB0hyBiNC1SSAN3aQVu81LGL5autPp5QAAgDIUSM0RaAhdmwTS0Pn9fj331Uq53B5d3ne6thw85vSSAABAGQmk5gg0hK5NAm3o0jJzdMO7c+Rye3TnRwuUncctxwAACAWB1hyBhNC1SSAO3eKth9S016lLGN6dttHp5QAAgDIQiM0RKAhdmwTq0MVP32her7vw10NOLwcAAJynQG2OQEDo2iRQhy4nz6e7hy2Qy+1R23e45RgAAMEuUJsjEBC6Ngnkoduaelwt+3LLMQAAQkEgN4fTCF2bBPrQTVyxy7yEYfziHU4vBwAAlFKgN4eTCF2bBPrQ+f1+vTBhlVxuj5r3ma6N+zKcXhIAACiFQG8OJxG6NgmGoUs/maub35srl9uj9kPm6WQOtxwDACDYBENzOIXQtUmwDF3SjiO65LdbjvWestbp5QAAgBIKluZwAqFrk2Aauo9m/2perztj3T6nlwMAAEogmJqjvBG6NgmmofP6/Oo6Yolcbo+ujkvQnqOZTi8JAAAUUzA1R3kjdG0SbEO3N+2kruk3Uy63R10+W6w8r8/pJQEAgGIItuYoT4SuTYJx6Gau329ewjAoYZPTywEAAMUQjM1RXghdmwTr0MV9v14ut0dNe3m0aAuPCAYAINAFa3OUB0LXJsE6dFm5Xt350alHBLd5Z7ZSj/GIYAAAAlmwNkd5IHRtEsxDtzX1uK5469Qjgp8Ys1w+H48IBgAgUAVzc9iN0LVJsA/df1fuMa/X/fe8rU4vBwAAWAj25rAToWuTUBi6Vyf+LJfbo0t7T9PKnUedXg4AADiHUGgOuxC6NgmFoTuRnafbByXK5fbo5vfmKi0zx+klAQCAs4RCc9iF0LVJqAzdhr0Zat5nulxuj/45Pkl+P9frAgAQSEKlOexA6NoklIZuwrJd5vW6oxZsc3o5AADgDKHUHGWN0LVJKA2d3+/Xi1+vlsvt0WVvTtPPu9OcXhIAAPhNKDVHWSN0bRJqQ3csK1cx//pJLrdHt7w/V+knc51eEgAAUOg1R1kK6dBNTk5Whw4dFBkZqYYNG6pnz57KySn6C1UpKSl6+OGHVbt2bdWsWVP33nuvtm/fXqJzh+LQrUtJN6/XfeZLrtcFACAQhGJzlJWQDd2jR4+qUaNGateunRISEjRmzBjVqVNH3bt3L/Q4r9erq6++WpdddpkmTZqkqVOn6rrrrlPTpk11/PjxYp8/VIfuy6U7zet1Ry8sWfwDAICyF6rNURZCNnTj4+NVs2ZNHTlyxNw3cuRIRUREaO/evZbHffPNNwoLC9P69evNfSkpKapataqGDh1a7POH6tD5/X51n7BKLrdHzXpP06pd3F8XAAAnhWpzlIWQDd3o6Gh17tw53760tDSFhYVp3Lhxlse98cYbuvjiiwvsb926tW6//fZinz+Uh+74GffXvSl+jo6e4P66AAA4JZSb43yFbOjWr19fffr0KbA/KipKbrfb8riXXnpJzZo1K7D/5ptv1h//+Mdinz/Uh27jvgy1+O163b+PXS6fj+t1AQBwQqg3x/kI2dCtVKmSBg0aVGB/q1at1K1bN8vjhg8frkqVKunAgQPmvuPHj6tOnTqqUqWK5XEZGRlKSUkxt6SkpJAfukkrdpvX636auMXp5QAAUCERutYI3bMcOXJEF1xwge644w5t375de/bs0SOPPKKIiAhVrVrV8ri4uDgZhlFgC+Wh8/v9em3SGrncHl3Sy6Ol2w47vSQAACocQtdayIZuaS9dkKRZs2bpoosuMmP1tttu01NPPaWmTZtaHlMRP9GVpMycPHUcOk8ut0dt3pmtg8eynF4SAAAVCqFrLWRDNzo6Wl26dMm3Lz09vcgvo53m9Xq1ceNG7dixQ5J011136a9//Wuxz1+Rhm7LweO64q0Zcrk9enjEEuV5fU4vCQCACqMiNUdJhWzoxsfHq1atWkpL+/1xtaNGjSry9mLnkpycrCpVqmju3LnFPqaiDd0Pa/aa1+vGT9/o9HIAAKgwKlpzlETIhu7pB0bExMRo5syZGjt2rOrWrVvggRGxsbGKjY3Nt+/111/XlClTNHfuXA0dOlT16tXT008/XaLzV8Shi/t+vRm7M9fvd3o5AABUCBWxOYorZENXkjZu3Kj27durevXqatCggXr06FHgEcAxMTGKiYnJt+/hhx9WgwYNVKVKFV1++eUaMmSIvF5vic5dEYcuJ8+nzp8uksvt0VVxCdp5+ITTSwIAIORVxOYorpAOXSdV1KHbl35SfxowSy63R3d8tEBZuSX7BwIAACiZitocxUHo2qQiD92CX1PVtNepSxh6fLtGfj8PkwAAwC4VuTmKQujapKIP3cdzfjWv152wbJfTywEAIGRV9OYoDKFrk4o+dD6fX/8Yt0Iut0fN35yun3enFX0QAAAosYreHIUhdG3C0EnpmbmK/uAnudwe3RQ/R4ePZzu9JAAAQg7NYY3QtQlDd8qGvRlq0We6XG6PHh21VF4f1+sCAFCWaA5rhK5NGLrf/W/VHvN63fdnJDu9HAAAQgrNYY3QtQlDl99b360zY3fGOh4mAQBAWaE5rBG6NmHo8svJ8+mB3x4m0ertBG05eNzpJQEAEBJoDmuErk0YuoL2p2ep9cDZcrk9un1woo5l5Tq9JAAAgh7NYY3QtQlDd27Ltx9Rs97T5HJ71O2LJPn4choAAOeF5rBG6NqEobM2fvEO83rdYXN+dXo5AAAENZrDGqFrE4bOmt/v16uTfpbL7VHTXh79tOmg00sCACBo0RzWCF2bMHSFy8r16u5hC+Rye3R1XIJ2HDrh9JIAAAhKNIc1QtcmDF3R9hzN1J8GzJLL7VHHofN0IjvP6SUBABB0aA5rhK5NGLriWbzlkC797ctpz365ki+nAQBQQjSHNULXJgxd8Y1euN38ctrHfDkNAIASoTmsEbo2YeiK78wvp7ncHs3acMDpJQEAEDRoDmuErk0YupLJyvXq3k8WnvHktGNOLwkAgKBAc1gjdG3C0JXcvvSTaj3w1JfTbhuUqPSTPDkNAICi0BzWCF2bMHSls2LHEV325qkvpz05drm8fDkNAIBC0RzWCF2bMHSl959lO83rdeOnbXR6OQAABDSawxqhaxOG7vz0nbrOjN0pq/c4vRwAAAIWzWGN0LUJQ3d+cr0+PTJyiVxuj5r3ma6fd6c5vSQAAAISzWGN0LUJQ3f+jpzI0S3vz5XL7dGf352tAxlZTi8JAICAQ3NYI3RtwtCVjeT9GbrirRlyuT26b/giZeV6nV4SAAABheawRujahKErOwnr95vX674y8Wf5/dyJAQCA02gOa4SuTRi6sjVszq9m7H6auMXp5QAAEDBoDmuErk0YurLl9/v14ter5XJ71LSXRzPX73d6SQAABASawxqhaxOGruxl5Xp132+PCb7irRnasDfD6SUBAOA4msMaoWsThs4eBzKydMO7c+Rye3Tze3OVeizb6SUBAOAomsMaoWsThs4+a/ek6/K+0+Vye9Tls8XciQEAUKHRHNYIXZswdPby/LKPOzEAACCaozCErk0YOvt9OHuzGbufzP3V6eUAAOAImsMaoWsThs5+fr9f/++3OzG43B55ftnn9JIAACh3NIe1kA7d5ORkdejQQZGRkWrYsKF69uypnJycIo/buXOn/vrXv+qPf/yjatasqTZt2mjKlCklOjdDVz6ycr164NNFcrk9atFnutbsTnN6SQAAlCuaw1rIhu7Ro0fVqFEjtWvXTgkJCRozZozq1Kmj7t27F3pcdna2WrZsqZYtW2rSpEmaNWuWnnjiCYWFhWnOnDnFPj9DV35Sj2Xr5vfmyuX2qM07s5WSdtLpJQEAUG5oDmshG7rx8fGqWbOmjhw5Yu4bOXKkIiIitHfvXsvjli5dKsMwlJiYaO7z+Xxq2rSpnnnmmWKfn6ErX5v2H1OrtxPkcnvU6cP5Op6d5/SSAAAoFzSHtZAN3ejoaHXu3DnfvrS0NIWFhWncuHGWxy1YsECGYWj16tX59l999dXq1q1bsc/P0JW/nzYd1CW9Tl2v+49xK+T1cScGAEDoozmshWzo1q9fX3369CmwPyoqSm632/K4vLw8tWrVSnfeeae2b9+utLQ0DRs2TFWrVtXy5cuLfX6GzhnjFm03v5wW9/16p5cDAIDtaA5rIRu6lSpV0qBBgwrsb9WqVZGfzB48eFBt27aVYRgyDEPVq1fX999/X+gxGRkZSklJMbekpCSGziFx3683Y3fsou1OLwcAAFsRutYI3bOcPHlS0dHR+vOf/6ypU6dq7ty5eu655xQZGamFCxdaHhcXF2eG8ZkbQ1f+vD6/nhq3Qi63R5f08mj2hgNOLwkAANsQutZCNnRLe+nCJ598omrVqunw4cP59nfo0EHR0dGWx/GJbmA5kZ2nuz5eIJfbo5Z9Z2hdSrrTSwIAwBaErrWQDd3o6Gh16dIl37709PQiv4z2/PPPq3nz5gX29+zZU1FRUcU+P0PnvP3pWbrh3TlyuT1q+85s7eW2YwCAEERzWAvZ0I2Pj1etWrWUlvb7AwRGjRpV5O3F3n//fVWtWlWHDh3Ktz82NlY333xzsc/P0AWGDXszdOVbM8zbjmVk5Tq9JAAAyhTNYS1kQ/f0AyNiYmI0c+ZMjR07VnXr1i3wwIjY2FjFxsaav969e7dq166tNm3a6L///a9mzpypp59+WoZh6Jtvvin2+Rm6wPFT8kFd2nuaXG6PHh+9TLlen9NLAgCgzNAc1kI2dCVp48aNat++vapXr64GDRqoR48eBR4BHBMTo5iYmHz7Vq1apTvvvFMNGjRQrVq11Lp1a02YMKFE52boAst/lu0078TQ49s18vu5xy4AIDTQHNZCOnSdxNAFnvemJwfF3VQAACAASURBVJux+/GcX51eDgAAZYLmsEbo2oShCzw+n18vfr3ajN3/rdrj9JIAADhvNIc1QtcmDF1gysr16qF/L5bL7dFlb07T4i2Hij4IAIAARnNYI3RtwtAFrqMncnT7oES53B5d9XaCNu7LcHpJAACUGs1hjdC1CUMX2HYdzlTrgbPkcnt0w7tzuMcuACBo0RzWCF2bMHSBb+2edF3x2z12Ow6dp/ST3GMXABB8aA5rhK5NGLrg8NOm3++x+/CIJcrO8zq9JAAASoTmsEbo2oShCx6TVuw278TQfcIq+XzcYxcAEDxoDmuErk0YuuDy4ezNZuwO+HGD08sBAKDYaA5rhK5NGLrg4vf75f7vL2bsfj5/m9NLAgCgWGgOa4SuTRi64JPn9empcSvM2J26mr87AEDgozmsEbo2YeiC08kcrzp/ukgut0fNek/T/M2pTi8JAIBC0RzWCF2bMHTB6+iJHMUOPvVAiSvfmqG1e9KdXhIAAJZoDmuErk0YuuC252im/vzubLncHrUeOEs7Dp1wekkAAJwTzWGN0LUJQxf8kvdn6Kq4BLncHkV/8JMOZmQ5vSQAAAqgOawRujZh6ELDsm2H1bzPdLncHt3x0QKengYACDg0hzVC1yYMXehIWL9fl/Q6dSeGriOWKCuXp6cBAAIHzWGN0LUJQxdaJq7YZd527OkvkpTn9Tm9JAAAJNEchSF0bcLQhZ5PE7eYsdtz8hr5/TwqGADgPJrDGqFrE4Yu9Pj9fg34cYMZu/HTNzq9JAAAaI5CELo2YehCk8/n16sTfzZj99/ztjq9JABABUdzWCN0bcLQha5cr0//HP/7o4K/Xr7L6SUBACowmsMaoWsThi60ZeV61XXEErncHl3Sy6Npa/c5vSQAQAVFc1gjdG3C0IW+jKxc3fXxArncHl325jQt+DXV6SUBACogmsMaoWsThq5iOHQ8W7cPSpTL7VHLvjO0cudRp5cEAKhgaA5rhK5NGLqKY8/RTN0YP0cut0dXxyVow94Mp5cEAKhAaA5rhK5NGLqKZcvB47p+wCy53B61HjhL21KPO70kAEAFQXNYI3RtwtBVPOtS0nVVXIJcbo9uip+jlLSTTi8JAFAB0BzWCF2bMHQV08qdR9Sy7wy53B7dNihRqceynV4SACDE0RzWHAvdzp07F3vr0qWLU8ssNYau4lrwa6qavzldLrdHnT6cr7TMHKeXBAAIYTSHNcdC97bbbivRFmwYuoptxrr9urT3NLncHt33yUIdy8p1ekkAgBBFc1jj0gWbMHSYsnqPmvY69fS0rv9eosycPKeXBAAIQTSHNULXJgwdJGnCsl3mo4IfH71MWblep5cEAAgxNIe1gAnd7OxsTZ06VR988IH69+9fYAs2DB1OG71wuxm7/xy/Qrlen9NLAgCEEJrDWkCE7s6dO9WkSRNVrlxZ4eHhioyMVFhYmMLCwhQZGakLLrjA6SWWGEOHM30y91czdrtPWCWvz+/0kgAAIYLmsBYQodu5c2e1b99eaWlpCgsL06pVq5SRkaERI0bo0ksv1dq1a0v1vsnJyerQoYMiIyPVsGFD9ezZUzk5hX8Dfty4cTIM45xbp06din1uhg5n8vv9en9Gshm7r01aIx+xCwAoAzSHtYAI3YYNG+r777+Xz+dTWFiYli5dav5syJAhateuXYnf8+jRo2rUqJHatWunhIQEjRkzRnXq1FH37t0LPS41NVVLly7Nt3355ZcyDEMfffRRsc/P0OFsfr9fcd+vN2O31/9+IXYBAOeN5rAWEKFbu3ZtzZ8/X5JUr149TZkyxfzZ3LlzFRkZWeL3jI+PV82aNXXkyBFz38iRIxUREaG9e/eW6L3i4uIUERGh/fv3F/sYhg7n4vf71XvKWjN23/punfx+YhcAUHo0h7WACN02bdroq6++kiTFxsaqY8eOyszMVHZ2tv72t7+pWbNmJX7P6Ohode7cOd++05dGjBs3rkTv1aJFC3Xs2LFExzB0sOLz+dXj2zVm7A78cQOxCwAoNZrDWkCE7pAhQ/T6669LkhYvXqxatWqpUqVKqly5siIiIjR+/PgSv2f9+vXVp0+fAvujoqLkdruL/T5JSUkyDKPEcczQoTBen18vf7PajN1/JSQTuwCAUqE5rAVE6J5t165dGjlypD7++ONSfxGtUqVKGjRoUIH9rVq1Urdu3Yr9Pq+++qqqVaumjIyMQl+XkZGhlJQUczsdyAwdrOR5fXr+PyvN2B06a7PTSwIABCFC11pAhm5ZKIvQ9fl8atSokR588MEiXxsXF3fOOzUwdChMrtenp79IMmP34zm/Or0kAECQIXStBUzo5uXlafHixZo0aZK++OKLAltJlcWlC3PmzJFhGPm+HGeFT3RRWjl5Pj01boUZu8N/2uL0kgAAQYTQtRYQobtixQo1btxY4eHh5oMiztzCw8NL/J7R0dHq0qVLvn3p6ekl+jLaU089pbp16yo7O7vE52foUBLZeV79fexyM3Y/S9zq9JIAAEGC5rAWEKHbpk0bXX/99Zo/f75SU1OVnp5eYCup+Ph41apVS2lpaea+UaNGFfv2YtnZ2apbt66eeuqpEp9bYuhQclm5Xv3fmN9jd+R8YhcAUDSaw1pAhG6NGjWUkJBQpu95+oERMTExmjlzpsaOHau6desWeGBEbGysYmNjCxw/ZcoUGYahOXPmlOr8DB1KIyvXq8dHLzNj9/P525xeEgAgwNEc1gIidG+44YZSXYdblI0bN6p9+/aqXr26GjRooB49ehR4BHBMTIxiYmIKHPvQQw+pUaNG8vl8pTo3Q4fSOpnj1aOjlhK7AIBioTmsBUTorl69Wtddd51++umnUodloGHocD7Ojl0uYwAAWKE5rAVE6NatW1fVq1dXeHi4KleurAsuuKDAFmwYOpyvkzlePTZqGbELACgUzWEtIEI3Li5O/fr1K3QLNgwdykJWbv7YHTGP2AUA5EdzWAuI0A1FDB3Kytmx+2ki99kFAPyO5rBG6NqEoUNZOvtuDDxBDQBwGs1hLSBC97rrrtOf/vSnc26tW7dWbGysXn/9dW3evNnppRYbQ4eydvZ9dofM2iy/3+/0sgAADqM5rAVE6P79739X48aNVb16dbVv315/+9vfzNuCNW7cWF26dFFUVJSqV6+uRYsWOb3cYmHoYIfsPG++xwV/MCOZ2AWACo7msBYQoTtixAi1bdtWhw4dyrc/NTVVbdu21bBhw3Ty5ElFR0erXbt2Dq2yZBg62CUnz6duXySZsfvutI3ELgBUYDSHtYAI3aZNm+r7778/58++++47NW7cWJI0efJk1ahRozyXVmoMHeyU6/Xpua9WmrEb9/16YhcAKiiaw1pAhG716tX1zTffnPNnEyZMUGRkpCQpMTGR0AV+k+f16cWvV5ux2+t/a+XzEbsAUNHQHNYCInTvuOMONWnSRMuXL8+3f+nSpWrSpInuvPNOSdLIkSN15ZVXOrHEEmPoUB68Pr9em7TGjN3XJq2Rl9gFgAqF5rAWEKG7e/duXXvttQoPD1eDBg3UqlUrNWjQQOHh4bruuuu0e/duSaeu5R0/frzDqy0ehg7lxefzq/eUtWbsvvj1auV6Q+NR2gCAotEc1gIidE/zeDzq16+fnnvuOfXr10/Tpk1zekmlxtChPPn9fvX7Yb0Zu898maTsPK/TywIAlAOaw1pAhW4oYehQ3vx+v96bnmzG7hNjlutkDrELAKGO5rDmWOhmZGSY3xLPyMgocgs2DB2c4Pf79dHsX83Y7TpiiY5l5Tq9LACAjWgOa46Fbnh4uPnls7CwMIWHhxe6BRuGDk76fP42M3bv+2Sh0jJznF4SAMAmNIc1x0J3/PjxOnz4sCRp3LhxGj9+fKFbsGHo4LT/LNuppr1OxW6nD+cr9Vi200sCANiA5rDGNbo2YegQCKas3qNLe0+Ty+3R7YMSlZJ20uklAQDKGM1hLSBC99ixYwUe/ztx4kTFxcVp/vz5Dq3q/DB0CBQz1u3XZW+eit2b4udoW+pxp5cEAChDNIe1gAjd+++/X88++6z564EDByosLEwXXnihIiIi9N///tfB1ZUOQ4dAMm9zqi7vO10ut0etB87S+r3pTi8JAFBGaA5rARG6jRo10pQpUySd+tZ4w4YN1adPH0nSq6++qrZt2zq5vFJh6BBoknYc0VVxCXK5PboqLkErdx5xekkAgDJAc1gLiNCtWrWqFixYIElKSkpSeHi4tm3bJklKTExU7dq1nVxeqTB0CETr96br+gGz5HJ71LLvDM3fnOr0kgAA54nmsBYQodukSRONGDFCktSvXz9deuml5s88Ho/q1q3r1NJKjaFDoNqaelw3xc+Ry+3RZW9O07S1+5xeEgDgPNAc1gIidHv06KHatWvroYceUmRkpPr372/+bODAgVy6AJSxlLSTun1Qolxuj5r28mjCsl1OLwkAUEo0h7WACN28vDz1799f99xzj95++23l5v7+JKcHHnhAgwcPdnB1pcPQIdAdOp6tu4ctMB8sMfynLebTCgEAwYPmsBYQoRuKGDoEg2NZuXpk5BIzdgf+uEE+H7ELAMGE5rBG6NqEoUOwyMr16ukvkszYfW3SGuV5fU4vCwBQTDSHNULXJgwdgkme16fXv11jxu5T41boZI7X6WUBAIqB5rBG6NqEoUOw8fv9enfaRjN2u3y2WGmZOU4vCwBQBJrDGqFrE4YOwWrk/K1m7HYcOk/70k86vSQAQCFoDmuErk0YOgSz/63ao2a9p8nl9uim+DnacvCY00sCAFigOaw5Frq//PJLibZgw9Ah2P206aAu7ztdLrdH1/afqZU7jzq9JADAOdAc1hwL3bCwMIWHhxe5nX5dsGHoEApW7Tqqa/vPlMvt0eV9p2vWhgNOLwkAcBaaw5pjoTtv3rwSbcGGoUOo2HLwuG5+b65cbo8u6eXR18t5ihoABBKaw1pIX6ObnJysDh06KDIyUg0bNlTPnj2Vk1O8b5EvXLhQt99+u2rUqKHatWvrlltu0ebNm4t9boYOoeRARpbu+Oj3p6gNnbWZp6gBQICgOayFbOgePXpUjRo1Urt27ZSQkKAxY8aoTp066t69e5HHzpo1S1WqVNGLL76o2bNny+Px6M0339SaNWuKfX6GDqEmIytXf/t8qRm77v/+woMlACAA0BzWAiZ0hw8frmuvvVbVq1c/57W6JRUfH6+aNWvqyJEj5r6RI0cqIiJCe/futTwuLy9PTZo0Ue/evUv1+ziNoUMoys7z6sWvV5ux+/exy3UiO8/pZQFAhUZzWAuI0P3ss89Uo0YNvfHGGwoLC9Mrr7yiV155RU2aNNEll1yiDz/8sMTvGR0drc6dO+fbl5aWprCwMI0bN87yuOnTp8swjEJjuDgYOoQqn8+vdzwbzNi9Z9hCHTyW5fSyAKDCojmsBUToXnXVVRo8eLC8Xq/CwsK0atUqSVJubq46duyot99+u8TvWb9+ffXp06fA/qioKLndbsvj3n77bV144YWaOnWqmjdvroiICLVs2VITJ04s0fkZOoS6sYu2q2mvU7F7y/tztTX1uNNLAoAKieawFhChGxkZad5ZoUqVKpo7d675sx9//FFRUVElfs9KlSpp0KBBBfa3atVK3bp1szzumWeeUbVq1VSvXj199tlnmjNnjh5//HEZhqFFixZZHpeRkaGUlBRzS0pKYugQ8mas26cWfX6/127SjiNFHwQAKFOErrWACN3GjRtr2rRpkqRLL71UQ4YMMX/21VdfqXbt2iV+z9KGbrdu3WQYhj799FNzn9/v1zXXXKN77rnH8ri4uDgZhlFgY+gQ6lbuPGLea7d5n+n68Zfzu+wHAFAyhK61gAjdRx55RP3795ck9enTRzVr1pTb7Vbfvn31hz/8QQ888ECJ37O0ly688cYbMgxDycnJ+fa//PLLuvTSSy2P4xNdVGTbUo8r+oOfzOt2/z1vK7cfA4ByQuhaC4jQTU5ONi9XyMrK0osvvqioqChdcMEFeuihh3TgQMmfxhQdHa0uXbrk25eenl7kl9HGjRt3ztB96aWXdNlllxX7/AwdKppDx7N1//BFZuy+OWUttx8DgHJAc1gLiNC1Q3x8vGrVqqW0tDRz36hRo4q8vdi+fftUqVIlffLJJ+Y+v9+vq6++Wg8++GCxz8/QoSLKyvXq2S9XmrH75NjlOs7txwDAVjSHtZAN3dMPjIiJidHMmTM1duxY1a1bt8ADI2JjYxUbG5tv38svv6yaNWtq2LBhSkhI0F//+ldFRETwwAigGM6+/dgdHy3QvvSTTi8LAEIWzWEtIEL39ttvL3IrjY0bN6p9+/aqXr26GjRooB49ehR4BHBMTIxiYmLy7cvLy1Pfvn0VFRWlKlWqqHXr1kpISCjRuRk6VHRfLNmhS367/Vjbd2ZrXUq600sCgJBEc1gLiNC9//779cADD+TbYmJiVLt2bTVu3LjAgx+CAUMHSHOTD+jKt2bI5faoZd8Zmrl+v9NLAoCQQ3NYC4jQtXLkyBHdeuutmjBhgtNLKTGGDjhlw94M3Rg/Ry63R017eTRqwTbuyAAAZYjmsBbQoStJ33//faG39QpUDB3wuwMZWbp72IJ8d2TI5Y4MAFAmaA5rAR+6EyZMUJ06dZxeRokxdEB+mTl5evqLJDN2Hx21VOmZuU4vCwCCHs1hLSBC9/vvvy+wTZ48WQMGDNAf/vAH3XfffU4vscQYOqAgr8+vd6dtNGP39sGJ2n7ohNPLAoCgRnNYC4jQDQsLO+dWpUoVde3aVQcPHnR6iSXG0AHWJq7YpWa9p8nl9uiafjO1ZOthp5cEAEGL5rAWEKG7c+fOAtuBAweC+gsrDB1QuCVbD+va/jPlcnvUrPc0TVyxy+klAUBQojmsBUTohiKGDijajkMndPvgRPNShgE/buCxwQBQQjSHNcdC91zX5Ra2BRuGDiie9MxcPT56mRm7/zdmuTKy+JIaABQXzWHNsdA9+3rc8PBwhYeHn3NfeHi4U8ssNYYOKL48r09x3683YzeWL6kBQLHRHNYcC90zr8ddsGCBXC6XXn31VS1cuFCbN2/WwoUL9corr8jlcmnBggVOLbPUGDqg5CYsy/8ltUVbDjm9JAAIeDSHtYC4RvfOO+/UgAEDzvmz/v37q1OnTuW8ovPH0AGls3TbYV3325fULu09TWMXbQ/qL6YCgN1oDmsBEbqRkZGaNWvWOX82c+ZMRUZGlvOKzh9DB5TersOZ+svQ+ealDD0nr1F2ntfpZQFAQKI5rAVE6DZp0kTPPPPMOX/29NNPq0mTJuW8ovPH0AHn53h2np758vcnqT3w6SIdzMhyelkAEHBoDmsBEbojRoxQWFiYOnTooOHDh+vbb7/V8OHD1b59e4WFhWnEiBFOL7HEGDrg/Pl8fg2dtdmM3T+/O1trdqc5vSwACCg0h7WACF1J+vHHH3XjjTeqcuXKCgsLU+XKlXXDDTfohx9+cHpppcLQAWVn+tp9uuKtGXK5PWreZ7q+Tdrt9JIAIGDQHNYCJnRP8/l8OnDggHy+4L5pPEMHlK2N+zJ06wdzzU93475fr1weLgEANEchAi50QwVDB5S9oydy9Nio3x8u8fCIJTp0PNvpZQGAo2gOa46F7quvvqrdu3eb/7uw7bXXXnNqmaXG0AH2yPP69I5ngxm7N8XP0do96U4vCwAcQ3NYcyx0mzZtqjVr1pj/u7DtkksucWqZpcbQAfb67ucUtegz3bxud/LKPU4vCQAcQXNY49IFmzB0gP3WpaTr5vd+v26379R1ysnjul0AFQvNYS2gQ9frDd4bxDN0QPk4ctZ1u50/XaQD3G8XQAVCc1gLiND98ssvNWzYMPPXGzZsUMuWLVWpUiW1b99ehw4F3/PuGTqg/Hh9fr0/I9mM3TbvzNaKHUecXhYAlAuaw1pAhO4111yjTz75xPz17bffriuvvFKffPKJWrRooW7dujm4utJh6IDyN33tPl352/12m/WepjELt8vv9zu9LACwFc1hLSBCt1atWpo7d64k6dChQ4qIiNCMGTMkSRMnTtTFF1/s5PJKhaEDnLHl4DHFDk40P93tPmGVTmTnOb0sALANzWEtIEK3Tp06mjlzpiTp22+/VbVq1ZSdferemPPnz1e1atWcXF6pMHSAc45n5+n5/6w0Y7fDkHnacvC408sCAFvQHNYCInTbtWun++67T+vXr9ett96qu+66y/zZV199pSZNmji4utJh6ABn+f1+jVqwTZf2niaX26Mr35qhaWv3Ob0sAChzNIe1gAjdRYsWqW7dugoPD1edOnWUlJRk/qxLly7q2rWrg6srHYYOCAzLth1W64GzzU93B/y4gUcHAwgpNIe1gAhdSTp27JhWrlyptLS0fPunTZumzZs3O7Sq0mPogMBxICNLXf+9xIzdLp8t1r70k04vCwDKBM1hLWBC90yZmZlB/01phg4ILHlen+KnbTRj908DZmnBr6lOLwsAzhvNYS1gQnfp0qXq1KmT6tWrp4iICK1atUqS9PLLL+t///ufw6srOYYOCEwJ6/frqrgEudweNe3l0dBZm+X1Bfc/rAFUbDSHtYAI3alTpyoiIkJ33XWXBg8erLCwMDN03333Xf3lL39xeIUlx9ABgWvn4RO686MF5qe7j45aqtRj2U4vCwBKheawFhChe9VVV6l79+6SpLy8vHyh+8MPP+iPf/yjk8srFYYOCGxZuV71nrI239PUlmw97PSyAKDEaA5rARG6VatW1Zw5cyRJXq83X+gmJiaqatWqTi6vVBg6IDh893OKrvjtaWqX9PLok7m/yselDACCCM1hLSBCt0mTJvrss88kFQzdYcOGqUWLFk4ur1QYOiB4bDl4XH8ZOt/8dPfx0cu4lAFA0KA5rAVE6L7xxhuqX7++5s2bZ4bu6tWrtXHjRl188cUaMGBAqd43OTlZHTp0UGRkpBo2bKiePXsqJyenyOMMwyiwXX755SU6N0MHBJeTOV71nLwm36UMi7cccnpZAFAkmsNaQIRudna27r77boWFhalhw4YKCwtTVFSUIiIidN999ykvr+TPqT969KgaNWqkdu3aKSEhQWPGjFGdOnXMa4ELYxiGXnnlFS1dutTc1qxZU6LzM3RAcJqyeo95KUPTXh4N4a4MAAIczWEtIEL3tDlz5qh3797q1q2b3G63Zs+eXer3io+PV82aNXXkyBFz38iRIxUREaG9e/cWeqxhGPrwww9LfW6JoQOC2dbU47rjjLsyPDxiiQ5kZDm9LAA4J5rDWkCFblmKjo5W586d8+1LS0tTWFiYxo0bV+ixhC6ArFyv+k5dl+8BE3OTDzi9LAAogOawFtChm5ubq1GjRqlly5YlPrZ+/frq06dPgf1RUVFyu92FHmsYhi688EJFREToggsu0GOPPab9+/eX6PwMHRAaPL/sMx8w4XJ71P+HDcrO8zq9LAAw0RzWHA3dbdu26b333tMLL7ygwYMH6/DhU/ewPHnypP71r3+pUaNGCgsLU3R0dInfu1KlSho0aFCB/a1atVK3bt0KPfbJJ5/U5MmTNX/+fA0bNkz169dXixYtlJmZaXlMRkaGUlJSzC0pKYmhA0LE7iOZeuDTRWbs3j1sgbYfOuH0sgBAEqFbGMdCd/HixapRo4bCwsLM7ZJLLtHKlSvVrFkzhYWFqX379po3b16p3v98QvdsS5YskWEYGj16tOVr4uLiznm3BoYOCA25Xp/en5Gspr1Oxe4Vb83Q5JV75PfzRTUAziJ0rTkWuu3bt1eLFi20bNkyZWVladOmTerQoYNq1KihP/zhD0pISDiv9z+fSxfOJSoqSs8//7zlz/lEF6gYFvyaqtYDZ5uf7v6/r1crIyvX6WUBqMAIXWuOhW6DBg00adKkfPu2bt2qsLAwjR8//rzfPzo6Wl26dMm3Lz09vVhfRjuXqKgovfDCC8V+PUMHhK7UY9l6cuxyM3ZveX+uVu48UvSBAGADmsOaY6EbFham5cuX59t3+mERSUlJ5/3+8fHxqlWrltLS0sx9o0aNKtbtxc62aNEiGYahMWPGFPsYhg4IbX6/X2MWblfzN6fL5fbo0t7T9NHsX5Xn9Tm9NAAVDM1hzdHQXbFiRb59Zz4V7XydfmBETEyMZs6cqbFjx6pu3boFHhgRGxur2NhY89cjR47UP//5T02cOFE//fSTPvzwQ1144YVq2bKlTp48WezzM3RAxbBhb4ZiByean+4++Nli7T5i/cVVAChrNIc1R0O3du3auuCCC/JtVvtLY+PGjWrfvr2qV6+uBg0aqEePHgUeARwTE6OYmBjz13PmzNHNN9+sevXqqVKlSoqKilK3bt2UmppaonMzdEDFcTLHq95T1pqxe9XbCZq6mv/bB1A+aA5rjoVuv379SrQFG4YOqHhmrt+v6/rPzPdFtfSTfFENgL1oDmsB/cCIYMbQARXTwYwsPT56mRm7N783V0u3HXZ6WQBCGM1hjdC1CUMHVFw+329fVOtz6otqTXt5FD9tI09UA2ALmsMaoWsThg5A8v4Mdfpwvvnp7h0fLdCm/cecXhaAEENzWCN0bcLQAZCk7Dyv4qdtNJ+o1vzN6Rq1YJt8Pp6oBqBs0BzWCF2bMHQAzrR022Hd/N5c89PdR0Yu0Z6j3IYMwPmjOawRujZh6ACcLSMrV69O/NmM3VZvJ+jbpN3y+/l0F0Dp0RzWCF2bMHQArExfuy/fbcie/iJJh45nO70sAEGK5rBG6NqEoQNQmIPHsvTP8SvM2P3TgFmasW6f08sCEIRoDmuErk0YOgBF8fv9mrRit1q9nWAG70vfrFZaZk7RBwPAb2gOa4SuTRg6AMW1+0im/jpyqRm7bd+ZrbnJB5xeFoAgQXNYI3RtwtABKAmfz69xi7br8r7TzeDt8e0aZWTxCGEAhaM5rBG6NmHoAJTG9kMn1PnTShyziwAAIABJREFURWbs3hg/R4mbDjq9LAABjOawRujahKEDUFpen18j5m01HyHscnvUczKf7gI4N5rDGqFrE4YOwPnacvCY7h/Op7sACkdzWCN0bcLQASgL5/p0t8e3a5Seyae7AE6hOawRujZh6ACUpS0Hj+uBM67dbfvObM1cv9/pZQEIADSHNULXJgwdgLLm9fk1asG2fHdmePHr1TrMU9WACo3msEbo2oShA2CXHYdO6OERS/I9Ve27n1Pk9/udXhoAB9Ac1ghdmzB0AOzk8/n15dKduvKtGWbwPjVuhfaln3R6aQDKGc1hjdC1CUMHoDykpJ3U/41ZbsZuq7cT9OXSnfL5+HQXqChoDmuErk0YOgDlxe/3a8rqPbq2/0wzeLv+e4m2HDzu9NIAlAOawxqhaxOGDkB5O3Q8Wy9+vdqM3eZvTtfHc35VTp7P6aUBsBHNYY3QtQlDB8Apszcc0I3xc8zg7TBknlbuPOL0sgDYhOawRujahKED4KRjWbl6+7t1atrrVOw27eVR36nreIwwEIJoDmuErk0YOgCBYNWuo/rL0Pnmp7t/fne2pq/dx63IgBBCc1gjdG3C0AEIFDl5Pn0y99d8jxF+atwK7Tma6fTSAJQBmsMaoWsThg5AoNlx6IQeG7XMjN2WfWfo8/nblOfly2pAMKM5rBG6NmHoAAQiv9+vqatTdP2AWWbw3vHRAq3addTppQEoJZrDGqFrE4YOQCBLy8zRG5N/MWO3aS+Pev1vrdIyc5xeGoASojmsEbo2YegABIOkHUfyfVnt+gGz9N+Ve/iyGhBEaA5rhK5NGDoAwSLX69OIeVvVsu8MM3gfHrFEmw8cc3ppAIqB5rBG6NqEoQMQbPYczdTTXySZsdus9zS9O22jTmTnOb00AIWgOawRujZh6AAEqzkbD+iW9+eawXtj/BxN4967QMCiOawRujZh6AAEs5M5Xg2euUnN3/z93ruPj16mLQePO700AGehOayFdOgmJyerQ4cOioyMVMOGDdWzZ0/l5JTsG8Uvv/yyDMPQyy+/XKLjGDoAoWBb6nE9Pvr3e+9e9uY0xU/ncgYgkNAc1kI2dI8ePapGjRqpXbt2SkhI0JgxY1SnTh1179692O+xdu1a1apVS7Vr1yZ0AVRYfr9f09fu003xc8zgveHdOfphzV4uZwACAM1hLWRDNz4+XjVr1tSRI0fMfSNHjlRERIT27t1brPdo166d3n77bblcLkIXQIWXmZOnQQn5L2d4ZOQSJe/PcHppQIVGc1gL2dCNjo5W586d8+1LS0tTWFiYxo0bV+Tx//nPf9S4cWNlZmYSugBwhu2HTujJscvN2L2kl0dvf7dO6Zm5Ti8NqJBoDmshG7r169dXnz59CuyPioqS2+0u9Nhjx46pUaNGmjhxoiQRugBwFr/fr9kbDij6g5/M4L2u/0xNWLZLXh+XMwDlieawFrKhW6lSJQ0aNKjA/latWqlbt26FHvvaa68pOjra/HVxQjcjI0MpKSnmlpSUxNABCHlZuV4N/2lLvodN3PnRAi3ffqTogwGUCULXGqF7lvXr16tq1ar6+eefzX3FCd24uDgZhlFgY+gAVAR7006q+4RVZuy63B51n7BKKWknnV4aEPIIXWshG7qlvXShU6dOevzxx5WWlmZujRs31nPPPae0tDT5fL5zHscnugAgLd9+RHd9vMCM3cv7TtfQWZt1Msfr9NKAkEXoWgvZ0I2OjlaXLl3y7UtPTy/yy2gul+ucn8ye3nbs2FGs8zN0ACoqr8+vr5fv0p8GzMr3dLWpq1Pk4/pdoMzRHNZCNnTj4+NVq1YtpaWlmftGjRpV5O3Fli5dqsTExHxbw4YN9eCDDyoxMVFZWVnFOj9DB6CiSz+Zq4E/blCz3tPM4L1/+CKt2nXU6aUBIYXmsBayoXv6gRExMTGaOXOmxo4dq7p16xZ4YERsbKxiY2MLfS/uugAApbct9bj+OT4p3/W7/+/r1Vy/C5QRmsNayIauJG3cuFHt27dX9erV1aBBA/Xo0aPAI4BjYmIUExNT6PsQugBw/hb+ekh/GTrfjN0WfabrgxnJOpbF/XeB80FzWAvp0HUSQwcABeV5ffrPsp1qPfD363dbD5ylCct2Kc977i/7AigczWGN0LUJQwcA1o5l5eqDGclq3uf3xwl3HDpPP206KL+fL6wBJUFzWCN0bcLQAUDR9hzN1EvfrM53/e5jo5Zp/d50p5cGBA2awxqhaxOGDgCK7+fdaer67yVm7Dbt5dFrk9ZoXzpfWAOKQnNYI3RtwtABQMn4/X4lrN+v2wYl5vvC2vszkpV+ki+sAVZoDmuErk0YOgAonVyvT+MX78j3wIlr+8/UqAXblJ3HE9aAs9Ec/7+9e4+Lot7/Bz6rCCpLoKYoCYoKWQKWhXoy5XuUk3U81Q89lmV+6WupFX47pZYXVL6VgWmpmXnJC6l5Qcy8oVzkoqUk6yW1REVURPMKuytyh339/vAwx3EZRN3Zy/B6Ph7zeBw+zc58due9h5fDZz4feQy6CmHRERE9GGNpBWYlZuPRqf95YK3PzFSusEZ0B2YOeQy6CmHRERFZxiVDKSZuPALfSf95YO35eXuQzhkaiAAwc9SFQVchLDoiIss6dfmG2Qprry7Zh8Pn9Xd/MZGKMXPIY9BVCIuOiEgZurMFGLJwryTwjll1ADlXbti6a0Q2wcwhj0FXISw6IiLlmEwmpPxxGX+bkyGGXd9J2zF+w2/ILyy2dfeIrIqZQx6DrkJYdEREyquqNiH+QD6eiUkVA6/flB2I2vI7rhWV2bp7RFbBzCGPQVchLDoiIuspq6xC7C9n8NRn/5mS7LFpOzErMRuGYs7BS+rGzCGPQVchLDoiIuu7WVaJb1JPIWB6ohh4A6MSsSAtBzfLKm3dPSJFMHPIY9BVCIuOiMh29MXlmLkzG12n7hQDb49Pk7F0Ty5KK7joBKkLM4c8Bl2FsOiIiGzvyo1SRG35HX5T/rPoRM/PU7By31muskaqwcwhj0FXISw6IiL7cUFfgokbj6DT5AQx8D4Tk4q1+/NQUVVt6+4RPRBmDnkMugph0RER2Z9z129iXNxvklXWnv0iFXFZ5xl4yWExc8hj0FUIi46IyH6dvlqE/117CB1vC7z9ZqVhg+48Khl4ycEwc8hj0FUIi46IyP6dvHwD7605KFllLWRWGuIP5DPwksNg5pDHoKsQFh0RkePIvmTEO6sPSAJvv1lpiNNxSAPZP2YOeQy6CmHRERE5nj8uGjF6lU4SePt+kYb1WXxojewXM4c8Bl2FsOiIiBzX8T/N7/A+E5OK1ZnnOC0Z2R1mDnkMugph0REROb4Tl26N4b39obVen+9C7C9nuPAE2Q1mDnkMugph0RERqcepyzfw/rpDkmnJnp6Rgu9253JpYbI5Zg55DLoKYdEREanPmWs3MWHDb+h828IT3T9Jwte7TsFQUmHr7lEDxcwhj0FXISw6IiL1Ol9QjMmbjkqWFu42PREzd2bjWlGZrbtHDQwzhzwGXYWw6IiI1O+SoRSfbvsDXafuFAOvf+QOTN98DPmFxbbuHjUQzBzyGHQVwqIjImo4rheVYXbiCQRMTxQDb+fJCfgw7jBOXb5h6+6RyjFzyGPQVQiLjoio4TGWVuDb9Bw89VmyZGqyt1fqcDCv0NbdI5Vi5pDHoKsQFh0RUcNVWlGFVfvO4pmYVEngHbp4H9Kyr8BkMtm6i6QizBzyGHQVwqIjIqKKqmpsOpSP5+bslgTegXN3Y9OhfK62RhbBzCGPQVchLDoiIqphMpmQln0FQxftkwTev0TvwtI9uSjiXLz0AJg55DHoKoRFR0REtTlwrhCjVuokq60FRt2amuyKsdTW3SMHxMwhT9VBNzs7G6GhoWjevDk8PT3x0Ucfoby8vM7XXLp0CS+99BK8vb3h4uKCtm3bYsiQIThx4sQ9nZtFR0REdTl9tQiTfjwimYvXb8oOTNjwG05ypga6B8wc8lQbdAsLC9GuXTv069cPiYmJWL58Odzd3REREVHn63JzcxEeHo7vv/8e6enpWLduHQIDA+Hl5YWCgoJ6n59FR0RE9XHlRilmJWYjMCpRMqzhv5fvxy851/jgGt0VM4c81Qbd6OhoaLVaSThdsmQJGjdujIsXL97TsU6dOgVBEBAXF1fv17DoiIjoXtwsq0TsL2fw7BfSmRqen7cHGw/ko7ySD65R7Zg55Kk26Pbt2xdhYWGSNr1eD41Gg9jY2Hs6VkFBAQRBwOrVq+v9GhYdERHdj8qqamw/8ideWvCLJPAGz0jBgrQcFN6sewgeNTzMHPJUG3Rbt26NyMhIs3YvLy9MnDjxrq+vrq5GZWUl8vLy8MYbb6BDhw4wGo31Pj+LjoiIHoTJZELW2QKMWXVA8uDao1N3YMqmo8i5wnG8dAszhzzVBl0nJyfMnj3brL1bt24YNWrUXV8/ZswYCIIAQRDQqVMnnDx5ss79jUYjLly4IG46nY5FR0REFnHu+k1Ebfkdj03baTaON+PkVY7jbeAYdOUx6MrIy8tDVlYWNm3ahJCQEPj4+CA/P192/6ioKDEY376x6IiIyFIMJRVYsvu02YprA77KwOrMcygu53y8DRGDrjzVBt0HHbpwu5KSEnh5eWHs2LGy+/COLhERWUtlVTUSjv6JIQv3SgJvYFQiPk84jvMFxbbuIlkRg6481Qbdvn37YvDgwZI2g8FwXw+jAcCAAQMwcODAeu/PoiMiIms4fF6Pf607hC5TEsTA6ztpO0av0mHvaU5P1hAwc8hTbdCNjo6Gm5sb9Hq92LZ06dL7ml7sxo0baNOmDd555516v4ZFR0RE1nTFWIqvkk/iqc+SJXd5n5uzGz/8ymENasbMIU+1QbdmwYiQkBAkJSVhxYoV8PDwMFswon///ujfv7/485dffomIiAjExcUhIyMDq1atwtNPPw03N7d7Wh2NRUdERLZQVlmFjQfy8eI3P0sCb0BUIj7d9gfOXrtp6y6ShTFzyFNt0AWA48ePY8CAAWjWrBnatGmDCRMmmC0BHBISgpCQEPHnlJQU9O/fHw8//DBcXFzQqVMnjBgxAjk5Ofd0bhYdERHZkslkwsG8Qrx/x7CGmtkadh2/jKpqDmtQA2YOeaoOurbEoiMiIntx5UYp5qWcQvCMFEng7TMzFYsyTqOAi1A4NGYOeQy6CmHRERGRvan496prryzeJwm8fpE78MH6wzhwroAPrzkgZg55DLoKYdEREZE9O3HpBiJ/OorH71iE4vl5e/DDr+dQVMaH1xwFM4c8Bl2FsOiIiMgRFJVVYlXmOQycu1sSeB+fthNTNh3FHxeNtu4i3QUzhzwGXYWw6IiIyJGYTCbozhbg/XWH4DdlhyT0vrzgF2zQnUdJeZWtu0m1YOaQx6CrEBYdERE5qutFZViccRr9ZqWZTVE2ffMxZF/iXV57wswhj0FXISw6IiJydNXVJuw5dRWjV+nQabJ0irL/9+0viNOd50IUdoCZQx6DrkJYdEREpCZXjKVYkJaDPjNTpXd5pydiyqajOJpvsHUXGyxmDnkMugph0RERkRpVV5uw++RVjFl1AJ3vuMv796/3YNW+szCUVNi6mw0KM4c8Bl2FsOiIiEjtrt4ow6KM0/iv2emSwOv/73l5952+znl5rYCZQx6DrkJYdERE1FCYTCZk5l7Hv9Ydgn+kdMaGfrPSsCAtB5cMpbbupmoxc8hj0FUIi46IiBoiQ3EFVu07i0Hz90gCr++k7fjv5fux/cifKKvkNGWWxMwhj0FXISw6IiJq6H6/aMD0zccQ9H9JktDb/ZMkTN98DMcuGDi0wQKYOeQx6CqERUdERHRLaUUVtv52ESOW70fHSdsloXfg3N1YuicXV2+U2bqbDouZQx6DrkJYdEREROYu6kswf9cphNyxGEWnyQkYGZuFhKMc2nCvmDnkMegqhEVHREQkz2QyIetsAT6OP4Ju0xMloTfo/5IwZdNRHMwr5NCGemDmkMegqxAWHRERUf0Ul1fix4P5GL70V7OhDf81Ox1f7zqF8wXFtu6m3WLmkMegqxAWHRER0b27qC/Bt+k56P+ldG7eDhO345+L9mLNr3kwFHNBitsxc8hj0FUIi46IiOj+mUwm/HZej6gtv+PJT5Mlgddvyg6MWXUAO49d4nheMHPUhUFXISw6IiIiy6ioqkbKH5fx3g8H4XfHghSBUYmY9OMRZOZeR3V1wxzPy8whj0FXISw6IiIiyzOUVCAu6zyGLck0G8/bO3oXohOO4/eLDWt+XmYOeQy6CmHRERERKetPQwkWZ5zGwLm7zcbzhn6Vgfm7TuHc9Zu27qbimDnkMegqhEVHRERkPScu3cCsxGz0mZlqFnpf+uZnLN2Ti8vGUlt3UxHMHPIYdBXCoiMiIrI+k8mEA+cKMW3zMfS44yG2jpO245XF+7A68xyuF6lnJTZmDnkMugph0REREdlWZVU1dp+8ivEbfkPAHYtSdJqcgDeW/Yq4rPMOP10ZM4c8Bl2FsOiIiIjsR2lFFXYeu4SINQfRdepOSejtMiUB/xObhY0H8mEsdbzQy8whj0FXISw6IiIi+1RcXomtv13E6FU6s+nK/KbswFvf67DpkOOEXmYOeQy6CmHRERER2b8bpRXYdCgfb32fhS5TEmoJvVn48WA+DCX2G3qZOeQx6CqERUdERORYDCUV2HggH2+u2G8WemuGN8TpzkNfXG7rrkowc8hj0FUIi46IiMhxGYpvhd6RsVnwmyId3tD53w+yrfk1D9fsYPYGZg55DLoKYdERERGpg7G0Aj8ezMfbK83H9PpO2o6hi/dh+c9ncEFfYpP+MXPIY9BVCIuOiIhIfYrKbj3I9u4PB8xmb6hZnGJBWg5yrhRZrU/MHPIYdBXCoiMiIlK3kvJbU5Z9sP4wAqISzUJv/y/T8cXObPx2Xg+TyaRYP5g55DHoKoRFR0RE1HCUV1Yj4+RVTPrxCJ76LNks9PaO3qXYwhTMHPIYdBXCoiMiImqYqqpNyDpbgM+2/YFnv0hFh4nb8dyc3Yqdj5lDnqqDbnZ2NkJDQ9G8eXN4enrio48+Qnl53VOCnDhxAhEREXjsscfQrFkzdOzYEe+++y6uX79+T+dm0REREZHJZMLxP43Ye/qaYudg5pCn2qBbWFiIdu3aoV+/fkhMTMTy5cvh7u6OiIiIOl/3zTffICgoCPPmzUN6ejq+//57+Pj4oFu3bncNybdj0REREZE1MHPIU23QjY6OhlarRUFBgdi2ZMkSNG7cGBcvXpR93fXr180GjO/duxeCIGDz5s31Pj+LjoiIiKyBmUOeaoNu3759ERYWJmnT6/XQaDSIjY29p2OVlJRAEAQsXLiw3q9h0REREZE1MHPIU23Qbd26NSIjI83avby8MHHixHs6VnJyMgRBQHp6er1fw6IjIiIia2DmkKfaoOvk5ITZs2ebtXfr1g2jRo2q93FKS0sRFBSEp556qs458IxGIy5cuCBuOp2ORUdERESKY9CVx6B7F+Hh4XBzc8Pvv/9e535RUVEQBMFsY9ERERGRkhh05ak26Fpi6EJkZCScnJyQmJh41315R5eIiIhsgUFXnmqDbt++fTF48GBJm8FgqPfDaPPnz4dGo8HKlSvv6/wsOiIiIrIGZg55qg260dHRcHNzg16vF9uWLl161+nFAGDt2rXQaDSYOXPmfZ+fRUdERETWwMwhT7VBt2bBiJCQECQlJWHFihXw8PAwWzCif//+6N+/v/hzRkYGmjRpggEDBiAzM1Oy5efn1/v8LDoiIiKyBmYOeaoNugBw/PhxDBgwAM2aNUObNm0wYcIEs9XNQkJCEBISIv4s91CZIAiIioqq97lZdERERGQNzBzyVB10bYlFR0RERNbAzCGPQVchLDoiIiKyBmYOeQy6CmHRERERkTUwc8hj0FUIi46IiIisgZlDHoOuQvLy8iAIAnQ6nWQhCW7cuHHjxo0bN0tuNYtU5eXl2Tr+2B0GXYXUFB03bty4cePGjZs1Np1OZ+v4Y3cYdBVSXl4OnU6HvLw8xf8Fx7vG9vEvaV4H22+8Fvax8TrYx8brYD+b0tciLy8POp3ObApVYtB1aBcucEyOPeB1sB+8FvaB18E+8DrYD14L22HQdWD84tgHXgf7wWthH3gd7AOvg/3gtbAdBl0Hxi+OfeB1sB+8FvaB18E+8DrYD14L22HQdWBGoxFRUVEwGo227kqDxutgP3gt7AOvg33gdbAfvBa2w6BLRERERKrEoEtEREREqsSgS0RERESqxKBLRERERKrEoEtEREREqsSg64Cys7MRGhqK5s2bw9PTEx999BFXQ7GgDRs24KWXXsIjjzyC5s2bo3v37oiNjYXJZJLst2zZMvj5+cHFxQVBQUHYtm2b2bEMBgNGjhyJFi1aQKvVYsiQIfjzzz+t9VZUpaioCI888ggEQcDhw4cl/43XQnlVVVWYNWsW/P394ezsDC8vL4wZM0ayj8lkQkxMDLy9vdG0aVP07t0bmZmZZse6ePEiBg8eDK1WixYtWuCtt97i0+j1tGXLFvTs2RNarRZt27bFK6+8gjNnzpjtx++E5eTk5GD06NEIDAxEo0aNEBISUut+lvzM9+7di969e6Np06bw8fHBzJkzzX4HUf0w6DqYwsJCtGvXDv369UNiYiKWL18Od3d3RERE2LprqtG7d28MGzYM69evR2pqKiZNmoRGjRphxowZ4j7r1q2DRqPB1KlTkZaWhjFjxsDJycnsl/rAgQPRvn17xMXFYcuWLQgICED37t1RWVlp7bfl8D7++GN4enqaBV1eC+sIDw9Hu3btsHDhQmRkZGDNmjUYP368ZJ+YmBg4Oztjzpw52LVrF8LCwuDm5obc3Fxxn4qKCgQEBCAgIABbt27F+vXr0b59ewwaNMjab8nhpKeno1GjRnjzzTeRkpKC9evXw9/fH/7+/igtLRX343fCsjZv3gxvb28MHToU/v7+tQZdS37mOTk50Gq1CAsLw65duzBnzhw4Oztj9uzZSr9VVWLQdTDR0dHQarUoKCgQ25YsWYLGjRvj4sWLNuyZely7ds2sbdSoUWjRooX4s7+/P15//XXJPn/5y1/wwgsviD/v27cPgiAgOTlZbDtx4gQ0Gg3i4uIU6Ll6ZWdnw9XVFYsXLzYLurwWyktKSoKTkxP++OMP2X1KS0vx0EMPYcqUKWJbeXk5OnTogHfffVdsW7t2LTQaDU6cOCE5viAI2L9/vzJvQCXGjBkDX19fyZ29tLQ0CIKAffv2iW38TlhWdXW1+L9ffvnlWoOuJT/z0aNHo2PHjpK/1E6ePBkeHh4oKyuzxFtqUBh0HUzfvn0RFhYmadPr9dBoNIiNjbVNpxqAhQsXQhAElJSUIDc3F4IgYMuWLZJ9vv76azg7O4v/RzRt2jS0bNnS7M9NTz75JMLDw63VdVUIDQ3F+PHjkZ6eLgm6vBbW8corr+C5556rc5/U1FQIgoAjR45I2j/88EN06NBB/HnEiBF48sknJfuYTCa0bNkSUVFRluqyKo0cORJBQUGStoMHD0IQBOzduxcAvxNKqy3oWvoz9/b2xocffijZ58iRIxAEAenp6RZ7Lw0Fg66Dad26NSIjI83avby8MHHiRBv0qGF4/fXXxV/WCQkJEAQBOTk5kn2Sk5MhCAKys7MBAEOHDkWfPn1qPVavXr0U77NaxMfHw9PTE0aj0Szo8lpYh4+PD8aOHYv3338fDz30EJo2bYpBgwZJxoZ+++230Gg0Zs8LfPfdd9BoNCgpKQEABAcHY/jw4WbneOaZZ/Dqq68q+0Yc3J49e+Dk5IRvv/0WBoMBubm5GDhwIJ5++mnxriO/E8qqLeha8jO/efMmBEHA0qVLJfuUl5dDo9Fg0aJFFnw3DQODroNxcnKqdZxOt27dMGrUKBv0SP1+/vlnNGrUCPPnzwcA/PDDDxAEwWyIg06nk9xZCQ0NrXXcYUREBPz8/JTvuAoUFxfD29sby5cvBwCzoMtrYR3Ozs7QarXo3bs3EhISsGHDBnTu3BmPPfaYOLZwxowZcHV1NXttfHw8BEEQh1Z16dKl1mcKBg0ahL/97W/KvhEV2Lp1K7RaLQRBgCAIePLJJ3HlyhXxv/M7oazagq4lP/MLFy5AEATEx8eb7efq6orPP//cQu+k4WDQdTAMutaVn58PLy8vDBgwQLxjwl8k1jN58mQ8/fTT4p/6GHRto0mTJmjevLnkc675k/nGjRsBMOhawy+//AIPDw+MGzcOaWlpiI+PR1BQEIKDg8U/j/M7oSwGXcfDoOtgOHTBevR6PQICAhAYGAiDwSC280+D1nHu3Dk4OzsjISEBer0eer0e27ZtgyAI2LNnD4qKingtrKRNmza1fk7u7u749NNPAXDogjU89dRTGDp0qKQtPz8fGo1G/KsHvxPK4tAFx8Og62D69u2LwYMHS9oMBgMfRrOwkpIS9OnTB97e3rhw4YLkv9U8eLB161ZJ+/z58+Hs7Cz+op82bRpatWplduwePXrwYY96qLl7K7eFhITwWlhJSEiIbNCtmXav5mG0o0ePSvYZN26c2cNoPXr0kOxjMpnQqlUrPox2F82aNav1jl7r1q3F2S74nVBWXQ+jWeoz9/b2xrhx4yT7HD16lA+j3ScGXQcTHR0NNzc36PV6sW3p0qWcXsyCKisr8Y9//AMtW7aUnU7J398fb7zxhqStT58+tU4ls2vXLrHt5MmTnL6nnvR6PdLT0yXb3LlzxbsdNcMXeC2UN3v2bDRr1gxXr14V22r+LFszKX7N9GJTp04V96moqEDHjh1rnV7s1KlTYltKSgqnF6uHrl27mt3oOHfuHDQaDb777juxjd8J5dQ1vZilPvPRo0fD19cXFRUVYltkZCQ8PDy4ONR9YNB1MDULRoSEhCApKQkrVqyAh4cHF4ywoFGjRkEAO2nBAAAIHElEQVQQBHz11VfIzMyUbDXj4Gp+WU+fPh3p6el455134OTkJJnLErg1Obi3tzc2bNiArVu3IjAwkBOyP4A7x+gCvBbWYDQa4ePjg549e2LLli1Yt24dfH19ERwcLJkqKSYmBi4uLpg3bx5SU1MxZMgQ2QUjAgMDsW3bNsTFxcHb25sLRtTDvHnzIAgC3n//fXHBiICAALRr1w6FhYXifvxOWFZxcTHi4+MRHx+P4OBgPP744+LPNf/4s+RnnpOTA1dXVwwZMgSpqamYN28eF4x4AAy6Duj48eMYMGAAmjVrhjZt2mDChAn8V54FdejQQfbP5WfPnhX3W7ZsGbp06QJnZ2fxl/adapZ79PDwgFarxeDBg3nn/QHUFnQBXgtrOH36NP7+97/D1dUV7u7ueO2113D58mXJPiaTCdHR0Wjfvj1cXFzQq1cvs1/0wK0HbmqWAPbw8MDIkSO5BHA9mEwmLFq0CEFBQXB1dUXbtm0RFhYmWXyjBr8TlnP27FnZ3wm3DyWw5Ge+d+9e9OrVCy4uLmjfvj1iYmK4BPB9YtAlIiIiIlVi0CUiIiIiVWLQJSIiIiJVYtAlIiIiIlVi0CUiIiIiVWLQJSIiIiJVYtAlIiIiIlVi0CUiIiIiVWLQJSIiIiJVYtAlogYhKipKXM1Io9HA3d0dTzzxBMaPHy9Z8c7SSktLERUVZbaaW81qSz/99NN9HffSpUvQarU4fvy4JbpZLwMHDsQnn3xitfMRET0oBl0iahCioqKg1WqRmZmJzMxMJCUlISYmBj4+PtBqtUhKSlLkvHq9HoIgIDY2VtL+oEF37NixGDx4sAV6WH/p6elwd3dHYWGhVc9LRHS/GHSJqEGIioqCu7u7Wbter0dgYCBatGgBo9Fo8fMqEXSLiorg6uqKHTt2WKiX9efr64u5c+da/bxERPeDQZeIGgS5oAsAO3fuhCAIWLx4sdhWXV2NmTNnokuXLnB2dkbnzp0l/x0AwsPD0b17d2zduhVdu3aFi4sLevbsiYMHD4r71AyXuH07e/asGHTXrFmDt99+Gw899BDat2+Pzz///K7vJTY2Fu7u7qioqBDb5ILzyy+/jJCQELPPISsrC8HBwWjatCmCg4Nx/PhxGI1GDB8+HG5ubvD19cXatWvNzv3xxx/jiSeeuGsfiYjsAYMuETUIdQXd0tJSODk5ITw8XGx777330Lx5c0RHRyMlJQVTp05F48aNsW7dOnGf8PBwtGrVCr6+vli1ahV++ukn9OjRA61atYLBYAAApKSkQBAETJ06VRw2UVZWJgZTHx8fjB8/HsnJyfjggw8gCAISEhLqfC/Dhw9HaGiopO1egm7Tpk0REBCAlStXYuvWrfDz80NQUBBefPFFTJ8+HcnJyRg2bBicnJyQl5cnOd7mzZshCAKuXbtWZx+JiOwBgy4RNQh1BV0AaNu2LZ5//nkAQE5ODjQaDZYvXy7Z591334Wfn5/4c3h4OARBwO7du8W2K1euwMXFBTExMQDuPnRh+PDhYpvJZIKvry/eeuutOt+Lv78/Pvjgg1qPV5+gKwgCkpOTxbY1a9ZAEARERESIbUajEU2aNME333xT63lsMWyCiOheMegSUYNwt6Dr6ekpBt3FixejUaNG0Ov1qKysFLeNGzdCEATxbm3NHd079e/fH2FhYQDuHnTvbH/hhRcwcODAOt+Lm5sbZsyYUevx6hN0GzdujMrKSrFt//79tb7W29sbEydOlLQVFRVBEASsWLGizj4SEdkDBl0iahDqM3ThzTffBADMmDGj1rG1NduxY8cA3Aq6Xbt2NTvesGHD0Lt3bwD3/jDancG0Nk2aNMHs2bPv63i1fQ6HDx+GIAhIT0+XtHfu3Bn/+te/JG0VFRUQBAELFy6ss49ERPaAQZeIGoS6gu6OHTsgCAKWLFkCAFi4cCEaNWqEzMxM6HQ6s62kpASAZe7o3k/Q9fT0RGRkZK3Hi4+Pl7SHhoZaNOhevXoVgiBg/fr1dfaRiMgeMOgSUYMgF3QNBoPZ9GInT56ERqO56zjU+ozRLS4uhiAIWLRokeS1DxJ0+/Xrh9dff73W43322WdiW3l5OTw9PS0adLOysiAIAg4dOlRnH4mI7AGDLhE1CHcuGJGcnIyZM2eiQ4cOtS4YMXbsWLRo0UKcdSEhIQFffvklXnvtNXGf22ddWL16tTjrQsuWLcVxvMCtsa5//etf8fPPP0On06G8vPyBgu7kyZPx6KOPStpqjvfwww9j2bJlSEpKwj//+U80adIEPj4+SEtLEz+HBwm6CxcuhFarRVVVVZ19JCKyBwy6RNQg3LkEsJubG7p3745x48bVugSwyWTC/Pnz0a1bNzg7O6NVq1Z49tlnJXPp1syju3nzZjz66KNwdnZGcHAwdDqd5Fg7duxAQEAAXFxczObRvZ+ge+DAAQiCgDNnzohtNcebNm0aunbtiqZNm2LIkCHYuHEjtFotXn31VfFzeJCg++KLL2LEiBF19o+IyF4w6BIR3aeaoGsLTzzxhGRxiQddUrg+CgsL4ezsjIyMDMXOQURkSQy6RET3yZZB98cff8QjjzyC8vJyANYJujNmzEC/fv0UOz4RkaUx6BIR3SdbBl0AmDVrFnJzcwFYJ+guWLAAR44cUez4RESWxqBLRERERKrEoEtEREREqsSgS0RERESqxKBLRERERKrEoEtEREREqsSgS0RERESqxKBLRERERKrEoEtEREREqsSgS0RERESqxKBLRERERKrEoEtEREREqvT/AeON0+6lfci6AAAAAElFTkSuQmCC\" width=\"639.8333142648146\">"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"Text(0.5,1,'Silicon @ 17.8 keV')"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"depth = numpy.linspace(0, 1000, 100)\n",
"res = numpy.exp(-mu*depth*1e-6)\n",
"fig, ax = subplots()\n",
"ax.plot(depth, res, \"-\")\n",
"ax.set_xlabel(\"Depth (µm)\")\n",
"ax.set_ylabel(\"Residual signal\")\n",
"ax.set_title(\"Silicon @ 17.8 keV\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is consistent with:\n",
"http://henke.lbl.gov/optical_constants/filter2.html\n",
"\n",
"Now we can model the detector\n",
"\n",
"## Modeling of the detector:\n",
"\n",
"The detector is seen as a 2D array of voxel. Let vox, voy and voz be the dimention of the detector in the three dimentions.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Detector Pilatus 1M\t PixelSize= 1.720e-04, 1.720e-04 m\n",
"0.000172 0.000172 0.00045\n"
]
}
],
"source": [
"detector= pyFAI.detector_factory(calib[\"detector\"])\n",
"print(detector)\n",
"\n",
"vox = detector.pixel2 # this is not a typo\n",
"voy = detector.pixel1 # x <--> axis 2\n",
"voz = thickness\n",
"\n",
"print(vox, voy, voz)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The intensity grabbed in this voxel is the triple integral of the absorbed signal coming from this pixel or from the neighboring ones.\n",
"\n",
"There are 3 ways to perform this intergral:\n",
"* Volumetric analytic integral. Looks feasible with a change of variable in the depth\n",
"* Slice per slice, the remaining intensity depand on the incidence angle + pixel splitting between neighbooring pixels\n",
"* raytracing: the decay can be solved analytically for each ray, one has to throw many ray to average out the signal.\n",
"\n",
"For sake of simplicity, this integral will be calculated numerically using this raytracing algorithm.\n",
"http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.42.3443&rep=rep1&type=pdf\n",
"\n",
"Knowing the input position for a X-ray on the detector and its propagation vector, this algorithm allows us to calculate the length of the path in all voxel it crosses in a fairly efficient way.\n",
"\n",
"To speed up the calculation, we will use a few tricks:\n",
"* One ray never crosses more than 16 pixels, which is reasonable considering the incidance angle \n",
"* we use numba to speed-up the calculation of loops in python\n",
"* We will allocate the needed memory by chuncks of 1 million elements\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"from numba import jit \n",
"\n",
"BLOCK_SIZE = 1<<20 # 1 milion\n",
"BUFFER_SIZE = 16 \n",
"BIG = numpy.finfo(numpy.float32).max\n",
"\n",
"mask = numpy.load(\"mask.npy\").astype(numpy.int8)\n",
"from scipy.sparse import csr_matrix, csc_matrix, linalg"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(array([0, 0, 1, 1], dtype=int32), array([0, 1, 1, 2], dtype=int32), array([0.00029791, 0.00029791, 0.00059583, 0.00059583], dtype=float32))\n",
"The slowest run took 9.46 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
"10.1 µs ± 8.53 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
"10.3 µs ± 9.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
]
}
],
"source": [
"@jit\n",
"def calc_one_ray(entx, enty, \n",
" kx, ky, kz,\n",
" vox, voy, voz):\n",
" \"\"\"For a ray, entering at position (entx, enty), with a propagation vector (kx, ky,kz),\n",
" calculate the length spent in every voxel where energy is deposited from a bunch of photons comming in the detector \n",
" at a given position and and how much energy they deposit in each voxel. \n",
" \n",
" Direct implementation of http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.42.3443&rep=rep1&type=pdf\n",
" \n",
" :param entx, enty: coordinate of the entry point in meter (2 components, x,y)\n",
" :param kx, ky, kz: vector with the direction of the photon (3 components, x,y,z)\n",
" :param vox, voy, voz: size of the voxel in meter (3 components, x,y,z)\n",
" :return: coordinates voxels in x, y and length crossed when leaving the associated voxel\n",
" \"\"\"\n",
" array_x = numpy.empty(BUFFER_SIZE, dtype=numpy.int32)\n",
" array_x[:] = -1\n",
" array_y = numpy.empty(BUFFER_SIZE, dtype=numpy.int32)\n",
" array_y[:] = -1\n",
" array_len = numpy.empty(BUFFER_SIZE, dtype=numpy.float32)\n",
" \n",
" #normalize the input propagation vector\n",
" n = numpy.sqrt(kx*kx + ky*ky + kz*kz)\n",
" kx /= n\n",
" ky /= n\n",
" kz /= n\n",
" \n",
" assert kz>0\n",
" step_X = -1 if kx<0.0 else 1\n",
" step_Y = -1 if ky<0.0 else 1\n",
" \n",
" assert vox>0\n",
" assert voy>0\n",
" assert voz>0\n",
" \n",
" X = int(entx//vox)\n",
" Y = int(enty//voy)\n",
" \n",
" if kx>0.0:\n",
" t_max_x = ((entx//vox+1)*(vox)-entx)/ kx\n",
" elif kx<0.0:\n",
" t_max_x = ((entx//vox)*(vox)-entx)/ kx\n",
" else:\n",
" t_max_x = BIG\n",
"\n",
" if ky>0.0:\n",
" t_max_y = ((enty//voy+1)*(voy)-enty)/ ky\n",
" elif ky<0.0:\n",
" t_max_y = ((enty//voy)*(voy)-enty)/ ky\n",
" else:\n",
" t_max_y = BIG\n",
" \n",
" #Only one case for z as the ray is travelling in one direction only\n",
" t_max_z = voz / kz\n",
" \n",
" t_delta_x = abs(vox/kx) if kx!=0 else BIG\n",
" t_delta_y = abs(voy/ky) if ky!=0 else BIG\n",
" t_delta_z = voz/kz\n",
" \n",
" finished = False\n",
" last_id = 0\n",
" array_x[last_id] = X\n",
" array_y[last_id] = Y\n",
" \n",
" while not finished:\n",
" if t_max_x < t_max_y:\n",
" if t_max_x < t_max_z:\n",
" array_len[last_id] = t_max_x\n",
" last_id+=1\n",
" X += step_X\n",
" array_x[last_id] = X\n",
" array_y[last_id] = Y\n",
" t_max_x += t_delta_x\n",
" else:\n",
" array_len[last_id] = t_max_z\n",
" finished = True\n",
" else:\n",
" if t_max_y < t_max_z:\n",
" array_len[last_id] = t_max_y\n",
" last_id+=1\n",
" Y += step_Y\n",
" array_x[last_id] = X\n",
" array_y[last_id] = Y \n",
" t_max_y += t_delta_y\n",
" else:\n",
" array_len[last_id] = t_max_z\n",
" finished = True\n",
" if last_id>=array_len.size-1:\n",
" print(\"resize arrays\")\n",
" old_size = len(array_len)\n",
" new_size = (old_size//BUFFER_SIZE+1)*BUFFER_SIZE\n",
" new_array_x = numpy.empty(new_size, dtype=numpy.int32)\n",
" new_array_x[:] = -1\n",
" new_array_y = numpy.empty(new_size, dtype=numpy.int32)\n",
" new_array_y[:] = -1\n",
" new_array_len = numpy.empty(new_size, dtype=numpy.float32)\n",
" new_array_x[:old_size] = array_x\n",
" new_array_y[:old_size] = array_y\n",
" new_array_len[:old_size] = array_len\n",
" array_x = new_array_x\n",
" array_y = new_array_y\n",
" array_len = new_array_len\n",
" return array_x[:last_id], array_y[:last_id], array_len[:last_id]\n",
"\n",
"print(calc_one_ray(0.0,0.0, 1,1,1, 172e-6, 172e-6, 450e-6))\n",
"import random\n",
"%timeit calc_one_ray(10+random.random(),11+random.random(),\\\n",
" random.random()-0.5,random.random()-0.5,0.5+random.random(), \\\n",
" vox, voy, voz)\n",
"%timeit calc_one_ray.py_func(10+random.random(),11+random.random(),\\\n",
" random.random()-0.5,random.random()-0.5,0.5+random.random(), \\\n",
" vox, voy, voz)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we are able to perform raytracing for any ray comming in the detector, we can calculate the contribution to the neighboring pixels, using the absorption law (the length travelled is already known). \n",
"To average-out the signal, we will sample a few dozens of rays per pixel to get an approximatation of the volumic integrale. \n",
"\n",
"Now we need to store the results so that this transformation can be represented as a sparse matrix multiplication:\n",
"\n",
"b = M.a\n",
"\n",
"Where b is the recorded image (blurred) and a is the \"perfect\" signal. \n",
"M being the sparse matrix where every pixel of a gives a limited number of contribution to b.\n",
"\n",
"Each pixel in *b* is represented by one line in *M* and we store the indices of *a* of interest with the coefficients of the matrix.\n",
"So if a pixel i,j contributes to (i,j), (i+1,j), (i+1,j+1), there are only 3 elements in the line. \n",
"This is advantagous for storage.\n",
"\n",
"We will use the CSR sparse matrix representation:\n",
"https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_.28CSR.2C_CRS_or_Yale_format.29\n",
"where there are 3 arrays:\n",
"* data: containing the actual non zero values\n",
"* indices: for a given line, it contains the column number of the assocated data (at the same indice)\n",
"* idptr: this array contains the index of the start of every line.\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"from numba import jitclass, int8, int32, int64, float32, float64\n",
"spec = [(\"vox\",float64),(\"voy\",float64),(\"voz\",float64),(\"mu\",float64),\n",
" (\"dist\",float64),(\"poni1\",float64),(\"poni2\",float64),\n",
" (\"width\", int64),(\"height\", int64),(\"mask\", int8[:,:]),\n",
" (\"sampled\", int64), (\"data\", float32[:]),(\"indices\", int32[:]),(\"idptr\", int32[:]),\n",
" ]\n",
"@jitclass(spec)\n",
"class ThickDetector(object):\n",
" \"Calculate the point spread function as function of the geometry of the experiment\"\n",
" \n",
" def __init__(self, vox, voy, thickness, mask, mu, \n",
" dist, poni1, poni2):\n",
" \"\"\"Constructor of the class:\n",
" \n",
" :param vox, voy: detector pixel size in the plane\n",
" :param thickness: thickness of the sensor in meters\n",
" :param mask: \n",
" :param mu: absorption coefficient of the sensor material\n",
" :param dist: sample detector distance as defined in the geometry-file\n",
" :param poni1, poni2: coordinates of the PONI as defined in the geometry \n",
" \"\"\"\n",
" self.vox = vox\n",
" self.voy = voy\n",
" self.voz = thickness\n",
" self.mu = mu\n",
" self.dist=dist\n",
" self.poni1 = poni1\n",
" self.poni2 = poni2\n",
" self.width = mask.shape[-1]\n",
" self.height = mask.shape[0]\n",
" self.mask = mask\n",
" self.sampled = 0\n",
" self.data = numpy.zeros(BLOCK_SIZE, dtype=numpy.float32)\n",
" self.indices = numpy.zeros(BLOCK_SIZE,dtype=numpy.int32)\n",
" self.idptr = numpy.zeros(self.width*self.height+1, dtype=numpy.int32)\n",
" \n",
" def calc_one_ray(self, entx, enty):\n",
" \"\"\"For a ray, entering at position (entx, enty), with a propagation vector (kx, ky,kz),\n",
" calculate the length spent in every voxel where energy is deposited from a bunch of photons comming in the detector \n",
" at a given position and and how much energy they deposit in each voxel. \n",
"\n",
" Direct implementation of http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.42.3443&rep=rep1&type=pdf\n",
"\n",
" :param entx, enty: coordinate of the entry point in meter (2 components, x,y)\n",
" :return: coordinates voxels in x, y and length crossed when leaving the associated voxel\n",
" \"\"\"\n",
" array_x = numpy.empty(BUFFER_SIZE, dtype=numpy.int32)\n",
" array_x[:] = -1\n",
" array_y = numpy.empty(BUFFER_SIZE, dtype=numpy.int32)\n",
" array_y[:] = -1\n",
" array_len = numpy.empty(BUFFER_SIZE, dtype=numpy.float32)\n",
"\n",
" #normalize the input propagation vector\n",
" kx = entx - self.poni2\n",
" ky = enty - self.poni1\n",
" kz = self.dist\n",
" n = numpy.sqrt(kx*kx + ky*ky + kz*kz)\n",
" kx /= n\n",
" ky /= n\n",
" kz /= n\n",
"\n",
" step_X = -1 if kx<0.0 else 1\n",
" step_Y = -1 if ky<0.0 else 1\n",
"\n",
" X = int(entx/self.vox)\n",
" Y = int(enty/self.voy)\n",
"\n",
" if kx>0.0:\n",
" t_max_x = ((entx//self.vox+1)*(self.vox)-entx)/ kx\n",
" elif kx<0.0:\n",
" t_max_x = ((entx//self.vox)*(self.vox)-entx)/ kx\n",
" else:\n",
" t_max_x = BIG\n",
"\n",
" if ky>0.0:\n",
" t_max_y = ((enty//self.voy+1)*(self.voy)-enty)/ ky\n",
" elif ky<0.0:\n",
" t_max_y = ((enty//self.voy)*(self.voy)-enty)/ ky\n",
" else:\n",
" t_max_y = BIG\n",
"\n",
" #Only one case for z as the ray is travelling in one direction only\n",
" t_max_z = self.voz / kz\n",
"\n",
" t_delta_x = abs(self.vox/kx) if kx!=0 else BIG\n",
" t_delta_y = abs(self.voy/ky) if ky!=0 else BIG\n",
" t_delta_z = self.voz/kz\n",
"\n",
" finished = False\n",
" last_id = 0\n",
" array_x[last_id] = X\n",
" array_y[last_id] = Y\n",
"\n",
" while not finished:\n",
" if t_max_x < t_max_y:\n",
" if t_max_x < t_max_z:\n",
" array_len[last_id] = t_max_x\n",
" last_id+=1\n",
" X += step_X\n",
" array_x[last_id] = X\n",
" array_y[last_id] = Y\n",
" t_max_x += t_delta_x\n",
" else:\n",
" array_len[last_id] = t_max_z\n",
" last_id+=1\n",
" finished = True\n",
" else:\n",
" if t_max_y < t_max_z:\n",
" array_len[last_id] = t_max_y\n",
" last_id+=1\n",
" Y += step_Y\n",
" array_x[last_id] = X\n",
" array_y[last_id] = Y \n",
" t_max_y += t_delta_y\n",
" else:\n",
" array_len[last_id] = t_max_z\n",
" last_id+=1\n",
" finished = True\n",
" if last_id>=array_len.size-1:\n",
" print(\"resize arrays\")\n",
" old_size = len(array_len)\n",
" new_size = (old_size//BUFFER_SIZE+1)*BUFFER_SIZE\n",
" new_array_x = numpy.empty(new_size, dtype=numpy.int32)\n",
" new_array_x[:] = -1\n",
" new_array_y = numpy.empty(new_size, dtype=numpy.int32)\n",
" new_array_y[:] = -1\n",
" new_array_len = numpy.empty(new_size, dtype=numpy.float32)\n",
" new_array_x[:old_size] = array_x\n",
" new_array_y[:old_size] = array_y\n",
" new_array_len[:old_size] = array_len\n",
" array_x = new_array_x\n",
" array_y = new_array_y\n",
" array_len = new_array_len\n",
" return array_x[:last_id], array_y[:last_id], array_len[:last_id]\n",
"\n",
" def one_pixel(self, row, col, sample):\n",
" \"\"\"calculate the contribution of one pixel to the sparse matrix and populate it.\n",
"\n",
" :param row: row index of the pixel of interest\n",
" :param col: column index of the pixel of interest\n",
" :param sample: Oversampling rate, 10 will thow 10x10 ray per pixel\n",
"\n",
" :return: the extra number of pixel allocated\n",
" \"\"\"\n",
" if self.mask[row, col]:\n",
" return (numpy.empty(0, dtype=numpy.int32),\n",
" numpy.empty(0, dtype=numpy.float32))\n",
"\n",
" counter = 0\n",
" tmp_size = 0\n",
" last_buffer_size = BUFFER_SIZE\n",
" tmp_idx = numpy.empty(last_buffer_size, dtype=numpy.int32)\n",
" tmp_idx[:] = -1\n",
" tmp_coef = numpy.zeros(last_buffer_size, dtype=numpy.float32)\n",
"\n",
" pos = row * self.width + col\n",
" start = self.idptr[pos]\n",
" for i in range(sample):\n",
" posx = (col+1.0*i/sample)*vox\n",
" for j in range(sample):\n",
" posy = (row+1.0*j/sample)*voy\n",
" array_x, array_y, array_len = self.calc_one_ray(posx, posy)\n",
"\n",
" rem = 1.0\n",
" for i in range(array_x.size):\n",
" x = array_x[i]\n",
" y = array_y[i]\n",
" l = array_len[i]\n",
" if (x<0) or (y<0) or (y>=self.height) or (x>=self.width):\n",
" break\n",
" elif (self.mask[y, x]):\n",
" continue\n",
" idx = x + y*self.width\n",
" dos = numpy.exp(-self.mu*l)\n",
" value = rem - dos\n",
" rem = dos\n",
" for j in range(last_buffer_size):\n",
" if tmp_size >= last_buffer_size:\n",
" #Increase buffer size\n",
" new_buffer_size = last_buffer_size + BUFFER_SIZE\n",
" new_idx = numpy.empty(new_buffer_size, dtype=numpy.int32)\n",
" new_coef = numpy.zeros(new_buffer_size, dtype=numpy.float32)\n",
" new_idx[:last_buffer_size] = tmp_idx\n",
" new_idx[last_buffer_size:] = -1\n",
" new_coef[:last_buffer_size] = tmp_coef\n",
" last_buffer_size = new_buffer_size\n",
" tmp_idx = new_idx\n",
" tmp_coef = new_coef\n",
"\n",
" if tmp_idx[j] == idx:\n",
" tmp_coef[j] += value\n",
" break\n",
" elif tmp_idx[j] < 0:\n",
" tmp_idx[j] = idx\n",
" tmp_coef[j] = value\n",
" tmp_size +=1\n",
" break \n",
" return tmp_idx[:tmp_size], tmp_coef[:tmp_size]\n",
"\n",
" def calc_csr(self, sample):\n",
" \"\"\"Calculate the CSR matrix for the whole image\n",
" :param sample: Oversampling factor\n",
" :return: CSR matrix\n",
" \"\"\"\n",
" size = self.width * self.height\n",
" allocated_size = BLOCK_SIZE\n",
" idptr = numpy.zeros(size+1, dtype=numpy.int32) \n",
" indices = numpy.zeros(allocated_size, dtype=numpy.int32)\n",
" data = numpy.zeros(allocated_size, dtype=numpy.float32)\n",
" self.sampled = sample*sample\n",
" pos = 0\n",
" start = 0\n",
" for row in range(self.height):\n",
" for col in range(self.width): \n",
" line_idx, line_coef = self.one_pixel(row, col, sample)\n",
" line_size = line_idx.size\n",
" if line_size == 0:\n",
" new_size = 0\n",
" pos+=1\n",
" idptr[pos] = start\n",
" continue\n",
"\n",
" stop = start + line_size\n",
" \n",
" if stop >= allocated_size:\n",
" new_buffer_size = allocated_size + BLOCK_SIZE\n",
" new_idx = numpy.zeros(new_buffer_size, dtype=numpy.int32)\n",
" new_coef = numpy.zeros(new_buffer_size, dtype=numpy.float32)\n",
" new_idx[:allocated_size] = indices\n",
" new_coef[:allocated_size] = data\n",
" allocated_size = new_buffer_size\n",
" indices = new_idx\n",
" data = new_coef\n",
"\n",
" indices[start:stop] = line_idx\n",
" data[start:stop] = line_coef\n",
" pos+=1\n",
" idptr[pos] = stop\n",
" start = stop\n",
" \n",
" last = idptr[-1]\n",
" self.data = data\n",
" self.indices = indices\n",
" self.idptr = idptr\n",
" return (self.data[:last]/self.sampled, indices[:last], idptr)\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 2.04 s, sys: 24 ms, total: 2.06 s\n",
"Wall time: 2.06 s\n"
]
},
{
"data": {
"text/plain": [
"(array([0., 0., 0., ..., 0., 0., 0.], dtype=float32),\n",
" array([ 2, 2, 4, ..., 1023180, 1023181, 1023182],\n",
" dtype=int32),\n",
" array([ 0, 0, 0, ..., 1902581, 1902582, 1902583],\n",
" dtype=int32))"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"thick = ThickDetector(vox,voy, thickness=thickness, mu=mu, dist=dist, poni1=poni1, poni2=poni2, mask=mask)\n",
"%time thick.calc_csr(1)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 20.6 s, sys: 152 ms, total: 20.8 s\n",
"Wall time: 20.6 s\n"
]
}
],
"source": [
"thick = ThickDetector(vox,voy, thickness=thickness, mu=mu, dist=dist, poni1=poni1, poni2=poni2, mask=mask)\n",
"%time pre_csr = thick.calc_csr(8)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Validation of the CSR matrix obtained:\n",
"\n",
"For this we will build a simple 2D image with one pixel in a regular grid and calculate the effect of the transformation calculated previously on it. "
]
},
{
"cell_type": "code",
"execution_count": 11,
"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,iVBORw0KGgoAAAANSUhEUgAAA2gAAANoCAYAAAC1Fsk9AAAgAElEQVR4nOzdfbCfdX3n/yshRAkpKCogGAXFqEiwtgTbFczMTyvbRctqf/FXqC7dutyoW9d7WxVPraKj4l3qTbWi3R1BtumUu9URxcRaAZ2srR1LcYqCaKwaVsNZjMJg8v79AeeMx5ycfEKuc30/7+vzeM48Z+DLOddJXnN8vfM65iRdAAAAAACqoJv0DwAAAAAAcC8GGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQZaD9x8881x7rnnxpo1a2Lp0qWxbt26ed/uO9/5Tpx11lnxkIc8JB74wAfGCSecEJdffvmwP1gAwBx0OACgJgy0Hrjiiiti1apVsX79+li9evW8x33r1q1x9NFHxzOf+cy44oor4tprr413v/vdcemllw7/AwYAzKLDAQA1YaD1wM6dO2f/+Ywzzpj3uJ955plxyimnzHlbAMDk0eEAgJow0HpmvuN+xx13xIEHHugrrQBQOTocADBpDLSeme+4b9q0Kbqui8suuyxOOeWUWLZsWRx55JFxwQUXxM9//vPJ/EABALuhwwEAk8ZA65n5jvull14aXdfFIYccEq9+9atj06ZN8eY3vzmWLVsWF1544YLPm56ejq1bt856yy23xBe+8IW47bbb5rxOktm87bbbYsuWLXH33XcvYivvGzqcJPdujf09Jgy0npnvuF9yySXRdV2sX79+zut/9Ed/FIccckjs2rVrj8+bmpqKrutIcrRu2bJlMer4fqHDSbLcmvp7TBhoPTPfcf/0pz8dXdfFhz70oTmvX3755dF1XXznO9/Z4/N++auvX/7yl6PruljbPT1O6U4nyWrd21dgt2zZEl3XxW233bYYdXy/GKrDTz7+nDj1Sa/Yo08+44JenPTnwGK6UH774pN/54JenHQei5bzr76yF+U8XNa/9uwLejFbf48JA61n5jvut95667zH/W//9m+j67rYunVr8fO3bt0aXdfFKd3p8Ywl/y9JVmtpn+1LBy42Q3X4qU96RfzW2j/do08586JenPTnwGK6UH774lN+76JenHQei5bzyW/qRTkPl/Vv/H8X9WJJl9XU32PCQOuZPf0RzU984hPjd3/3d+e89pKXvCQe+tCH7tPzDTSSWSzts5oO/FAdbqDtvwbaQDkbaOmyNtDyY6D1wI4dO2Ljxo2xcePGWLt2bRx//PGz/75t27aIuPcrrUuWLImXv/zl8bnPfS6mpqbigAMOiPe///379LEMNJJZLO2zSR/4SXS4gbb/GmgD5WygpcvaQMuPgdYDM7/9ZT43b948+3aXXHJJPOEJT4jly5fHscceG+9+97v3+WMZaCSzWNpnkz7wk+hwA23/NdAGytlAS5e1gZYfAy0ZBhrJLJb2WUsH3kDrTwNtoJwNtHRZG2j5MdCSYaCRzGJpn7V04A20/jTQBsrZQEuXtYGWHwMtGaUDbef3j+vFSZdVBmUt57HZV9alfdbSgZ/5Oa879kVx2urX7NFHbbioFyf9ubSYnva41/airOVci31lfcz7LurFki5rqb+HxEBLhoFWn7KW89g00BYPA60/DQc5j00DDTMYaMkw0OpT1nIemwba4mGg9afhIOexaaBhBgMtGQZafcpazmPTQFs8DLT+NBzkPDYNNMxgoCXDQKtPWct5bBpoi4eB1p+Gg5zHpoGGGQy0ZBho9SlrOY9NA23xMND603CQ89g00DCDgZYMA60+ZS3nsWmgLR4GWn8aDnIemwYaZjDQkmGg1aes5Tw2DbTFw0DrT8NBzmPTQMMMBloyDLT6lLWcx6aBtngYaP1pOMh5bBpomMFAS4aBVp+ylvPYNNAWDwOtPw0HOY9NAw0zGGjJMNDqU9ZyHpsG2uJhoPWn4SDnsWmgYQYDLRkGWn3KWs5j00BbPAy0/jQc5Dw2DTTMYKAlw0CrT1nLeWwaaIuHgdafhoOcx6aBhhkMtGQYaPUpazmPTQNt8TDQ+tNwkPPYNNAwg4GWDAOtPmUt57FpoC0eBlp/Gg5yHpsGGmYw0JJhoNWnrOU8Ng20xcNA60/DQc5j00DDDAZaMgy0+pS1nMemgbZ4GGj9aTjIeWwaaJjBQEuGgVafspbz2DTQFg8DrT8NBzmPTQMNMxhoySgdaCQ5aUv7rKUDr8NJZrGky1rq7yEx0JLhuJPMYmmftXTgdTjJLJZ0WUv9PSQGWjIcd5JZLO2zlg68DieZxZIua6m/h8RAS4bjTjKLpX3W0oHX4SSzWNJlLfX3kBhoPXDzzTfHueeeG2vWrImlS5fGunXrFnz7yy+/PLquiyc96Un7/LEcd5JZLO2zSR94HU6Su1vSZZPu77FioPXAFVdcEatWrYr169fH6tWrFzzuP/3pT+OYY46JI444wnEnOWpL+2zSB16Hk+TulnTZpPt7rBhoPbBz587Zfz7jjDMWPO4XXHBBPO1pT4uzzz7bcSc5akv7bNIHXoeT5O6WdNmk+3usGGg9s9Bx/+Y3vxkrVqyIr33ta447ydFb2mc1HXgdTpL3WtJlNfX3mDDQemah43766afH+eefHxHhuJMcvaV9VtOB1+Ekea8lXVZTf48JA61n9nTcr7rqqnjwgx8ct99+e0SUH/fp6enYunXrrFu2bHHcSaZwb9R44HU4Sd7rQtTY32PCQOuZ+Y77z372s3j0ox8dGzZsmH2t9LhPTU1F13W76biTrN29UeOB1+Ekea8LUWN/jwkDrWfmO+5ve9vb4rjjjovbb789tm/fHtu3b48zzzwzTjjhhNi+fXvcfffde3yer76SzOreqPHA63CSvNeFqLG/x4SB1jPzHfezzz573q+gzvjxj3+8+Pm+f4FkFkv7rKYDr8NJ8l5Luqym/h4TBlrPzHfcb7rppti8efMcTzvttHjMYx4Tmzdvju9///vFz3fcSWaxtM9qOvA6nCTvtaTLaurvMWGg9cCOHTti48aNsXHjxli7dm0cf/zxs/++bdu2ed/HnwBGcuyW9tmkD7wOJ8ndLemySff3WDHQeuDWW2/d42992bx587zv47iTHLulfTbpA6/DSXJ3S7ps0v09Vgy0ZDjuJLNY2mctHXgdTjKLJV3WUn8PiYGWDMedZBZL+6ylA6/DSWaxpMta6u8hMdCS4biTzGJpn7V04HU4ySyWdFlL/T0kBloyHHeSWSzts5YOvA4nmcWSLmupv4fEQEuG404yi6V91tKB1+Eks1jSZS3195AYaMlw3ElmsbTPWjrwOpxkFku6rKX+HhIDLRmOO8kslvZZSwdeh5PMYkmXtdTfQ2KgJcNxJ5nF0j5r6cDrcJJZLOmylvp7SAy0ZDjuJLNY2mctHXgdTjKLJV3WUn8PiYGWDMedZBZL+6ylA6/DSWaxpMta6u8hMdCS4biTzGJpn7V04HU4ySyWdFlL/T0kBloyHHeSWSzts5YOvA4nmcWSLmupv4fEQEuG404yi6V91tKB1+Eks1jSZS3195AYaMlw3ElmsbTPWjrwOpxkFku6rKX+HhIDLRmOO8kslvZZSwdeh5PMYkmXtdTfQ2KgJcNxJ5nF0j5r6cDrcJJZLOmylvp7SAy0ZDjuJLNY2mctHXgdTjKLJV3WUn8PiYGWDMedZBZL+6ylA6/DSWaxpMta6u8hMdCS4biTzGJpn7V04HU4ySyWdFlL/T0kBloyHHeSWSzts5YOvA4nmcWSLmupv4fEQEuG404yi6V91tKB1+Eks1jSZS3195AYaMlw3ElmsbTPWjrwOpxkFku6rKX+HhIDLRmOO8kslvZZSwdeh5PMYkmXtdTfQ2KgJcNxJ5nF0j5r6cDrcJJZLOmylvp7SAy0Hrj55pvj3HPPjTVr1sTSpUtj3bp1c/779PR0TE1Nxdq1a+PQQw+Nww8/PJ797GfHP//zP+/zx3LcSWaxtM8mfeB1OEnubkmXTbq/x4qB1gNXXHFFrFq1KtavXx+rV6/e7bh//etfjyOPPDJe//rXxzXXXBNXXnllnHrqqXHwwQfHN77xjX36WI47ySyW9tmkD7wOJ8ndLemySff3WDHQemDnzp2z/3zGGWfsdtx/8pOfxI4dO+a8duedd8Zhhx0WL3vZy/bpYznuJLNY2meTPvA6nCR3t6TLJt3fY8VA65n5jvueOPnkk+N5z3vePj3fcSeZxdI+q+nA63CSvNeSLqupv8eEgdYzpcd9+/btsWLFipiamtqn5zvuJLNY2mc1HXgdTpL3WtJlNfX3mDDQeqb0uJ9zzjmxcuXK+N73vrfg201PT8fWrVtn3bJli+NOMoV7o8YDr8NJ8l4Xosb+HhMGWs+UHPePfexj0XVd/NVf/dVenzc1NRVd1+2m406ydvdGjQdeh5PkvS5Ejf09Jgy0ntnbcf/0pz8dy5YtiwsuuKDoeb76SjKre6PGA6/DSfJeF6LG/h4TBlrPLHTcb7jhhlixYkW88IUvvN/P9/0LJLNY2mc1HXgdTpL3WtJlNfX3mDDQemZPx/3GG2+Mww47LJ71rGfFPffcc7+f77iTzGJpn9V04HU4Sd5rSZfV1N9jwkDrgR07dsTGjRtj48aNsXbt2jj++ONn/33btm3xwx/+MB7xiEfE0UcfHZ///OfjhhtumPXGG2/cp4/luJPMYmmfTfrA63CS3N2SLpt0f48VA60Hbr311nm/Cbzruti8eXNs3rx5j/+99O/bmcFxJ5nF0j6b9IHX4SS5uyVdNun+HisGWjIcd5JZLO2zlg68DieZxZIua6m/h8RAS4bjTjKLpX3W0oHX4SSzWNJlLfX3kBhoyXDcSWaxtM9aOvA6nGQWS7qspf4eEgMtGY47ySyW9llLB37m57xu9R/FaU983R79f059Sy9O+nNgMV0ov31R1nvJ+YTX9+LTT3lLL046j5ayLumylvp7SAy0ZBhoJLNY2mctHXgDrT8NtIFyrmw0TDqPlrIu6bKW+ntIDLRkGGgks1jaZy0deAOtPw20gXKubDRMOo+Wsi7pspb6e0gMtGQYaCSzWNpnLR14A60/DbSBcq5sNEw6j5ayLumylvp7SAy0ZBhoJLNY2mctHXgDrT8NtIFyrmw0TDqPlrIu6bKW+ntIDLRkGGgks1jaZy0deAOtPw20gXKubDRMOo+Wsi7pspb6e0gMtGQYaCSzWNpnLR14A60/DbSBcq5sNEw6j5ayLumylvp7SAy0ZBhoJLNY2mctHXgDrT8NtIFyrmw0TDqPlrIu6bKW+ntIDLRkGGgks1jaZy0deAOtPw20gXKubDRMOo+Wsi7pspb6e0gMtGQYaCSzWNpnLR14A60/DbSBcq5sNEw6j5ayLumylvp7SAy0ZBhoJLNY2mctHXgDrT8NtIFyrmw0TDqPlrIu6bKW+ntIDLRkGGgks1jaZy0deAOtPw20gXKubDRMOo+Wsi7pspb6e0gMtGQYaCSzWNpnLR14A60/DbSBcq5sNEw6j5ayLumylvp7SAy0ZBhoJLNY2mctHXgDrT8NtIFyrmw0TDqPlrIu6bKW+ntIDLRkGGgks1jaZy0deAOtPw20gXKubDRMOo+Wsi7pspb6e0gMtGQYaCSzWNpnLR14A60/DbSBcq5sNEw6j5ayLumylvp7SAy0ZBhoJLNY2mctHXgDrT8NtIFyrmw0TDqPlrIu6bKW+ntIDLRkGGgks1jaZy0deAOtPw20gXKubDRMOo+Wsi7pspb6e0gMtGSUDrSd3z+uFyddVhmUtZzHZl9Zl/ZZSwd+dqAd+6I4bfVr9uijNlzUi5P+XFpMF8pvXzzmfRf14qTzWLScH/faXpRzvqxLuqyl/h4SAy0ZBlp9ylrOY9NAWzwMtP400AbKubLRMOk8Wsq6pMta6u8hMdCSYaDVp6zlPDYNtMXDQOtPA22gnCsbDZPOo6WsS7qspf4eEgMtGQZafcpazmPTQFs8DLT+NNAGyrmy0TDpPFrKuqTLWurvITHQeuDmm2+Oc889N9asWRNLly6NdevWzft2H/3oR+Oxj31sPOABD4gTTzwxrr766n3+WAZafcpazmOztYE2iQ430PZfA22gnCsbDZPOo6WsS7ps0v09Vgy0Hrjiiiti1apVsX79+li9evW8x/2Tn/xkLFmyJN7whjfEpk2b4rzzzotly5bFDTfcsE8fy0CrT1nLeWy2NtAm0eEG2v5roA2Uc2WjYdJ5tJR1SZdNur/HioHWAzt37pz95zPOOGPe47569eo466yz5rz2m7/5m/Hbv/3b+/SxDLT6lLWcx2ZrA20SHW6g7b8G2kA5VzYaJp1HS1mXdNmk+3usGGg9M99x/9a3vhVd18WVV1455/X3ve99sXz58rjrrruKn2+g1aes5Tw2Wxtov8hQHW6g7b8G2kA5VzYaJp1HS1mXdFlN/T0mDLSeme+4f+pTn4qu6+Lmm2+e8/pnP/vZ6LoubrrppuLnG2j1KWs5j00Dbd2c1xajww20/ddAGyjnykbDpPNoKeuSLqupv8eEgdYz8x33T3ziE9F1Xdx+++1zXt+yZUt0XRfXXXfdHp83PT0dW7dunXXmfQy0epS1nMemgbZuzmuL0eEG2v5roA2Uc2WjYdJ5tJT1QtTY32PCQOuZvo/71NRUdF23mwZaPcpazmPTQFs357XF6HADbf810AbKubLRMOk8Wsp6IWrs7zFhoPVM3789xv+DVr+ylvPYNNDWzXltMTrcQNt/DbSBcq5sNEw6j5ayXoga+3tMGGg9s9A3mF911VVzXt+wYUMsX7487r777uLn+x60+pS1nMemgbZuzmuL0eEG2v5roA2Uc2WjYdJ5tJR1SZfV1N9jwkDrmYX+iObnP//5c1576lOf6o/ZH4GylvPYNNDW7fZ63x1uoO2/BtpAOVc2GiadR0tZl3RZTf09Jgy0HtixY0ds3LgxNm7cGGvXro3jjz9+9t+3bdsWERGXXnppLFmyJN74xjfG5s2b4/zzz49ly5bF9ddfv08fy0CrT1nLeWy2NtAm0eEG2v5roA2Uc2WjYdJ5tJR1SZdNur/HioHWA7feeuu83wTedV1s3rx59u0++tGPxnHHHRfLly+PNWvWxNVXX73PH8tAq09Zy3lstjbQJtHhBtr+a6ANlHNlo2HSebSUdUmXTbq/x4qBlgwDrT5lLeex2dpAGxIDrT8NtIFyrmw0TDqPlrIu6bKW+ntIDLRkGGj1KWs5j00DbfEw0PrTQBso58pGw6TzaCnrki5rqb+HxEBLRulAI8lJW9pnLR14HU4yiyVd1lJ/D4mBlgzHnWQWS/uspQOvw0lmsaTLWurvITHQkuG4k8xiaZ+1dOB1OMkslnRZS/09JAZaMhx3klks7bOWDrwOJ5nFki5rqb+HxEBLhuNOMoulfdbSgdfhJLNY0mUt9feQGGjJcNxJZrG0z1o68DqcZBZLuqyl/h4SAy0ZjjvJLJb2WUsHXoeTzGJJl7XU30NioCXDcSeZxdI+a+nA63CSWSzpspb6e0gMtGQ47iSzWNpnLR14HU4yiyVd1lJ/D4mBlgzHnWQWS/uspQOvw0lmsaTLWurvITHQkuG4k8xiaZ+1dOB1OMkslnRZS/09JAZaMhx3klks7bOWDrwOJ5nFki5rqb+HxEBLhuNOMoulfdbSgdfhJLNY0mUt9feQGGjJcNxJZrG0z1o68DqcZBZLuqyl/h4SAy0ZjjvJLJb2WUsHXoeTzGJJl7XU30NioCXDcSeZxdI+a+nA63CSWSzpspb6e0gMtGQ47iSzWNpnLR14HU4yiyVd1lJ/D4mBlgzHnWQWS/uspQOvw0lmsaTLWurvITHQkuG4k8xiaZ+1dOB1OMkslnRZS/09JAZaMhx3klks7bOWDrwOJ5nFki5rqb+HxEBLhuNOMoulfdbSgdfhJLNY0mUt9feQGGjJcNxJZrG0z1o68DqcZBZLuqyl/h4SA21Arrzyyjj55JNj5cqVceSRR8bznve8uOWWW/bpGY47ySyW9lmWA6/DSbZkSZdl6e9sGGgDsXnz5li6dGn8wR/8QXzuc5+Lyy67LFavXh2rV6+On/3sZ8XPcdxJZrG0zzIceB1OsjVLuixDf2fEQBuI8847L4499tjYtWvX7GubNm2Kruvi+uuvL36O404yi6V9luHA63CSrVnSZRn6OyMG2kD84R/+YZx44olzXvvqV78aXdfFddddV/wcx51kFkv7LMOB1+EkW7OkyzL0d0YMtIH44he/GMuWLYsPfOADcccdd8S3vvWtOO200+Kkk06KnTt3Fj/HcSeZxdI+y3DgdTjJ1izpsgz9nREDbUCuuuqqWLlyZXRdF13XxZOf/OT44Q9/uOD7TE9Px9atW2fdsmWL404yhXsj24HX4SRbciGy9Xc2DLSB+NKXvhQPetCD4hWveEVs2rQpNm7cGCeeeGKsXbs27rrrrj2+39TU1OwvBn5Rx51k7e6NTAdeh5NszYXI1N8ZMdAG4td//ddj/fr1c1777ne/G0uWLImLL754j+/nq68ks7o3Mh14HU6yNRciU39nxEAbiIMOOiguvPDC3V5/2MMeFq973euKn+P7F0hmsbTPMhx4HU6yNUu6LEN/Z8RAG4jHP/7x8dznPnfOa9/+9rdjyZIl8ZGPfKT4OY47ySyW9lmGA6/DSbZmSZdl6O+MGGgD8d73vje6rouXvvSls3/J6QknnBAPf/jD48c//nHxcxx3klks7bMMB16Hk2zNki7L0N8ZMdAGYteuXfGhD30oTjzxxDj44IPjyCOPjOc85znxjW98Y5+e47iTzGJpn2U48DqcZGuWdFmG/s6IgZYMx51kFkv7rKUDr8NJZrGky1rq7yEx0JLhuJPMYmmftXTgdTjJLJZ0WUv9PSQGWjIcd5JZLO2zlg68DieZxZIua6m/h8RAS4bjTjKLpX3W0oHX4SSzWNJlLfX3kBhoyXDcSWaxtM9aOvA6nGQWS7qspf4eEgMtGY47ySyW9llLB16Hk8xiSZe11N9DYqAlw3EnmcXSPmvpwOtwklks6bKW+ntIDLRkOO4ks1jaZy0deB1OMoslXdZSfw+JgZYMx51kFkv7rKUDr8NJZrGky1rq7yEx0JLhuJPMYmmftXTgdTjJLJZ0WUv9PSQGWjIcd5JZLO2zlg68DieZxZIua6m/h8RAS4bjTjKLpX3W0oHX4SSzWNJlLfX3kBhoyXDcSWaxtM9aOvA6nGQWS7qspf4eEgMtGY47ySyW9llLB16Hk8xiSZe11N9DYqAlw3EnmcXSPmvpwOtwklks6bKW+ntIDLRkOO4ks1jaZy0deB1OMoslXdZSfw+JgZYMx51kFkv7rKUDr8NJZrGky1rq7yEx0JLhuJPMYmmftXTgdTjJLJZ0WUv9PSQGWjIcd5JZLO2zlg68DieZxZIua6m/h8RAS4bjTjKLpX3W0oHX4SSzWNJlLfX3kBhoyXDcSWaxtM9aOvAzP+dTf/WV8Vsnv2mPnnzWRb046c+BxXSh/PZFWe/Fp7ypF59y5kW9OPE8Gsq6pMta6u8hMdCSYaCRzGJpn7V04A20/jTQBrKy0TDxPBrKuqTLWurvITHQkmGgkcxiaZ+1dOANtP400AaystEw8Twayrqky1rq7yEx0JJhoJHMYmmftXTgDbT+NNAGsrLRMPE8Gsq6pMta6u8hMdAG5Oc//3m84x3viNWrV8fy5cvjqKOOivPOO2+fnmGgkcxiaZ9lOfB9driBtv8aaANZ2WiYeB4NZV3SZVn6OxsG2oCcffbZ8fCHPzw++MEPxhe+8IW45JJL4pWvfOU+PcNAI5nF0j7LcuD77HADbf810AaystEw8Twayrqky7L0dzYMtIG45pprYtmyZXHjjTfu13MMNJJZLO2zDAe+7w430PZfA20gKxsNE8+joaxLuixDf2fEQBuI5z3vefHMZz5zv59joJHMYmmfZTjwfXe4gbb/GmgDWdlomHgeDWVd0lInuUoAACAASURBVGUZ+jsjBtpAPPKRj4z/+l//a7z0pS+NQw45JB74wAfG6aefHrfccss+PcdAI5nF0j7LcOD77nADbf810AaystEw8TwayrqkyzL0d0YMtIFYvnx5rFy5Mn7jN34jPvWpT8Vf//Vfx2Me85h4whOeEPfcc88e3296ejq2bt0665YtWww0kincG5kOfN8dbqDtvwbaQFY2GiaeR0NZL0Sm/s6IgTYQBx54YKxYsSJuv/322de++tWvRtd18Td/8zd7fL+pqanoum43DTSStbs3Mh34vjvcQNt/DbSBrGw0TDyPhrJeiEz9nREDbSAOP/zweMpTnrLb64ceemj82Z/92R7fz/+DRjKreyPTge+7ww20/ddAG8jKRsPE82go64XI1N8ZMdAGYt26dXs87m95y1uKn+N70EhmsbTPMhz4vjvcQNt/DbSBrGw0TDyPhrIu6bIM/Z0RA20g3vnOd8ZBBx0U27Ztm31t5iupV199dfFzDDSSWSztswwHvu8ON9D2XwNtICsbDRPPo6GsS7osQ39nxEAbiOnp6XjkIx8ZJ598clx55ZXxyU9+Mo499thYu3Zt7Nq1q/g5BhrJLJb2WYYD33eHG2j7r4E2kJWNhonn0VDWJV2Wob8zYqANyDe/+c34D//hP8TBBx8chx56aJx55pnxgx/8YJ+eYaCRzGJpn2U58H12uIG2/xpoA1nZaJh4Hg1lXdJlWfo7GwZaMgw0klks7bOWDryB1p8G2kBWNhomnkdDWZd0WUv9PSQGWjIMNJJZLO2zlg68gdafBtpAVjYaJp5HQ1mXdFlL/T0kBloySgfazu8f14sTL6sEylrOY7OvrEv7rKUDP/NzXnfsi+K01a/Zo8e876JenPTn0mK6UH77oqzlXIunPe61vXjMey/qxZIua6m/h8RAS4aBVp+ylvPYNNAWDwOtPw0HOY9NAw0zGGjJMNDqU9ZyHpsG2uJhoPWn4SDnsWmgYQYDLRkGWn3KWs5j00BbPAy0/jQc5Dw2DTTMYKAlw0CrT1nLeWwaaIuHgdafhoOcx6aBhhkMtGQYaPUpazmPTQNt8TDQ+tNwkPPYNNAwg4GWDAOtPmUt57FpoC0eBlp/Gg5yHpsGGmYw0JJhoNWnrOU8Ng20xcNA60/DQc5j00DDDAZaMgy0+pS1nMemgbZ4GGj9aTjIeWwaaJjBQEuGgVafspbz2DTQFg8DrT8NBzmPTQMNMxhoyTDQ6lPWch6bBtriYaD1p+Eg57FpoGEGAy0ZBlp9ylrOY9NAWzwMtP40HOQ8Ng00zGCgJcNAq09Zy3lsGmiLh4HWn4aDnMemgYYZDLRkGGj1KWs5j00DbfEw0PrTcJDz2DTQMIOBlgwDrT5lLeexaaAtHgZafxoOch6bBhpmMNCSYaDVp6zlPDYNtMXDQOtPw0HOY9NAwwwGWjIMtPqUtZzHpoG2eBho/Wk4yHlsGmiYwUBLhoFWn7KW89g00BYPA60/DQc5j00DDTMYaMkoHWgkOWlL+6ylA6/DSWaxpMta6u8hMdCS4biTzGJpn7V04HU4ySyWdFlL/T0kBloyHHeSWSzts5YOvA4nmcWSLmupv4fEQEuG404yi6V91tKB1+Eks1jSZS3195AYaBPizjvvjKOPPjq6rot//Md/LH4/x51kFkv7LOOB1+Ekx25Jl2Xs7wwYaBPiNa95TRxxxBGOO8nRWtpnGQ+8Dic5dku6LGN/Z8BAmwA33XRTHHzwwfEXf/EXjjvJ0VraZ9kOvA4n2YIlXZatv7NgoE2AZzzjGfHKV74yNm/e7LiTHK2lfZbtwOtwki1Y0mXZ+jsLBtrAbNy4MY444oiYnp523EmO2tI+y3TgdTjJVizpskz9nQkDbUB27NgRq1atiosvvjgioui4T09Px9atW2fdsmWL404yhXsj24HX4SRbciGy9Xc2DLQB+ZM/+ZM46aSTYteuXRFRdtynpqai67rddNxJ1u7eyHbgdTjJllyIbP2dDQNtIL797W/H8uXL41Of+lRs3749tm/fHldffXV0XRdf/OIX484775z3/Xz1lWRW90amA6/DSbbmQmTq74wYaAMx85XWPblu3bqi5/j+BZJZLO2zDAdeh5NszZIuy9DfGTHQBmL79u2xefPmOb7nPe+JruviL//yL4u/ydxxJ5nF0j7LcOB1OMnWLOmyDP2dEQNtgvgTwEiO2dI+y3rgdTjJMVvSZVn7u3YMtAniuJMcs6V9lvXA63CSY7aky7L2d+0YaMlw3ElmsbTPWjrwOpxkFku6rKX+HhIDLRmOO8kslvZZSwdeh5PMYkmXtdTfQ2KgJcNxJ5nF0j5r6cDrcJJZLOmylvp7SAy0ZDjuJLNY2mctHXgdTjKLJV3WUn8PiYGWDMedZBZL+6ylA6/DSWaxpMta6u8hMdCS4biTzGJpn7V04HU4ySyWdFlL/T0kBloyHHeSWSzts5YOvA4nmcWSLmupv4fEQEuG404yi6V91tKB1+Eks1jSZS3195AYaMlw3ElmsbTPWjrwOpxkFku6rKX+HhIDLRmOO8kslvZZSwdeh5PMYkmXtdTfQ2KgJcNxJ5nF0j5r6cDrcJJZLOmylvp7SAy0ZDjuJLNY2mctHXgdTjKLJV3WUn8PiYGWDMedZBZL+6ylA6/DSWaxpMta6u8hMdCS4biTzGJpn7V04HU4ySyWdFlL/T0kBloyHHeSWSzts5YOvA4nmcWSLmupv4fEQEuG404yi6V91tKB1+Eks1jSZS3195AYaMlw3ElmsbTPWjrwOpxkFku6rKX+HhIDLRmOO8kslvZZSwdeh5PMYkmXtdTfQ2KgJcNxJ5nF0j5r6cDrcJJZLOmylvp7SAy0ZDjuJLNY2mctHXgdTjKLJV3WUn8PiYGWDMedZBZL+6ylA6/DSWaxpMta6u8hMdCS4biTzGJpn7V04HU4ySyWdFlL/T0kBloyHHeSWSzts5YOvA4nmcWSLmupv4fEQEuG404yi6V91tKB1+Eks1jSZS3195AYaAPx13/91/E7v/M7cfTRR8eKFSviSU96Unz84x+PXbt27dNzHHeSWSztswwHXoeTbM2SLsvQ3xkx0AbiN37jN+L3fu/34rLLLovPf/7z8cd//MexdOnSeMtb3rJPz3HcSWaxtM8yHHgdTrI1S7osQ39nxEAbiNtvv323184555x48IMfvE/PcdxJZrG0zzIceB1OsjVLuixDf2fEQJsgH/zgB6PruvjpT39a/D6OO8kslvZZ1gOvw0mO2ZIuy9rftWOgTZCzzjorHvWoR+3T+zjuJLNY2mdZD7wOJzlmS7osa3/XjoE2If7+7/8+li5dGhs2bFjw7aanp2Pr1q2zbtmyxXEnmcK9kfnA63CSY3chMvd3Bgy0CfDd7343jjrqqHj6058eO3fuXPBtp6amouu63XTcSdbu3sh64HU4yRZciKz9nQUDbWC2b98eJ5xwQqxZsybuuOOOvb69r76SzOreyHjgdTjJVlyIjP2dCQNtQH7605/GU5/61Fi1atX9/oT2/Qsks1jaZ1kOvA4n2ZIlXZalv7NhoA3EPffcE8961rPisMMOixtvvPF+P8dxJ5nF0j7LcOB1OMnWLOmyDP2dEQNtIM4555zoui7e9a53xQ033DDHu+66q/g5jjvJLJb2WYYDr8NJtmZJl2Xo74wYaAPxqEc9at5vFO+6Lm699dbi5zjuJLNY2mcZDrwOJ9maJV2Wob8zYqAlw3EnmcXSPmvpwOtwklks6bKW+ntIDLRkOO4ks1jaZy0deB1OMoslXdZSfw+JgZYMx51kFkv7rKUDr8NJZrGky1rq7yEx0JLhuJPMYmmftXTgZ37Op/7qK+O3Tn7THj359y/qxUl/DiyqT3lTL5581kW9OPE85Jzf3/yzXnzKmRf1YkmXtdTfQ2KgJcNAI5nF0j5r6cAbaD1qOMh5bBpouA8DLRkGGskslvZZSwfeQOtRw0HOY9NAw30YaMkw0EhmsbTPWjrwBlqPGg5yHpsGGu7DQEuGgUYyi6V91tKBN9B61HCQ89g00HAfBloyDDSSWSzts5YOvIHWo4aDnMemgYb7MNCSYaCRzGJpn7V04A20HjUc5Dw2DTTch4GWDAONZBZL+6ylA2+g9ajhIOexaaDhPgy0ZBhoJLNY2mctHXgDrUcNBzmPTQMN92GgJcNAI5nF0j5r6cAbaD1qOMh5bBpouA8DLRkGGskslvZZSwfeQOtRw0HOY9NAw30YaMkw0EhmsbTPWjrwBlqPGg5yHpsGGu7DQEuGgUYyi6V91tKBN9B61HCQ89g00HAfBloyDDSSWSzts5YOvIHWo4aDnMemgYb7MNCSYaCRzGJpn7V04A20HjUc5Dw2DTTch4GWDAONZBZL+6ylA2+g9ajhIOexaaDhPgy0ZBhoJLNY2mctHXgDrUcNBzmPTQMN92GgJcNAI5nF0j5r6cAbaD1qOMh5bBpouA8DLRmlA23n94/rxYmXVQJlLeex2VfWpX3W0oGf+TmvO/b8OO2xr96jx7zvol6c9OfSYnra6tf04jHvvagXJ52HnPNbW9YlXdZSfw+JgZYMA60+ZS3nsWmgLR4GWn/W9ovZSech5/zWlnVJl7XU30NioCXDQKtPWct5bBpoi4eB1p+1/WJ20nnIOb+1ZV3SZS3195AYaANy0003xTOe8YxYsWJFHHHEEfHqV7867r777n16hoFWn7KW89g00Oanzw430Pbf2n4xO+k85Jzf2rIu6bIs/Z0NA20gfvzjH8fDH/7weNrTnhaf+cxn4uKLL45DDz00XvKSl+zTcwy0+pS1nMemgbY7fXe4gbb/1vaL2UnnIef81pZ1SZdl6O+MGGgD8da3vjVWrlwZP/rRj2Zf+/CHPxwHHHBAfO973yt+joFWn7KW89g00Han7w430Pbf2n4xO+k85Jzf2rIu6bIM/Z0RA20gTj311HjOc54z57Xt27fHkiVL4uMf/3jxcwy0+pS1nMemgbY7fXe4gbb/1vaL2UnnIef81pZ1SZdl6O+MGGgD8bCHPSxe//rX7/b6UUcdFa997WuLn2Og1aes5Tw2DbTd6bvDDbT9t7ZfzE46Dznnt7asS7osQ39nxEAbiGXLlsU73/nO3V5/4hOfGOecc84e3296ejq2bt0665e//OXoui7Wdk+PU7rT9+h3/uGYXlzoY1DWch6nfWX9i901n1u2bImu6+K2225bzPrthb47/DdWPT/WHXv+Hn3Emy7oxUl/Li2m6459US8+4k8v6MVJ5yHn/NaW9Vj6OyMG2kDc3+M+NTUVXdeR5GjdsmXLYtZvL+hwktzdDP2dEQNtIO7vb4/55a++3nLLLfGFL3whbrvttr1+dbrkKx9btmzZr+dQ1rUo53w533bbbbFly5Z9/qPqJ4EOb1dZy3lM9pVzpv7OiIE2EKeeemo897nPnfPaHXfcsc/fYN4XW7f6vcNDIethkPMwtJqzDm8XWQ+DnIdBzjkw0AbirW99a/zKr/xKbN++ffa1v/zLv9znP6K5L/wPdDhkPQxyHoZWc9bh7SLrYZDzMMg5BwbaQMz8Jafr1q2La665Jj72sY/Fgx70oH3+S077wv9Ah0PWwyDnYWg1Zx3eLrIeBjkPg5xzYKANyL/8y7/E05/+9DjooIPi8MMPj1e96lUT+72709PTMTU1FdPT0xP5+C0h62GQ8zC0nLMObxNZD4Och0HOOTDQAAAAAKASDDQAAAAAqAQDDQAAAAAqwUADAAAAgEow0AAAAACgEgy0EXH11VfHk5/85Fi+fHk84hGPiDe+8Y3x85//fLe3+9GPfhQvetGL4sgjj4wHPOAB8djHPjb+4i/+Yva/33rrrdF13W6edtppQ/50qqUk5/nym/Hf/u3fZt/u7rvvjle96lVxxBFHxIoVK+IZz3hGfOMb3xj6p1QtfWY9339/3OMeN/RPqUpKcv75z38eb3/72+Nxj3tcHHTQQXHsscfGq1/96vjJT34y5+18Tt9/dPgw6PBh0N/DocPHh4E2Em644YZYunRpvOAFL4jPfOYz8a53vSsOOuigeOUrXznn7e6888444YQTYu3atfE//+f/jE2bNsUHP/jBeP/73z/7NjPH/e1vf3vccMMNs950001D/7SqozTnX8xtxsc+9rHxq7/6q3Pe7rzzzotDDz00Lr744vjMZz4Tp556ahx99NFxxx13DPnTqpK+s+66Ll72spfNebuvfe1rQ/6UqqQ05ze96U2xbNmyeNvb3habNm2KDRs2xMqVK+M//af/NOftfE7fP3T4MOjwYdDfw6HDx4mBNhJOO+20OPnkk+e8dtFFF8WBBx4YP/jBD2Zf+5M/+ZN4zGMeEz/96U/3+KyZ43755Zcv2o83K6U5/zIzmb7jHe+Yfe273/1uHHDAAfGRj3xk9rUf/ehHcfDBB8fb3/72/n/wyegz64h7D/x73vOeRfmxZqY058c97nFx9tlnz3m7N77xjXHQQQfNfqXW5/T9R4cPgw4fBv09HDp8nBhoI+Hwww+PN7zhDXNe+/rXvx5d18V//+//ffa1I444It761rcu+CzHfc+U5vzLXHjhhbFkyZL47ne/O/vaxRdfHEuWLInt27fPedvnPOc5sW7dul5/3BnpM+sIB35PlOb86Ec/Ol760pfOebt3vetd8YAHPCB27twZET6n9wcdPgw6fBj093Do8HFioI2EQw45JN785jfPee1f//Vfo+u6+OM//uOIiLjlllui67r40Ic+FKeffnosX748DjvssHjxi18856uxM8f9oQ99aCxdujQOP/zwePGLXxz/9//+30F/TjVSkvN8nHDCCbuV26tf/eo4+uijd3vb173udXHEEUf08uPNTJ9ZR9x74B/ykIfEAQccEA9+8IPj93//9+P73/9+3z/sdJTm/Na3vjUe/OAHx+c///m488474ytf+Uocc8wx8fKXv3z2bXxO3390+DDo8GHQ38Ohw8eJgTYSTjrppHj2s58957X/8T/+R3RdF+eee25ERFx//fXRdV2sXLky/vN//s9x7bXXxoYNG+Lggw+Oc845Z/b9/u3f/i1e/OIXx5VXXhmbN2+OCy+8MA4++OB42tOeFrt27Rr051UbJTn/Mv/0T/8UXdfFhz/84Tmv/5f/8l/iiU984m5v/853vjMOPPDA/n7QSekz64iIs88+OzZu3Bh/93d/Fxs2bIiHPexhsXr16tixY8ei/PizsC85X3DBBXO+Sf8FL3jB7FdeI3xO7w86fBh0+DDo7+HQ4ePEQBsJH//4x2PJkiXx53/+5/GjH/0o/v7v/z5WrVoVBxxwQJx33nkREXHddddF13W7/V7ld73rXbF06dLYtm3bHp9/6aWXRtd1ce211y7qz6N2SnL+ZV7zmtfEgQceGD/60Y/mvK4IF6bPrOdj5he7H/3oR/v+oaeiNOcNGzbEgx70oHjve98bf/d3fxcf/OAH4yEPecic3zLjc/r+o8OHQYcPg/4eDh0+Tgy0kbBr1654xSteEcuWLYuu62L58uXx1re+NR72sIfFn/7pn0ZExL/8y79E13Xx2te+ds77/uM//mN0XRdf/OIX9/j8n/3sZ3HAAQc0/02iJTn/8ts/8pGPjN/5nd/Z7b/5rQQL02fWe+Koo46KF73oRX3+sNNRkvP/+T//Jx7wgAfEBz7wgTnv+4lPfCKWLFkS3/zmNyPC5/T+oMOHQYcPg/4eDh0+Tgy0kXHHHXfEP/3TP8Udd9wR27Zti67r4jOf+UxE3Pt3Wyxfvny34/4P//AP0XVdfOlLX9rjc2eO+y//yUqtslDOv8gXv/jF6LouLrvsst3+28w34/7yH1373Oc+1zfj/gJ9ZL0njjrqqHjxi1/c5w83LQvl/JWvfCW6rovrrrtuzvvceOON0XVdfPazn40In9N9oMOHQYcPg/4eDh0+Lgy0EXPBBRfEscceO+cvKzz99NPjpJNOmvN2M//X9UK/reCSSy6Jruvi85///KL9eLMyX84znH/++bFy5cp5/0jsmT/O9hd/i8aPf/zjWLlyZfNf5d4T9zfr+fjSl74UXdfFxRdf3PcPMz2/nPMPfvCD6LouNmzYMOftZr7P4V//9V8jwud03+jwYdDhw6C/h0OH58dAGwnXX399vOMd74jPfvazceWVV8YLX/jCWL58+W7H+Ctf+UoceOCB8YIXvCCuueaaeM973hMrVqyIV73qVbNvMzU1FS972cvib/7mb+Laa6+NP/uzP4sVK1bEunXrmv8G89KcIyLuueeeeOhDHxrPf/7z9/i88847Lx70oAfFxz72sbjmmmti3bp1/kLI++gz6w9/+MPxwhe+MC677LLYtGlTvOc974mHPOQh8fjHP774FwRjpTTn//gf/2P8yq/8Slx00UWxadOm+PM///M47LDD4t//+38/5+18Tt8/dPgw6PBh0N/DocPHiYE2Ev73//7fcfLJJ8fBBx8cK1eujKc//elx/fXXz/u2n/3sZ+PXfu3XYvny5XHUUUfF6173urjnnntm//snP/nJOOmkk+LQQw+NZcuWxTHHHBOvetWr4ic/+clQP51q2Zec/9f/+l/RdV18+tOf3uPz7rrrrnjlK18Zhx9+eBx00EHxjGc8I2666abF+uGnos+sr7322vh3/+7fxWGHHRbLli2Lo446Ks4555wF/1CFVijNeXp6Ol71qlfFYx7zmHjgAx8Yxx57bPy3//bfdvv7cnxO3z90+DDo8GHQ38Ohw8eJgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCgAQAAAEAlGGgAAAAAUAkGGgAAAABUgoEGAAAAAJVgoAEAAABAJRhoAAAAAFAJBhoAAAAAVIKBBgAAAACVYKABAAAAQCUYaAAAAABQCQYaAAAAAFSCgQYAAAAAlWCg9cDNN98c5557bqxZsyaWLl0a69atm/ftvvOd78RZZ50VD3nIQ+KBD3xgnHDCCXH55ZcP+4MFAMxBhwMAasJA64ErrrgiVq1aFevXr4/Vq1fPe9y3bt0aRx99dDzzmc+MK664Iq699tp497vfHZdeeunwP2AAwCw6HABQEwZaD+zcuXP2n88444x5j/uZZ54Zp5xyypy3BQBMHh0OAKgJA61n5jvud9xxRxx44IG+0goAlaPDAQCTxkDrmfmO+6ZNm6LrurjsssvilFNOiWXLlsWRRx4ZF1xwQfz85z+fzA8UALAbOhwAMGkMtJ6Z77hfeuml0XVdHHLIIfHqV786Nm3aFG9+85tj2bJlceGFFy74vOnp6di6deust9xyS3zhC1+I2267bc7rJJnN2267LbZs2RJ33333IrbyvqHDSXLv1tjfY8JA65n5jvsll1wSXdfF+vXr57z+R3/0R3HIIYfErl279vi8qamp6LqOJEfrli1bFqOO7xc6nCTLram/x4SB1jPzHfdPf/rT0XVdfOhDH5rz+uWXXx5d18V3vvOdPT7vl7/6+uUvfzm6rou13dPjlO50kqzWvX0FdsuWLdF1Xdx2222LUcf3Cx1Okvearb/HhIHWM/Md91tvvXXe4/63f/u30XVdbN26tfj5W7duja7r4pTu9HjGkv+XJKu1tM/2pQMXGx1Okvda0mU19feYMNB6Zk9/RPMTn/jE+N3f/d05r73kJS+Jhz70ofv0fMedZBZL+6ymA6/DSfJeS7qspv4eEwZaD+zYsSM2btwYGzdujLVr18bxxx8/++/btm2LiHu/0rpkyZJ4+ctfHp/73OdiamoqDjjggHj/+9+/Tx/LcSeZxdI+m/SB1+EkubslXTbp/h4rBloPzPz2l/ncvHnz7Ntdcskl8YQnPCGWL18exx57bLz73e/e54/luJPMYmmfTfrA63CS3N2SLpt0f48VAy0ZjjvJLJb2WUsHXoeTzGJJl7XU30NioCXDcSeZxdI+a+nA63CSWSzpspb6e0gMtGSUHved3z+uFyddDhmUtZzHZl9Zl/ZZSwe+tMMfteGiXjxt9Wt6cdKfk4tpb1k/7rW9OOk8Wsl5zFkf876LerGky1rq7yEx0JJhoNWnrOU8Ng20xcNAq8/ahsOk82gl5zFnbaDlx0BLhoFWn7KW89g00BYPA60+axsOk86jlZzHnLWBlh8DLRkGWn3KWs5j00BbPAy0+qxtOEw6j1ZyHnPWBlp+DLRkGGj1KWs5j00DbfEw0OqztuEw6TxayXnMWRto+THQkmGg1aes5Tw2DbTFw0Crz9qGw6TzaCXnMWdtoOXHQEuGgVafspbz2DTQFg8DrT5rGw6TzqOVnMectYGWHwMtGQZafcpazmPTQFs8DLT6rG04TDqPVnIec9YGWn4MtGQYaPUpazmPTQNt8TDQ6rO24TDpPFrJecxZG2j5MdCSYaDVp6zlPDYNtMXDQKvP2obDpPNoJecxZ22g5cdAS4aBVp+ylvPYNNAWDwOtPmsbDpPOo5Wcx5y1gZYfAy0ZBlp9ylrOY9NAWzwMtPqsbThMOo9Wch5z1gZafgy0ZBho9SlrOY9NA23xMNDqs7bhMOk8Wsl5zFkbaPkx0JJhoNWnrOU8Ng20xcNAq8/ahsOk82gl5zFnbaDlx0BLhoFWRx88hQAAIABJREFUn7KW89g00BYPA60+axsOk86jlZzHnLWBlh8DLRkGWn3KWs5j00BbPAy0+qxtOEw6j1ZyHnPWBlp+DLRkGGj1KWs5j00DbfEw0OqztuEw6TxayXnMWRto+THQkmGg1aes5Tw2DbTFw0Crz9qGw6TzaCXnMWdtoOXHQEtG6XEnyUlb2mctHfjSDn/KmRf14m+t/dNenPTn0mL6lN+7qBd/6+Q39eKk82glZ1nv3ZIua6m/h8RAS4aBRjKLpX3W0oE30OqztuEw6TxayVnWBlrNGGjJMNBIZrG0z1o68AZafdY2HCadRys5y9pAqxkDLRkGGskslvZZSwfeQKvP2obDpPNoJWdZG2g1Y6D1wM033xznnnturFmzJpYuXRrr1q1b8O0vv/zy6LounvSkJ+3zxzLQSGaxtM8mfeBr7HADbThrGw6TzqOVnGVtoNWMgdYDV1xxRaxatSrWr18fq1evXvC4//SnP41jjjkmjjjiCAON5Kgt7bNJH/gaO9xAG87ahsOk82glZ1kbaDVjoPXAzp07Z//5jDPOWPC4X3DBBfG0pz0tzj77bAON5Kgt7bNJH/gaO9xAG87ahsOk82glZ1kbaDVjoPXMQsf9m9/8ZqxYsSK+9rWvGWgkR29pn9V04GvpcANtOGsbDpPOo5WcZW2g1YyB1jMLHffTTz89zj///IgIA43k6C3ts5oOfC0dbqANZ23DYdJ5tJKzrA20mjHQemZPx/2qq66KBz/4wXH77bdHRPlxn56ejq1bt866ZcsWA41kCvdGjQe+lg430IaztuEw6TxayVnWBlrNGGg9M99x/9nPfhaPfvSjY8OGDbOvlR73qamp6LpuNw00krW7N2o88LV0uIE2nLUNh0nn0UrOsjbQasZA65n5jvvb3va2OO644+L222+P7du3x/bt2+PMM8+ME044IbZv3x533333Hp/n/0EjmdW9UeOBr6XDDbThrG04TDqPVnKWtYFWMwZaz8x33M8+++x5v4I648c//vHi5/seNJJZLO2zmg58LR1uoA1nbcNh0nm0krOsDbSaMdB6Zr7jftNNN8XmzZvneNppp8VjHvOY2Lx5c3z/+98vfr6BRjKLpX1W04GvpcMNtOGsbThMOo9Wcpa1gVYzBloP7NixIzZu3BgbN26MtWvXxvHHHz/779u2bZv3ffwpjiTHbmmfTfrA19jhBtpw1jYcJp1HKznL2kCrGQOtB2699dY9/taXzZs3z/s+BhrJsVvaZ5M+8DV2uIE2nLUNh0nn0UrOsjbQasZAS4aBRjKLpX3W0oE30OqztuEw6TxayVnWBlrNGGjJMNBIZrG0z1o68AZafdY2HCadRys5y9pAqxkDLRkGGskslvZZSwdeh5PMYkmXtdTfQ2KgJcNxJ5nF0j5r6cDrcJJZLOmylvp7SAy0ZDjuJLNY2mctHXgdTjKLJV3WUn8PiYGWDMedZBZL+6ylA6/DSWaxpMta6u8hMdCS4biTzGJpn7V04HU4ySyWdFlL/T0kBloyHHeSWSzts5YOvA4nmcWSLmupv4fEQEuG404yi6V91tKB1+Eks1jSZS3195AYaMlw3ElmsbTPWjrwOpxkFku6rKX+HhIDLRmOO8kslvZZSwdeh5PMYkmXtdTfQ2KgJcNxJ5nF0j5r6cDrcJJZLOmylvp7SAy0ZDjuJLNY2mctHXgdTjKLJV3WUn8PiYGWDMedZBZL+6ylA6/DSWaxpMta6u8hMdCS4biTzGJpn7V04HU4ySyWdFlL/T0kBloyHHeSWSzts5YOvA4nmcWSLmupv4fEQEuG404yi6V91tKB1+Eks1jSZS3195AYaMlw3ElmsbTPWjrwOpxkFku6rKX+HhIDLRmOO8kslvZZSwdeh5PMYkmXtdTfQ2KgJcNxJ5nF0j5r6cDrcJJZLOmylvp7SAy0ZDjuJLNY2mctHXgdTjKLJV3WUn8PiYGWDMedZBZL+6ylA6/DSWaxpMta6u8hMdCS4biTzGJpn7V04HU4ySyWdFlL/T0kBloyHHeSWSzts5YOvA4nmcWSLmupv4fEQOuBm2++Oc4999xYs2ZNLF26NNatWzfnv09PT8fU1FSsXbs2Dj300Dj88MPj2c9+dvzzP//zPn8sx51kFkv7bNIHXoeT5O6WdNmk+3usGGg9cMUVV8SqVati/fr1sXr16t2O+9e//vU48sgj4/Wvf31cc801ceWVV8app54aBx98cHzjG9/Yp4/luJPMYmmfTfrA63CS3N2SLpt0f48VA60Hdu7cOfvPZ5xxxm7H/Sc/+Uns2LFjzmt33nlnHHbYYfGyl71snz6W404yi6V9NukDr8NJcndLumzS/T1WDLSeme+474mTTz45nve85+3T8x13klks7bOaDrwOJ8l7Lemymvp7TBhoPVN63Ldv3x4rVqyIqampfXq+404yi6V9VtOB1+Ekea8lXVZTf48JA61nSo/7OeecEytXrozvfe97C77d9PR0bN26ddYtW7Y47iRTuDdqPPA6nCTvdSFq7O8xYaD1TMlx/9jHPhZd18Vf/dVf7fV5U1NT0XXdbjruJGt3b9R44HU4Sd7rQtTY32PCQOuZvR33T3/607Fs2bK44IILip7nq68ks7o3ajzwOpwk73UhauzvMWGg9cxCx/2GG26IFStWxAtf+ML7/Xzfv0Ayi6V9VtOB1+Ekea8lXVZTf48JA61n9nTcb7zxxjjssMPiWc96Vtxzzz33+/mOO8kslvZZTQdeh5PkvZZ0WU39PSYMtB7YsWNHbNy4MTZu3Bhr166N448/fvbft23bFj/84Q/jEY94RBx99NHx+c9/Pm644YZZb7zxxn36WI47ySyW9tmkD7wOJ8ndLemySff3WDHQeuDWW2+d95vAu66LzZs3x+bNm/f430v/vp0ZHHeSWSzts0kfeB1Okrtb0mWT7u+xYqAlw3EnmcXSPmvpwOtwklks6bKW+ntIDLRkOO4ks1jaZy0deB1OMoslXdZSfw+JgZYMx51kFkv7rKUDr8NJZrGky1rq7yEx0JLhuJPMYmmftXTgdTjJLJZ0WUv9PSQGWjIcd5JZLO2zlg68DieZxZIua6m/h8RAS4bjTjKLpX3W0oHX4SSzWNJlLfX3kBhoyXDcSWaxtM9aOvA6nGQWS7qspf4eEgMtGY47ySyW9llLB16Hk8xiSZe11N9DYqAlw3EnmcXSPmvpwOtwklks6bKW+ntIDLRkOO4ks1jaZy0deB1OMoslXdZSfw+JgZYMx51kFkv7rKUDr8NJZrGky1rq7yEx0JLhuJPMYmmftXTgdTjJLJZ0WUv9PSQGWjIcd5JZLO2zlg68DieZxZIua6m/h8RAS4bjTjKLpX3W0oHX4SSzWNJlLfX3kBhoyXDcSWaxtM9aOvA6nGQWS7qspf4eEgMtGY47ySyW9llLB16Hk8xiSZe11N9DYqAlw3EnmcXSPmvpwOtwklks6bKW+ntIDLRkOO4ks1jaZy0deB1OMoslXdZSfw+JgZYMx51kFkv7rKUDr8NJZrGky1rq7yEx0JLhuJPMYmmftXTgdTjJLJZ0WUv9PSQGWjIcd5JZLO2zlg68DieZxZIua6m/h8RAS0bpcd/5/eN6cdLlkEFZy3ls9pV1aZ+1dOBLO/xRGy7qxdNWv6YXJ/05uZge876LelHWA+X8uNf25qQzqT3rki5rqb+HxEBLhoFWn7KW89g00BYPA60+DbRkORtog2Vd0mUt9feQGGjJMNDqU9ZyHpsG2uJhoNWngZYsZwNtsKxLuqyl/h4SA60Hbr755jj33HNjzZo1sXTp0li3bt28b/fRj340HvvYx8YDHvCAOPHEE+Pqq6/e549loNWnrOU8NlsbaDV2uIE2nAZaspwNtMGyLumySff3WDHQeuCKK66IVatWxfr162P16tXzHvdPfvKTsWTJknjDG94QmzZtivPOOy+WLVsWN9xwwz59LAOtPmUt57HZ2kCrscMNtOE00JLlbKANlnVJl026v8eKgdYDO3funP3nM844Y97jvnr16jjrrLPmvPabv/mb8du//dv79LEMtPqUtZzHZmsDrcYON9CG00BLlrOBNljWJV026f4eKwZaz8x33L/1rW9F13Vx5ZVXznn9fe97Xyxfvjzuuuuu4ucbaPUpazmPzdYG2i9SS4cbaMNpoCXL2UAbLOuSLqupv8eEgdYz8x33T33qU9F1Xdx8881zXv/sZz8bXdfFTTfdVPx8A60+ZS3nsWmgrZvz2iQ63EAbTgMtWc4G2mBZl3RZTf09Jgy0npnvuH/iE5+Iruvi9ttvn/P6li1bouu6uO666/b4vOnp6di6deusM+9joNWjrOU8Ng20dXNem0SHG2jDaaAly9lAGyzrhaixv8eEgdYzfR/3qamp6LpuNw20epS1nMemgbZuzmuT6HADbTgNtGQ5G2iDZb0QNfb3mDDQeqbv3x7j/0GrX1nLeWwaaOvmvDaJDjfQhtNAS5azgTZY1gtRY3+PCQOtZxb6BvOrrrpqzusbNmyI5cuXx9133138fN+DVp+ylvPYNNDWzXltEh1uoA2ngZYsZwNtsKxLuqym/h4TBlrPLPRHND//+c+f89pTn/pUf8z+CJS1nMemgbZut9eH7nADbTgNtGQ5G2iDZV3SZTX195gw0Hpgx44dsXHjxti4cWOsXbs2jj/++Nl/37ZtW0REXHrppbFkyZJ44xvfGJs3b47zzz8/li1bFtdff/0+fSwDrT5lLeex2dpAq7HDDbThNNCS5WygDZZ1SZdNur/HioHWA7feeuu83wTedV1s3rx59u0++tGPxnHHHRfLly+PNWvWxNVXX73PH8tAq09Zy3lstjbQauxwA204DbRkORtog2Vd0mWT7u+xYqAlw0CrT1nLeWy2NtCGxECrTwMtWc4G2mBZl3RZS/09JAZaMgy0+pS1nMemgbZ4GGj1aaAly9lAGyzrki5rqb+HxEBLhoFWn7KW89g00BYPA60+DbRkORtog2Vd0mUt9feQGGjJKD3uJDlpS/uspQNf2uEnn3VRL/7W2j/txUl/Li2mTznzol78rZPf1IuTzqOVnEed9e9d1IslXdZSfw+JgZYMA41kFkv7rKUDb6DVZ23DYdJ5tJLzqLM20NJjoCXDQCOZxdI+a+nAG2j1WdtwmHQereQ86qwNtPQYaMkw0EhmsbTPWjrwBlp91jYcJp1HKzmPOmsDLT0GWjIMNJJZLO2zlg68gVaftQ2HSefRSs6jztpAS4+BlgwDjWQWS/uspQNvoNVnbcNh0nm0kvOoszbQ0mOgJcNAI5nF0j5r6cAbaPVZ23CYdB6t5DzqrA209BhoyTDQSGaxtM9aOvAGWn3WNhwmnUcrOY86awMtPQZaMgw0klks7bOWDryBVp+1DYdJ59FKzqPO2kBLj4GWDAONZBZL+6ylA2+g1Wdtw2HSebSS86izNtDSY6Alw0AjmcXSPmvpwBto9VnbcJh0Hq3kPOqsDbT0GGjJMNBIZrG0z1o68AZafdY2HCadRys5jzprAy09BloyDDSSWSzts5YOvIFWn7UNh0nn0UrOo87aQEuPgZYMA41kFkv7rKUDb6DVZ23DYdJ5tJLzqLM20NJjoCXDQCOZxdI+a+nAG2j1WdtwmHQereQ86qwNtPQYaMkw0EhmsbTPWjrwBlp91jYcJp1HKzmPOmsDLT0GWjIMNJJZLO2zlg68gVaftQ2HSefRSs6jztpAS4+BlgwDjWQWS/uspQNvoNVnbcNh0nm0kvOoszbQ0mOgJcNAI5nF0j5r6cDrcJJZLOmylvp7SAy0ZDjuJLNY2mctHXgdTjKLJV3WUn8PiYGWDMedZBZL+6ylA6/DSWaxpMta6u8hMdAG5Morr4yTTz45Vq5cGUceeWQ873nPi1tuuWWfnuG4k8xiaZ9lOfA6nGRLlnRZlv7OhoE2EJs3b46lS5fGH/zBH8TnPve5uOyyy2L16tWxevXq+NnPflb8HMedZBZL+yzDgdfhJFuzpMsy9HdGDLSBOO+88+LYY4+NXbt2zb62adOm6Lourr/++uLnOO4ks1jaZxkOvA4n2ZolXZahvzNioA3EH/7hH8aJJ54457WvfvWr0XVdXHfddcXPcdxJZrG0zzIceB1OsjVLuixDf2fEQBuIL37xi7Fs2bL4wAc+EHfccUd861vfitNOOy1OOumk2LlzZ/FzHHeSWSztswwHXoeTbM2SLsvQ3xkx0AbkqquuipUrV0bXddF1XTz5yU+OH/7whwu+z/T0dGzdunXWLVu2OO4kU7g3sh14HU6yJRciW39nw0AbiC996UvxoAc9KF7xilfEpk2bYuPGjXHiiSfG2rVr46677trj+01NTc3+YuAXddxJ1u7eyHTgdTjJ1lyITP2dEQNtIH7913891q9fP+e17373u7FkyZK4+OKL9/h+vvpKMqt7I9OB1+EkW3MhMvV3Rgy0gTjooIPiwgsv3O31hz3sYfG6172u+Dm+f4FkFkv7LMOB1+EkW7OkyzL0d0YMtIF4/OMfH8997nPnvPbtb387lixZEh/5yEeKn+O4k8xiaZ9lOPA6nGRrlnRZhv7OiIE2EO9973uj67p46UtfOvuXnJ5wwgnx8Ic/PH784x8XP8dxJ5nF0j7LcOB1OMnWLOmyDP2dEQNtIHbt+v/bu9NYP8tq4cNPpw0dtIwtFhkqB9CAYwCNvtAP9IgRNcIH44TEYEXQOFFQUKyg9IQpeIoWUCvRqCCQGIIYIMyjpCoSAxhDQF6qErS0WxA1QNf7wbeNm9L2Lvz3vZ/1v68r+X3Z7O7Truyz7rN6aFkXF1xwQbzuda+LmTNnxi677BJHHHFE/O53v9uqr+Nxl5Sl0n2W4YG3wyW1Vskuy7C/M3KgJeNxl5Sl0n3W0gNvh0vKUskua2l/1+RAS8bjLilLpfuspQfeDpeUpZJd1tL+rsmBlozHXVKWSvdZSw+8HS4pSyW7rKX9XZMDLRmPu6Qsle6zlh54O1xSlkp2WUv7uyYHWjIed0lZKt1nLT3wdrikLJXsspb2d00OtGQ87pKyVLrPWnrg7XBJWSrZZS3t75ocaMl43CVlqXSftfTA2+GSslSyy1ra3zU50JLxuEvKUuk+a+mBt8MlZalkl7W0v2tyoCXjcZeUpdJ91tIDb4dLylLJLmtpf9fkQEvG4y4pS6X7rKUH3g6XlKWSXdbS/q7JgZaMx11Slkr3WUsPvB0uKUslu6yl/V2TAy0Zj7ukLJXus5YeeDtcUpZKdllL+7smB1oyHndJWSrdZy098Ha4pCyV7LKW9ndNDrRkPO6SslS6z1p64O1wSVkq2WUt7e+aHGjJeNwlZal0n7X0wNvhkrJUssta2t81OdCS8bhLylLpPmvpgbfDJWWpZJe1tL9rcqAl43GXlKXSfdbSA2+HS8pSyS5raX/X5EBLxuMuKUul+6ylB94Ol5Slkl3W0v6uyYGWjMddUpZK91lLD7wdLilLJbuspf1dkwMtGY+7pCyV7rOWHng7XFKWSnZZS/u7JgdaMh53SVkq3WctPfB2uKQsleyylvZ3TQ60ZDzukrJUus9aeuDtcElZKtllLe3vmhxoyXjcJWWpdJ+19MDb4ZKyVLLLWtrfNTnQkvG4S8pS6T5r6YG3wyVlqWSXtbS/a3KgVfTss8/GWWedFfvss0+MjIzEvHnz4thjj92qr+Fxl5Sl0n2W5YG3wyW1VMkuy7K/s3GgVXT00UfHK17xili+fHncfPPN8aMf/ShOOOGErfoaHndJWSrdZ1keeDtcUkuV7LIs+zsbB1ol1157bUydOjXuu+++l/R1PO6SslS6zzI88Ha4pNYq2WUZ9ndGDrRK3ve+98Xb3/72l/x1PO6SslS6zzI88Ha4pNYq2WUZ9ndGDrRKdt999/jUpz4Vn/70p+PlL395bLvttnH44YfHQw89tFVfx+MuKUul+yzDA2+HS2qtkl2WYX9n5ECrZGRkJGbNmhVvectb4uqrr47LLrss9tprr3jNa14TzzzzzCZ/3OjoaKxatWpDK1eu9LhLStGWZHrg7XBJrbU5mfZ3Rg60SqZNmxYzZsyIv/zlLxs+9qtf/Sq6rosrrrhikz9uyZIl0XXdRnncJfW9Lcn0wNvhklprczLt74wcaJXMmTMn3vzmN2/08dmzZ8fpp5++yR/nd18lZW1LMj3wdrik1tqcTPs7IwdaJQsWLNjk4/71r3+9+Ov48wuSslS6zzI88Ha4pNYq2WUZ9ndGDrRKzj777Jg+fXo8/vjjGz62/ndSr7rqquKv43GXlKXSfZbhgbfDJbVWyS7LsL8zcqBVMjo6GrvvvnscdNBBceWVV8Yll1wS8+fPjwMPPDDWrVtX/HU87pKyVLrPMjzwdrik1irZZRn2d0YOtIoefPDBeOc73xkzZ86M2bNnxwc+8IF47LHHtupreNwlZal0n2V54O1wSS1Vssuy7O9sHGjJeNwlZal0n7X0wNvhkrJUssta2t81OdCS8bhLylLpPmvpgbfDJWWpZJe1tL9rcqAl43GXlKXSfdbSA2+HS8pSyS5raX/X5EBLpvRxf+7P/zWQJno5ZMiszXnYGtSsS/dZSw986Q7f83/PGUiH7XPSQJro78nxzKzbnPNQz/ob5wykkl3W0v6uyYGWjAOtf5m1OQ9bDrTx40DrX2bd5pyHetYOtPQcaMk40PqXWZvzsOVAGz8OtP5l1m3Oeahn7UBLz4GWjAOtf5m1OQ9bDrTx40DrX2bd5pyHetYOtPQcaMk40PqXWZvzsOVAGz8OtP5l1m3Oeahn7UBLz4GWjAOtf5m1OQ9bDrTx40DrX2bd5pyHetYOtPQcaMk40PqXWZvzsOVAGz8OtP5l1m3Oeahn7UBLz4GWjAOtf5m1OQ9bDrTx40DrX2bd5pyHetYOtPQcaMk40PqXWZvzsOVAGz8OtP5l1m3Oeahn7UBLz4GWjAOtf5m1OQ9bDrTx40DrX2bd5pyHetYOtPQcaMk40PqXWZvzsOVAGz8OtP5l1m3Oeahn7UBLz4GWjAOtf5m1OQ9bDrTx40DrX2bd5pyHetYOtPQcaMk40PqXWZvzsOVAGz8OtP5l1m3Oeahn7UBLz4GWjAOtf5m1OQ9bDrTx40DrX2bd5pyHetYOtPQcaMk40PqXWZvzsOVAGz8OtP5l1m3Oeahn7UBLz4GWjAOtf5m1OQ9bDrTx40DrX2bd5pyHetYOtPQcaMk40PqXWZvzsOVAGz8OtP5l1m3Oeahn7UBLz4GWjAOtf5m1OQ9bDrTx40DrX2bd5pyHetYOtPQcaMmUPu6SNNGV7rOWHvjSHX7QB88ZSP990GkDaaK/l8Yzs67Tmz9wzkBa+ObTBlcP5tLnWZfsspb2d00OtGQcaJKyVLrPWnrgHWj9y6zr5EDLN+uSXdbS/q7JgZaMA01Slkr3WUsPvAOtf5l1nRxo+WZdssta2t81OdCScaBJylLpPmvpgXeg9S+zrpMDLd+sS3ZZS/u7JgfaBHnyySdj1113ja7r4p577in+cQ40SVkq3WcZH/jx3uGOhnqZdZ0caPlmXbLLMu7vDBxoE+Skk06KuXPnOtAkDW2l+yzjAz/eO9zRUC+zrpMDLd+sS3ZZxv2dgQNtAjzwwAMxc+bMuPDCCx1okoa20n2W7YGvscMdDfUy6zo50PLNumSXZdvfWTjQJsDChQvjhBNOiJtuusmBJmloK91n2R74Gjvc0VAvs66TAy3frEt2Wbb9nYUDrbLLL7885s6dG6Ojow40SUNd6T7L9MDX2uGOhnqZdZ0caPlmXbLLMu3vTBxoFf3973+P3XbbLVasWBERUfS4j46OxqpVqza0cuVKB5qkFG1Jtge+5g53NNTLrOvkQMs3683Jtr+zcaBVdPLJJ8cBBxwQ69ati4iyx33JkiXRdd1GOdAk9b0tyfbA19zhjoZ6mXWdHGj5Zr052fZ3Ng60Sv7whz/EyMhIXH311bFmzZpYs2ZNXHXVVdF1Xdx6663x5JNPvuCP8/9Bk5S1Lcn0wNfe4Y6Gepl1nRxo+Wa9OZn2d0YOtErW/07rplqwYEHR1/Fn0CRlqXSfZXjga+9wR0O9zLpODrR8sy7ZZRn2d0YOtErWrFkTN91005jOO++86LouvvOd7xT/IXMHmqQsle6zDA987R3uaKiXWdfJgZZv1iW7LMP+zsiBNoH8LY6ShrnSfZb1gfe3OA5HZl0nB1q+WZfssqz7u+8caBPIgSZpmCvdZ1kfeAfacGTWdXKg5Zt1yS7Lur/7zoGWjANNUpZK91lLD7wDrX+ZdZ0caPlmXbLLWtrfNTnQknGgScpS6T5r6YF3oPUvs66TAy3frEt2WUv7uyYHWjIONElZKt1nLT3wDrT+ZdZ1cqDlm3XJLmtpf9fkQEvGgSYpS6X7rKUH3g6XlKWSXdbS/q7JgZaMx11Slkr3WUsPvB0uKUslu6yl/V2TAy0Zj7ukLJXus5YeeDtcUpZKdllL+7smB1oyHndJWSrdZy098Ha4pCyV7LKW9ndNDrRkPO6SslS6z1p64O1wSVkq2WUt7e+aHGjJeNwlZal0n7X0wNvhkrJUssta2t81OdCS8bhLylLpPmvpgbfDJWWpZJe1tL9rcqAl43GXlKXSfdbSA2+HS8pSyS5raX/X5EBLxuMuKUul+6ylB94Ol5Slkl3W0v6uyYGWjMddUpZK91lLD7wdLilLJbuspf1dkwMtGY+7pCyV7rOWHng7XFKWSnZZS/u7JgdaMh53SVkq3WctPfB2uKQsleyylvZ3TQ60ZDzukrJUus9aeuDtcElZKtllLe3vmhxoyXjcJWWpdJ+19MDb4ZKyVLLLWtrfNTnQkvG4S8pS6T5r6YG3wyVlqWSXtbS/a3KgJeNxl5Sl0n3W0gNvh0vKUskua2l/1+RAS8bjLilLpfuspQfeDpeUpZJd1tL+rskrKXtVAAASm0lEQVSBlozHXVKWSvdZSw+8HS4pSyW7rKX9XZMDLRmPu6Qsle6zlh54O1xSlkp2WUv7uyYHWjIed0lZKt1nLT3wdrikLJXsspb2d00OtGQ87pKyVLrPWnrg7XBJWSrZZS3t75ocaJVcdtll8Z73vCd23XXXmDFjRrz+9a+Piy++ONatW7dVX8fjLilLpfsswwNvh0tqrZJdlmF/Z+RAq+Qtb3lLvP/9749LL700brjhhvjiF78YkydPjq9//etb9XU87pKyVLrPMjzwdrik1irZZRn2d0YOtEr+8pe/bPSxRYsWxfbbb79VX8fjLilLpfsswwNvh0tqrZJdlmF/Z+RAm0DLly+Pruvi6aefLv4xHndJWSrdZ1kfeDtc0jBXssuy7u++c6BNoA9+8IOxxx57bNWP8bhLylLpPsv6wNvhkoa5kl2WdX/3nQNtgtx2220xefLkWLZs2WY/b3R0NFatWrWhlStXetwlpWhLMj/wdrikYW9zMu/vDBxoE+DRRx+NefPmxaGHHhrPPffcZj93yZIl0XXdRnncJfW9Lcn6wNvhklpoc7Lu7ywcaJWtWbMm9t9//3jta18ba9eu3eLn+91XSVnbkowPvB0uqZU2J+P+zsSBVtHTTz8db3vb22K33XZ70d/Q/vyCpCyV7rMsD7wdLqmlSnZZlv2djQOtkmeeeSbe9a53xQ477BD33Xffi/46HndJWSrdZxkeeDtcUmuV7LIM+zsjB1olixYtiq7r4txzz4277rprTP/85z+Lv47HXVKWSvdZhgfeDpfUWiW7LMP+zsiBVskee+zxgn9QvOu6ePjhh4u/jsddUpZK91mGB94Ol9RaJbssw/7OyIGWjMddUpZK91lLD7wdLilLJbuspf1dkwMtGY+7pCyV7rOWHng7XFKWSnZZS/u7JgdaMh53SVkq3WctPfB2uKQsleyylvZ3TQ60ZDzukrJUus9aeuDtcElZKtllLe3vmhxoyXjcJWWpdJ+19MDb4ZKyVLLLWtrfNTnQkvG4S8pS6T5r6YG3wyVlqWSXtbS/a3KgJeNxl5Sl0n3W0gNvh0vKUskua2l/1+RAS8bjLilLpfuspQfeDpeUpZJd1tL+rsmBlozHXVKWSvdZSw+8HS4pSyW7rKX9XZMDLRmPu6Qsle6zlh54O1xSlkp2WUv7uyYHWjIed0lZKt1nLT3wdrikLJXsspb2d00OtGQ87pKyVLrPWnrg7XBJWSrZZS3t75ocaMl43CVlqXSftfTA2+GSslSyy1ra3zU50JLxuEvKUuk+a+mBt8MlZalkl7W0v2tyoCXjcZeUpdJ91tIDb4dLylLJLmtpf9fkQEvG4y4pS6X7rKUH3g6XlKWSXdbS/q7JgZaMx11Slkr3WUsPvB0uKUslu6yl/V2TAy0Zj7ukLJXus5YeeDtcUpZKdllL+7smB1oyHndJWSrdZy098Ha4pCyV7LKW9ndNDrRkPO6SslS6z1p64O1wSVkq2WUt7e+aHGjJeNwlZal0n7X0wNvhkrJUssta2t81OdCSKX3cn/vzfw2kiV4OGTJrcx62BjXr0n3W0gNfusP3/N9zBtJhe584kCb6e3I82/Mb5wykw/Y5aSBN9DxambNZb7mSXdbS/q7JgZaMA61/mbU5D1sOtPHjQOtffTscJnoerczZrB1ofeZAS8aB1r/M2pyHLQfa+HGg9a++HQ4TPY9W5mzWDrQ+c6BV9MADD8TChQtjxowZMXfu3DjxxBPjX//611Z9DQda/zJrcx62HGgvrOYOd6DVq2+Hw0TPo5U5m7UDrc8caJU88cQT8YpXvCIOOeSQuOaaa2LFihUxe/bs+OQnP7lVX8eB1r/M2pyHLQfaxmrvcAdavfp2OEz0PFqZs1k70PrMgVbJ0qVLY9asWbF69eoNH7voootiypQp8cc//rH46zjQ+pdZm/Ow5UDbWO0d7kCrV98Oh4meRytzNmsHWp850Co5+OCD44gjjhjzsTVr1sSkSZPi4osvLv46DrT+ZdbmPGw50DZWe4c70OrVt8NhoufRypzN2oHWZw60Snbeeef40pe+tNHH582bF1/4wheKv44DrX+ZtTkPWw60jdXe4Q60evXtcJjoebQyZ7N2oPWZA62SqVOnxtlnn73Rx/fbb79YtGjRJn/c6OhorFq1akO/+MUvouu6OLA7NP5Pd/gm+7+/3nMgbe5/hszanIezQc36P3fXC7Vy5croui4eeeSR8Vy/A1F7h7/ytFMH0oL5nxhIE/09OZ698qunDqQF848bSBM9j1bmbNZbblj2d0YOtEpe7OO+ZMmS6LpOkoa2lStXjuf6HQg7XJI2LsP+zsiBVsmL/ddjnv+7rw899FDcfPPN8cgjj2zxd6dLfudj5cqVL+nryKz7kjnnm/MjjzwSK1eu3Oq/qn4i2OHtZtbmPEwNas6Z9ndGDrRKDj744DjyyCPHfGzt2rVb/QfMB2XVKv/ucC1mXYc519HqnO3wdpl1HeZchznn4ECrZOnSpfGyl70s1qxZs+Fj3/nOd7b6r2geFP8LWo9Z12HOdbQ6Zzu8XWZdhznXYc45ONAqWf8fOV2wYEFce+218b3vfS+22267rf6PnA6K/wWtx6zrMOc6Wp2zHd4us67DnOsw5xwcaBXdf//9ceihh8b06dNjzpw5sXjx4gn7d3dHR0djyZIlMTo6OiH/81ti1nWYcx0tz9kOb5NZ12HOdZhzDg40AACAnnCgAQAA9IQDDQAAoCccaAAAAD3hQAMAAOgJB9oQueqqq+KNb3xjjIyMxCtf+cr4yle+Es8+++xGn7d69eo47rjjYpdddoltttkm9t5777jwwgs3/POHH344uq7bqMMOO6zmL6e3Sub8QvNb35/+9KcNn/evf/0rFi9eHHPnzo0ZM2bEwoUL43e/+13tX1JvDXLWL/TP991339q/pF4qmfOzzz4bZ555Zuy7774xffr0mD9/fpx44onx1FNPjfk839Mvnh1ehx1eh/1djx0+fBxoQ+Kuu+6KyZMnx1FHHRXXXHNNnHvuuTF9+vQ44YQTxnzek08+Gfvvv38ceOCB8ZOf/CRuvPHGWL58eXzzm9/c8DnrH/czzzwz7rrrrg098MADtX9ZvVM65/+c2/r23nvveMMb3jDm84499tiYPXt2rFixIq655po4+OCDY9ddd421a9fW/GX10qBn3XVdfPaznx3zeb/5zW9q/pJ6qXTOp512WkydOjX+53/+J2688cZYtmxZzJo1Kz7ykY+M+Tzf0y+OHV6HHV6H/V2PHT6cHGhD4rDDDouDDjpozMfOOeecmDZtWjz22GMbPnbyySfHXnvtFU8//fQmv9b6x/2nP/3puP18syqd8/Otn+lZZ5214WOPPvpoTJkyJb797W9v+Njq1atj5syZceaZZw7+J5/MIGcd8e8H/rzzzhuXn2tmpXPed9994+ijjx7zeV/5yldi+vTpG36n1vf0i2eH12GH12F/12OHDycH2pCYM2dOfPnLXx7zsd/+9rfRdV18//vf3/CxuXPnxtKlSzf7tTzum1Y65+c744wzYtKkSfHoo49u+NiKFSti0qRJsWbNmjGfe8QRR8SCBQsG+vPOaJCzjvDAb0rpnF/1qlfFpz/96TGfd+6558Y222wTzz33XET4nn4p7PA67PA67O967PDh5EAbEi9/+cvja1/72piP/f73v4+u6+KLX/xiREQ89NBD0XVdXHDBBXH44YfHyMhI7LDDDnH88ceP+d3Y9Y/7TjvtFJMnT445c+bE8ccfH3/729+q/pr6qGTOL2T//fffaLmdeOKJseuuu270uaecckrMnTt3ID/fzAY564h/P/A77rhjTJkyJbbffvv40Ic+FH/+858H/dNOp3TOS5cuje233z5uuOGGePLJJ+Puu++OPffcMz73uc9t+Bzf0y+eHV6HHV6H/V2PHT6cHGhD4oADDoh3v/vdYz72gx/8ILqui49//OMREXHnnXdG13Uxa9as+OhHPxrXX399LFu2LGbOnBmLFi3a8OP+9Kc/xfHHHx9XXnll3HTTTXHGGWfEzJkz45BDDol169ZV/XX1Tcmcn+/ee++NruvioosuGvPxj33sY7Hffvtt9Plnn312TJs2bXA/6aQGOeuIiKOPPjouv/zyuOWWW2LZsmWx8847xz777BN///vfx+Xnn8XWzPnUU08d84f0jzrqqA2/8xrhe/qlsMPrsMPrsL/rscOHkwNtSFx88cUxadKkOP/882P16tVx2223xW677RZTpkyJY489NiIi7rjjjui6bqN/V/ncc8+NyZMnx+OPP77Jr//jH/84uq6L66+/flx/HX1XMufnO+mkk2LatGmxevXqMR+3CDdvkLN+Iev/j93vfve7g/6pp1I652XLlsV2220X3/jGN+KWW26J5cuXx4477jjmX5nxPf3i2eF12OF12N/12OHDyYE2JNatWxef//znY+rUqdF1XYyMjMTSpUtj5513jq9+9asREXH//fdH13XxhS98YcyPveeee6Lrurj11ls3+fX/8Y9/xJQpU5r/Q6Ilc37+5+++++7xnve8Z6N/5l8l2LxBznpT5s2bF8cdd9wgf9rplMz5r3/9a2yzzTbxrW99a8yP/eEPfxiTJk2KBx98MCJ8T78Udngddngd9nc9dvhwcqANmbVr18a9994ba9eujccffzy6rotrrrkmIv7937YYGRnZ6HH/9a9/HV3Xxe23377Jr7v+cX/+36zUqs3N+T/deuut0XVdXHrppRv9s/V/GPf5f3XtkUce6Q/j/odBzHpT5s2bF8cff/wgf7ppbW7Od999d3RdF3fccceYH3PfffdF13Vx3XXXRYTv6UGww+uww+uwv+uxw4eLA22InXrqqTF//vwx/7HCww8/PA444IAxn7f+/3W9uX+t4Ec/+lF0XRc33HDDuP18s3qhOa/3iU98ImbNmvWCfyX2+r/O9j//FY0nnngiZs2a1fzvcm/Ki531C7n99tuj67pYsWLFoH+a6T1/zo899lh0XRfLli0b83nr/5zD73//+4jwPT1odngddngd9nc9dnh+DrQhceedd8ZZZ50V1113XVx55ZVxzDHHxMjIyEaP8d133x3Tpk2Lo446Kq699to477zzYsaMGbF48eINn7NkyZL47Gc/G1dccUVcf/31cfrpp8eMGTNiwYIFzf8B89I5R0Q888wzsdNOO8WHP/zhTX69Y489Nrbbbrv43ve+F9dee20sWLDAfxDy/xvkrC+66KI45phj4tJLL40bb7wxzjvvvNhxxx3j1a9+dfH/QTCsSuf83ve+N172spfFOeecEzfeeGOcf/75scMOO8Q73vGOMZ/ne/rFscPrsMPrsL/rscOHkwNtSPzyl7+Mgw46KGbOnBmzZs2KQw89NO68884X/Nzrrrsu3vSmN8XIyEjMmzcvTjnllHjmmWc2/PNLLrkkDjjggJg9e3ZMnTo19txzz1i8eHE89dRTtX45vbU1c/7Zz34WXdfFz3/+801+vX/+859xwgknxJw5c2L69OmxcOHCeOCBB8brp5/KIGd9/fXXx1vf+tbYYYcdYurUqTFv3rxYtGjRZv9ShVaUznl0dDQWL14ce+21V2y77bYxf/78+MxnPrPRfy/H9/SLY4fXYYfXYX/XY4cPJwcaAABATzjQAAAAesKBBgAA0BMONAAAgJ5woAEAAPSEAw0AAKAnHGgAAAA94UADAADoCQcaAABATzjQAAAAesKBBgAA0BMONAAAgJ5woAEAAPSEAw0AAKAnHGgAAAA94UADAADoCQcaAABATzjQAAAAesKBBgAA0BMONAAAgJ5woAEAAPSEAw0AAKAnHGgAAAA94UADAADoCQcaAABATzjQAAAAesKBBgAA0BMONAAAgJ5woAEAAPSEAw0AAKAnHGgAAAA94UADAADoCQcaAABATzjQAAAAesKBBgAA0BMONAAAgJ5woAEAAPSEAw0AAKAnHGgAAAA94UADAADoCQcaAABATzjQAAAAesKBBgAA0BMONAAAgJ5woAEAAPSEAw0AAKAnHGgAAAA94UADAADoCQcaAABATzjQAAAAesKBBgAA0BMONAAAgJ5woAEAAPSEAw0AAKAnHGgAAAA94UADAADoCQcaAABATzjQAAAAesKBBgAA0BMONAAAgJ5woAEAAPSEAw0AAKAnHGgAAAA94UADAADoCQcaAABATzjQAAAAesKBBgAA0BMONAAAgJ5woAEAAPSEAw0AAKAnHGgAAAA94UADAADoCQcaAABATzjQAAAAesKBBgAA0BMONAAAgJ5woAEAAPSEAw0AAKAnHGgAAAA94UADAADoCQcaAABAT/w/UzI6uDAnbzQAAAAASUVORK5CYII=\" width=\"799.3333095113444\">"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x7f71963acf98>"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dummy_image = numpy.zeros(mask.shape, dtype=\"float32\")\n",
"dummy_image[::5,::5] = 1\n",
"#dummy_image[mask] = -1\n",
"csr = csr_matrix(pre_csr)\n",
"dummy_blurred = csr.T.dot(dummy_image.ravel()).reshape(mask.shape)\n",
"fix, ax = subplots(2,2, figsize=(8,8))\n",
"ax[0,0].imshow(dummy_image)\n",
"ax[0,1].imshow(csr.dot(dummy_image.ravel()).reshape(mask.shape))\n",
"ax[1,1].imshow(csr.T.dot(dummy_image.ravel()).reshape(mask.shape))\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0, 16)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ax[0,0].set_xlim(964,981)\n",
"ax[0,0].set_ylim(0,16)\n",
"ax[0,1].set_xlim(964,981)\n",
"ax[0,1].set_ylim(0,16)\n",
"ax[1,1].set_xlim(964,981)\n",
"ax[1,1].set_ylim(0,16)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 672 ms, sys: 28 ms, total: 700 ms\n",
"Wall time: 699 ms\n",
"(1, 31, 0.00036844599621312167, 4.0083979131992555e-05, 2.1620353817016356, 4.873728452382846, 195.49168529213352)\n"
]
}
],
"source": [
"blured = csr.T.dot(dummy_image.ravel())\n",
"\n",
"# Invert this matrix: see https://arxiv.org/abs/1006.0758\n",
"\n",
"%time res = linalg.lsmr(csr.T, blured)\n",
"\n",
"restored = res[0].reshape(mask.shape)\n",
"ax[1,0].imshow(restored)\n",
"ax[1,0].set_xlim(964,981)\n",
"ax[1,0].set_ylim(0,16)\n",
"\n",
"print(res[1:])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conclusion of the raytracing part:\n",
"\n",
"We are able to simulate the path and the absorption of the photon in the thickness of the detector. \n",
"Numba helped substentially to make the raytracing calculation much faster. \n",
"The signal of each pixel is indeed spread on the neighboors, depending on the position of the PONI and this effect can be inverted using sparse-matrix pseudo-inversion.\n",
"\n",
"We will now save this sparse matrix to file in order to be able to re-use it in next notebook. But before saving it, it makes sense to spend some time in generating a high quality sparse matrix in throwing thousand rays per pixel in a grid of 32x32."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 5min 33s, sys: 1.68 s, total: 5min 34s\n",
"Wall time: 5min 32s\n"
]
}
],
"source": [
"%time pre_csr = thick.calc_csr(32)\n",
"hq_csr = csr_matrix(pre_csr)\n",
"from scipy.sparse import save_npz\n",
"save_npz(\"csr.npz\",hq_csr)"
]
}
],
"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.5.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}