diff --git a/examples/.ipynb_checkpoints/run-model-from-bmi-checkpoint.ipynb b/examples/.ipynb_checkpoints/run-model-from-bmi-checkpoint.ipynb new file mode 100644 index 0000000..c4ffdd0 --- /dev/null +++ b/examples/.ipynb_checkpoints/run-model-from-bmi-checkpoint.ipynb @@ -0,0 +1,386 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Run the `Heat` model through its BMI" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`Heat` models the diffusion of temperature on a uniform rectangular plate with Dirichlet boundary conditions. View the source code for the [model](https://github.com/csdms/bmi-example-python/blob/master/heat/heat.py) and its [BMI](https://github.com/csdms/bmi-example-python/blob/master/heat/bmi_heat.py) on GitHub." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Start by importing `os`, `numpy` and the `Heat` BMI:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "\n", + "from heat import BmiHeat" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create an instance of the model's BMI." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "x = BmiHeat()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What's the name of this model?" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The 2D Heat Equation\n" + ] + } + ], + "source": [ + "print(x.get_component_name())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Start the `Heat` model through its BMI using a configuration file:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "# Heat model configuration\n", + "shape:\n", + " - 6\n", + " - 8\n", + "spacing:\n", + " - 1.0\n", + " - 1.0\n", + "origin:\n", + " - 0.0\n", + " - 0.0\n", + "alpha: 1.0\n" + ] + } + ], + "source": [ + "cat heat.yaml" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "x.initialize(\"heat.yaml\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Check the time information for the model." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Start time: 0.0\n", + "End time: 1.7976931348623157e+308\n", + "Current time: 0.0\n", + "Time step: 0.25\n", + "Time units: s\n" + ] + } + ], + "source": [ + "print(\"Start time:\", x.get_start_time())\n", + "print(\"End time:\", x.get_end_time())\n", + "print(\"Current time:\", x.get_current_time())\n", + "print(\"Time step:\", x.get_time_step())\n", + "print(\"Time units:\", x.get_time_units())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Show the input and output variables for the component (aside on [Standard Names](https://csdms.colorado.edu/wiki/CSDMS_Standard_Names)):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(x.get_input_var_names())\n", + "print(x.get_output_var_names())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, get the identifier for the grid on which the temperature variable is defined:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "grid_id = x.get_var_grid(\"plate_surface__temperature\")\n", + "print(\"Grid id:\", grid_id)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then get the grid attributes:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Grid type:\", x.get_grid_type(grid_id))\n", + "\n", + "rank = x.get_grid_rank(grid_id)\n", + "print(\"Grid rank:\", rank)\n", + "\n", + "shape = np.ndarray(rank, dtype=int)\n", + "x.get_grid_shape(grid_id, shape)\n", + "print(\"Grid shape:\", shape)\n", + "\n", + "spacing = np.ndarray(rank, dtype=float)\n", + "x.get_grid_spacing(grid_id, spacing)\n", + "print(\"Grid spacing:\", spacing)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These commands are made somewhat un-Pythonic by the generic design of the BMI." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Through the model's BMI, zero out the initial temperature field, except for an impulse near the middle.\n", + "Note that *set_value* expects a one-dimensional array for input." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "temperature = np.zeros(shape)\n", + "temperature[3, 4] = 100.0\n", + "x.set_value(\"plate_surface__temperature\", temperature)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Check that the temperature field has been updated. Note that *get_value* expects a one-dimensional array to receive output." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "temperature_flat = np.empty_like(temperature).flatten()\n", + "x.get_value(\"plate_surface__temperature\", temperature_flat)\n", + "print(temperature_flat.reshape(shape))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now advance the model by a single time step:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x.update()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "View the new state of the temperature field:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x.get_value(\"plate_surface__temperature\", temperature_flat)\n", + "print(temperature_flat.reshape(shape))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There's diffusion!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Advance the model to some distant time:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "distant_time = 2.0\n", + "while x.get_current_time() < distant_time:\n", + " x.update()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "View the final state of the temperature field:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.set_printoptions(formatter={\"float\": \"{: 5.1f}\".format})\n", + "x.get_value(\"plate_surface__temperature\", temperature_flat)\n", + "print(temperature_flat.reshape(shape))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that temperature isn't conserved on the plate:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(temperature_flat.sum())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "End the model:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x.finalize()" + ] + } + ], + "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.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/.ipynb_checkpoints/run-model-from-bmi.-wth-different-states-checkpoint.ipynb b/examples/.ipynb_checkpoints/run-model-from-bmi.-wth-different-states-checkpoint.ipynb new file mode 100644 index 0000000..28b2b7e --- /dev/null +++ b/examples/.ipynb_checkpoints/run-model-from-bmi.-wth-different-states-checkpoint.ipynb @@ -0,0 +1,713 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Run the `Heat` model through its BMI" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`Heat` models the diffusion of temperature on a uniform rectangular plate with Dirichlet boundary conditions. View the source code for the [model](https://github.com/csdms/bmi-example-python/blob/master/heat/heat.py) and its [BMI](https://github.com/csdms/bmi-example-python/blob/master/heat/bmi_heat.py) on GitHub." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Start by importing `os`, `numpy` and the `Heat` BMI:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from heat import BmiHeat" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create an instance of the model's BMI." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "x = BmiHeat()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What's the name of this model?" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The 2D Heat Equation\n" + ] + } + ], + "source": [ + "print(x.get_component_name())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Start the `Heat` model through its BMI using a configuration file:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "# Heat model configuration\n", + "shape:\n", + " - 6\n", + " - 8\n", + "spacing:\n", + " - 1.0\n", + " - 1.0\n", + "origin:\n", + " - 0.0\n", + " - 0.0\n", + "alpha: 1.0\n" + ] + } + ], + "source": [ + "cat heat.yaml" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "x.initialize(\"heat.yaml\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Check the time information for the model." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{\n", + " \"time\": 0.0,\n", + " \"plate_surface__temperature\": {\n", + " \"value\": [\n", + " 0.7444227320828752,\n", + " 0.34022663250676577,\n", + " 0.8432555496724361,\n", + " 0.060846504441758764,\n", + " 0.8446697787107923,\n", + " 0.27627788779009277,\n", + " 0.5877718799413219,\n", + " 0.5737654859660629,\n", + " 0.34645112782286935,\n", + " 0.37414828806554334,\n", + " 0.7846696450611258,\n", + " 0.5304173216068899,\n", + " 0.0792535018455699,\n", + " 0.8793127450388358,\n", + " 0.15703439775879213,\n", + " 0.3912156458741114,\n", + " 0.9707058465142714,\n", + " 0.5511025163784199,\n", + " 0.1470283563347624,\n", + " 0.05645959012302082,\n", + " 0.9089090697986834,\n", + " 0.4495395294486181,\n", + " 0.011832496810315396,\n", + " 0.4171122020390934,\n", + " 0.33814771386922027,\n", + " 0.8466251708240787,\n", + " 0.6636025903978577,\n", + " 0.342616329035892,\n", + " 0.9154478512615142,\n", + " 0.14537348896585434,\n", + " 0.545772555177156,\n", + " 0.8827706254573504,\n", + " 0.2851195288529734,\n", + " 0.47448421167380395,\n", + " 0.2783120501807711,\n", + " 0.28066701069377975,\n", + " 0.841622847809508,\n", + " 0.7140017886819985,\n", + " 0.9744899247558164,\n", + " 0.45762952067167606,\n", + " 0.22664959884228209,\n", + " 0.7177555043888103,\n", + " 0.9900324558495176,\n", + " 0.7858216954853029,\n", + " 0.11380231880676006,\n", + " 0.16366335988935754,\n", + " 0.7884014952466945,\n", + " 0.14975089428180044\n", + " ],\n", + " \"type\": \"float64\",\n", + " \"itemsize\": 8,\n", + " \"nbytes\": 384\n", + " }\n", + "}\n" + ] + } + ], + "source": [ + "stateOut = x.get_state()\n", + "print(stateOut)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "stateOutDict = json.loads(stateOut)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'time': 0.0, 'plate_surface__temperature': {'value': [0.7444227320828752, 0.34022663250676577, 0.8432555496724361, 0.060846504441758764, 0.8446697787107923, 0.27627788779009277, 0.5877718799413219, 0.5737654859660629, 0.34645112782286935, 0.37414828806554334, 0.7846696450611258, 0.5304173216068899, 0.0792535018455699, 0.8793127450388358, 0.15703439775879213, 0.3912156458741114, 0.9707058465142714, 0.5511025163784199, 0.1470283563347624, 0.05645959012302082, 0.9089090697986834, 0.4495395294486181, 0.011832496810315396, 0.4171122020390934, 0.33814771386922027, 0.8466251708240787, 0.6636025903978577, 0.342616329035892, 0.9154478512615142, 0.14537348896585434, 0.545772555177156, 0.8827706254573504, 0.2851195288529734, 0.47448421167380395, 0.2783120501807711, 0.28066701069377975, 0.841622847809508, 0.7140017886819985, 0.9744899247558164, 0.45762952067167606, 0.22664959884228209, 0.7177555043888103, 0.9900324558495176, 0.7858216954853029, 0.11380231880676006, 0.16366335988935754, 0.7884014952466945, 0.14975089428180044], 'type': 'float64', 'itemsize': 8, 'nbytes': 384, 15: 0}}\n" + ] + } + ], + "source": [ + "stateOutDict['plate_surface__temperature'][15]=0\n", + "print(stateOutDict)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "x.set_state(json.dumps(stateOutDict))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Start time: 0.0\n", + "End time: 1.7976931348623157e+308\n", + "Current time: 0.0\n", + "Time step: 0.25\n", + "Time units: s\n" + ] + } + ], + "source": [ + "print(\"Start time:\", x.get_start_time())\n", + "print(\"End time:\", x.get_end_time())\n", + "print(\"Current time:\", x.get_current_time())\n", + "print(\"Time step:\", x.get_time_step())\n", + "print(\"Time units:\", x.get_time_units())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Show the input and output variables for the component (aside on [Standard Names](https://csdms.colorado.edu/wiki/CSDMS_Standard_Names)):" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "('plate_surface__temperature',)\n", + "('plate_surface__temperature',)\n" + ] + } + ], + "source": [ + "print(x.get_input_var_names())\n", + "print(x.get_output_var_names())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, get the identifier for the grid on which the temperature variable is defined:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Grid id: 0\n" + ] + } + ], + "source": [ + "grid_id = x.get_var_grid(\"plate_surface__temperature\")\n", + "print(\"Grid id:\", grid_id)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then get the grid attributes:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Grid type: uniform_rectilinear\n", + "Grid rank: 2\n", + "Grid shape: [6 8]\n", + "Grid spacing: [1. 1.]\n" + ] + } + ], + "source": [ + "print(\"Grid type:\", x.get_grid_type(grid_id))\n", + "\n", + "rank = x.get_grid_rank(grid_id)\n", + "print(\"Grid rank:\", rank)\n", + "\n", + "shape = np.ndarray(rank, dtype=int)\n", + "x.get_grid_shape(grid_id, shape)\n", + "print(\"Grid shape:\", shape)\n", + "\n", + "spacing = np.ndarray(rank, dtype=float)\n", + "x.get_grid_spacing(grid_id, spacing)\n", + "print(\"Grid spacing:\", spacing)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These commands are made somewhat un-Pythonic by the generic design of the BMI." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Through the model's BMI, zero out the initial temperature field, except for an impulse near the middle.\n", + "Note that *set_value* expects a one-dimensional array for input." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "temperature = np.zeros(shape)\n", + "temperature[3, 4] = 100.0\n", + "x.set_value(\"plate_surface__temperature\", temperature)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Check that the temperature field has been updated. Note that *get_value* expects a one-dimensional array to receive output." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 100. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0.]]\n" + ] + } + ], + "source": [ + "temperature_flat = np.empty_like(temperature).flatten()\n", + "x.get_value(\"plate_surface__temperature\", temperature_flat)\n", + "print(temperature_flat.reshape(shape))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now advance the model by a single time step:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "x.update()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "View the new state of the temperature field:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 12.5 0. 0. 0. ]\n", + " [ 0. 0. 0. 12.5 50. 12.5 0. 0. ]\n", + " [ 0. 0. 0. 0. 12.5 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. ]]\n" + ] + } + ], + "source": [ + "x.get_value(\"plate_surface__temperature\", temperature_flat)\n", + "print(temperature_flat.reshape(shape))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There's diffusion!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Generate an ensemble of models using get_state and set_state\n", + "We will demonstrate (some of) the use of get_state and set_state by generating an ensemble of models starting from the state of the one we just initialized, adding noise to each ensemble member and than running all of them forward in time. While running forward, we save the temperature at a single point of interest and plot a graph of all ensemble members to show the ensemble spread over time.\n", + "\n", + "First, we save the current state of the model." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "state_out = x.get_state()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "State is always returned as a string. In this case the string is formated as JSON, so Let's see what keys are in this state:" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dict_keys(['time', 'plate_surface__temperature'])\n" + ] + } + ], + "source": [ + "dictState = json.loads(state_out)\n", + "print(dictState.keys())" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dict_keys(['value', 'type', 'itemsize', 'nbytes'])\n" + ] + } + ], + "source": [ + "print(dictState['plate_surface__temperature'].keys())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we create an ensemble of BMI model objects. Each model gets initzialized in the same way as ```x``` was above." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "nEnsemble = 25\n", + "ensemble = []\n", + "\n", + "for ensembleMember in range(nEnsemble):\n", + " ensemble.append(BmiHeat())\n", + " ensemble[ensembleMember].initialize(\"heat.yaml\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, for each ensemble member, we take the state of x, add some noise, and set the state of the member" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "plate_surface_temp_from_x = np.array(dictState['plate_surface__temperature']['value']).reshape(shape)\n", + "\n", + "for ensembleMember in range(nEnsemble):\n", + " #for each ensemble, start with the same state as x.\n", + " ensembleMemberDictState = dictState\n", + " \n", + " noise = np.zeros_like(plate_surface_temp_from_x)\n", + " noise[2:5,3:6] = 2 * np.random.randn(3,3)\n", + " plate_surface_temp = plate_surface_temp_from_x + noise\n", + " \n", + " #this array we want to set in the state is a numpy array, which has to be transformed \n", + " #to a list otherwise we can not turn it into a json string.\n", + " ensembleMemberDictState['plate_surface__temperature']['value'] = plate_surface_temp.tolist()\n", + " \n", + " #Turn the state with noise added back into a json string\n", + " ensembleStateJSON = json.dumps(ensembleMemberDictState)\n", + " \n", + " #finally set the state of the ensemble member.\n", + " ensemble[ensembleMember].set_state(ensembleStateJSON)\n", + " \n", + "#Note that the above loop is needlessly verbose for educational purposes.\n", + "#This does use a lot of memmory. The code below is less readable, but uses less memmory.\n", + "#When using bigger models than this example model, this can make a big difference in \n", + "#performance.\n", + "\n", + "\n", + "plate_surface_temp_from_x = dictState['plate_surface__temperature']['value']\n", + "noise = np.zeros(shape)\n", + " \n", + "for ensembleMember in range(nEnsemble):\n", + "\n", + " noise[2:5,3:6] = 2 * np.random.randn(3,3)\n", + " \n", + "\n", + " dictState['plate_surface__temperature']['value'] = (plate_surface_temp_from_x + \n", + " noise).flatten().tolist()\n", + " ensemble[ensembleMember].set_state(json.dumps(dictState))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now run the entire ensemble forward in time. Every timestep, we are saving the temperature of one location of the plate of interest, for each ensemble member.\n", + "\n", + "Note that for bigger models, the for loop in this step can be run in parallel for all the models. " + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "loc_of_interest = [14]\n", + "\n", + "distant_time = 15.0\n", + "output = pd.DataFrame(columns = ['x'])\n", + "\n", + "outputValue = []\n", + "while x.get_current_time() < distant_time:\n", + " x.update()\n", + " x.get_value_at_indices('plate_surface__temperature',outputValue, loc_of_interest)\n", + " output.loc[x.get_current_time(),'x'] = outputValue[0]\n", + " \n", + " for ensembleMember in range(nEnsemble):\n", + " ensemble[ensembleMember].update()\n", + " ensemble[ensembleMember].get_value_at_indices('plate_surface__temperature',outputValue, loc_of_interest)\n", + " output.loc[ensemble[ensembleMember].get_current_time(),'ensemble' + str(ensembleMember)] = outputValue[0]\n", + " \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "View the final state of the temperatures of the entire ensemble at the locaiton of interest:" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'temperature')" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd4VEX7/u85W7JppCekQiAQkBBKQm+hhl4iSPlKr6+iwssPQQSVJoIgFhBQiqCoICiCdEW6dCnSQuhNAoEUEkKS3fv3x+4eUijxNdkEmM91zcXu2XPmzJzo3GfmeeZ5BElIJBKJRAIASlE3QCKRSCTFBykKEolEIlGRoiCRSCQSFSkKEolEIlGRoiCRSCQSFSkKEolEIlGRoiCRSCQSFSkKEolEIlGRoiCRSCQSFW1RN+Cf4unpydKlSxd1MyQSieSp4uDBg7dIej3pvKdOFEqXLo0DBw4UdTMkEonkqUIIcTE/58nlI4lEIpGoSFGQSCQSiYoUBYlEIpGoSFGQSCQSiYoUBYlEIpGoSFGQSCQSiYoUBYlEIpGoPHX7FJ4Hdu3ahf379yMsLAxVq1aFp6dnUTdJIpE8J0hRKEYkJSXhzTffxBdffJHjeEBAAKpWrYoaNWrgjTfegIuLSxG1UCKRPOvI5aNiwurVq/HCCy9g/vz5GDFiBC5duoRff/0V06dPR1RUFC5cuIDx48ejZs2aOHHiRFE3VyKRPKPImUIRc+PGDbz++utYvnw5wsPDsWrVKtSoUQMAEBgYiKZNm6rnbt++HV26dEGtWrXw1Vdf4cUXXyyqZkskkmcUOVMoQpKTkxEZGYlVq1Zh0qRJOHDggCoID6Nhw4Y4dOgQwsLC0LlzZ4waNQpZWVk2bLFEInnWkaJQhEyePBlXrlzBli1b8Pbbb0On0z3xGn9/f2zduhVDhgzBtGnT0LJlS9y6dcsGrZVIJM8DUhSKiLNnz+Ljjz9G7969Ua9evX90rZ2dHebMmYOFCxdi586daNWqFdLS0gqppRKJ5Hmi0ERBCLFQCBEvhPjrEb//nxDiqKXsFkJUKay2FEdGjhwJnU6H999//3+uo2/fvli+fDkOHjyIfv36gWQBtlAikTyPFOZM4SsALR/z+3kAjUiGA5gI4IvHnPtM8fvvv+Onn37CmDFj4Ofn96/qat++PaZMmYJly5Zh4sSJBdRCiUTyvCIK8+1SCFEawC8kw55wnhuAv0j6P6nOyMhIPs1JdoxGIyIiIpCYmIiTJ0/C3t7+X9dJEn369MGSJUuwfPlydOnSpQBaKpFIniWEEAdJRj7pvOLiktofwPqiboQtWLhwIY4cOYLly5cXiCAAgBACX3zxBeLi4tC7d2+UKVMGERERBVK3RCJ5vijymYIQojGAzwHUJ5nwiHMGARgEAEFBQREXL+Yrq1yxIykpCeXKlUOFChWwbds2CCEKtP74+HjUqFEDWVlZ2L9//79empJIJM8O+Z0pFKn3kRAiHMB8AB0eJQgAQPILkpEkI728nph3utgyefJk3Lp1CzNnzixwQQAAb29vrFmzBklJSYiJiUFmZmaB30MikTzbFJkoCCGCAPwIoCfJ2KJqh604f/48Pv74Y/Tp06dQl3bCw8OxaNEi7N27F1OmTCm0+0gkkmeTwnRJ/Q7AHwBChRBXhBD9hRBDhBBDLKe8A8ADwOdCiMNCiKfXepwPvvnmG2RlZWHChAmFfq8uXbqgR48emDhxIg4ePFjo95NIJM8OhWpTKAyeVu+jatWqwdHRETt37nzseQcOHMD777+PuLg4REVFoXnz5oiKioKzs/M/ut+dO3cQFhYGV1dXHDx4EAaD4d80XyKRPOU8FTaF54Vz587h8OHDiImJeeQ5u3btQqtWrVCjRg1s3boV3t7emD9/Ptq3bw93d3c0aNAA77//PlJTU/N1Tzc3NyxYsAAnTpzAuHHjCqorEonkGUeKgg346aefAACdOnXK89vOnTvRuHFj1K9fHwcPHsQHH3yACxcu4Ndff8WdO3ewZcsWjBw5Eunp6Xj77bdRt25dnD9/Pl/3bdmyJQYPHowZM2Zgx44dBdoniUTybCKXj2xA3bp1kZ6ejkOHDuU4fvToUURGRsLT0xNvvvkmBg4cCEdHx0fWs3HjRnTr1g2KomD58uU5wmo/irt376JKFXMEkSNHjsDJyenfdUYikTyVyOWjYsK1a9fwxx9/5Ml9kJmZiT59+sDNzQ1Hjx7FsGHD4OjoiPj4eIwZMwZdu3bFrFmzcPz4cTWmUXR0NPbv3w9fX19ER0fj448/fmK8IycnJ3z11Vc4f/48Ro4cWWj9lEgkzwgkn6oSERHBp4nZs2cTAE+cOJHj+IQJEwiAK1euJEleunSJr732Gg0GAwHQ29ubAAiABoOB3t7e9Pf355QpU3jt2jV26tSJANizZ0+mpaU9sR3/7//9PwLgr7/+Wij9lEgkxRsAB5iPMVYuHxUyzZo1w9WrV3Hy5En1mHXZ6MUXX8TkyZMxZcoULF68GCTRrVs3nDp1Ctn7aDAYYDKZkJGRAcBsRB4xYgTu3buHyZMno2XLllizZg202kdHLUlPT0flypUhhMDRo0elN5JE8pyR3+WjIn/z/6flaZop3Lp1ixqNhmPGjFGPZWRksFq1avT29uaBAwfo7OxMOzs7vvrqqzxx4gQbNmxIRVE4c+ZM/vnnn0xOTiZJmkwmfvvtt7Szs6NerycAurm5sX379gTAIUOG0GQyPbY9GzduJACOHz++UPstkUiKH8jnTKHIB/l/Wp4mUVi4cCEB8MCBA+qx8ePHEwBXrFjBRo0a0dnZmWfOnGFqaiobN25MRVH47bffPrLOv/76i+XLl6eiKKxUqRIB0NPTkwA4Y8aMJ7apW7dutLOzY2xsbIH0USKRPB1IUSgGtG3blkFBQeob/OHDh6nVatm9e3fOnDmTALhgwQKmpqaySZMmFELw66+/fmK9SUlJqk2hYcOGdHBwoKOjI4UQ/Omnnx577bVr11iiRAk2b978iTMLiUTy7CBFoYhJTk6mXq/nsGHDSOZcNtq9ezcNBgPbtm3L1NRUNmvWjEIILl68ON/1m0wmTp06lYqisEqVKnRxcaFOp6PBYOD+/fsfe+2sWbMIgN99992/6qNEInl6kKJQxHz//fcEwB07dpAkFy9eTABctmwZa9SoQXd3d167do1t27alEIKLFi36n+7z3XffEQBbtGhBLy8vKopCDw8PXrx48ZHXZGVlMTIykiVLlmRiYuL/dF+JRPJ0kV9RkPsUCokff/wRPj4+qFOnDgBzQLzSpUvj1KlT2L9/P+bMmYOTJ0/il19+wdSpU9GnT5//6T7dunXDhx9+iE2bNqFDhw7w9fVFQkICmjRpgvT09Ideo9FoMHfuXMTHx2Ps2LH/axclEsmzSH6UoziVp2GmkJaWRkdHRw4ePJgkefXqVSqKwv79+1Or1bJbt24kyRYtWtDHx4f37t37V/czmUwcOnSo6lkUGBhIAOzevftjr3vttdcohOC+ffv+1f0lEknxB3L5qOj4+eefCYAbN24kSc6YMYMAGBISQl9fXyYkJPDQoUMEwClTphTIPbOystixY0cKITh//nw6OjoSwGM9mRITE+nr68saNWrQaDQWSDskEknxRIpCEdK/f3+6uLjw/v37JMlq1aqxdOnSBMA1a9aQNLuGOjs7886dO/mqMzMzk5cuXXqsx1BaWhrr1KlDg8HAOXPmUAhBnU7HS5cuPfKar7/+mgA4f/78f9BDiUTytCFFoQipWLEi27RpQ9K8rwAAy5Qpw9DQUJpMJp49e5aKonDkyJGPrOPGjRtctWoVR48ezUaNGtHBwYEAGB4ezlmzZj3SQHzz5k2WK1eOnp6efOuttwiAgYGBj5wJmEwm1qtXj15eXvkWKIlE8vQhRaGISExMpBCCEyZMIEmOHj2aGo2GAPj++++TJF955RXq9XpevXo1z/X3799nlSpV1LhHsMQ+cnV1pZ+fH4OCggiA9vb27Nu3L/fs2ZNn9nDy5Ek6ODiwcePGbNq0KQGwc+fOj2zzoUOHKITg66+/XoBPQiKRFCekKBQRmzdvVu0JRqORgYGBDAkJoRCCly9f5o0bN2gwGNi/f/8818bHx9Pf3199u2/cuDFjYmLYtWtXdu3aldWrVycAhoaGsm3btqrdoFOnTkxJSclR16JFiwiA77zzjhpc73FLREOGDKFGo+GxY8cK/JlIJJKiR4pCETFp0iQC4J07d7h161YCoJeXF5s1a0aSHDt2LIUQPHXqVI7rDh48SBcXFwJgdHT0Q20H1vhHVu+itm3bcsSIEeoGtux7E0wmE19++WUqisKlS5dSo9FQo9Hkua+VW7du0d3dnVFRUXKns0TyDCJFoYho164dQ0NDSZIDBgygvb09AXDJkiVMSUmhm5sbO3XqlOOaJUuWUKvVEgDbt2/PtLQ0bt++nd988w0nT57MQYMGMTo6mrVr1+bMmTN548YNTpo0iY6OjtTpdBw4cCBLlChBHx8f7t69W603OTmZ5cqVo7+/Pz/55BMCoL+/P7Oysh7a9jlz5qgb7CQSybNFkYsCgIUA4gH89YjfBYBPAcQBOAqgen7qLc6iYDKZ6O3tzV69evHevXt0cXFhSEgInZycePfuXX700UcEwD179qjXvPnmm6rtoG3btrx8+XIem4KXlxcbVq3K/ytXjjqArq6uHD16NP/880927NiRAPjKK6+wbNmy1Ov1OeInHTp0iHq9nm3btmVUVBQBqPsncpOVlcVq1aoxICCAd+/eLfTnJZFIbEdxEIWGAKo/RhRaA1hvEYfaAPbmp97iLArnz58nAM6ePZsrVqxQDcJ9+vTh/fv36e/vz6ioKPX8bdu2EQCFEGzevDmPHz/O4OBgOjg4sEe7dpxQuzaXlCzJQxoNswAS4G0huMLTk1FCUK/Vsk+fPuzcubMqDNaBf+zYseoy0KeffkoAnDRpkjpz2b59+0P7sHPnTgLg22+/bZNnJpFIbEORi4K5DSj9GFGYB6B7tu+nAfg+qc7iLArWeEf79+9nx44dVRvB77//zpUrVxIA165dS9I8q4iIiKAQgjVr1uSuXbvo7e3NEs7O/EirZaZFBO4D3AZwkqJwSIkS/AbgXctv1zQaztBoGO7mpgrDwIED2a9fPwLg1KlT1Xt17NiROp2OU6dOJSyzDWuuhty8/PLL1Ov1jIuLs9mzk0gkhcvTIAq/AKif7ftvACKfVGdxFoXhw4fTzs6O169fp06nY6lSpViqVCkajUZ1Q1tmZiZJ8scff1SXhxYtWkRnZ2d6enpyohAkwG8Uhf+tUoVTxo3j9u3b+e2337J9+/Z88cUXWbNSJb6s0XANwAyAtwC21WrVhDsvv/wyu3btqtoySDIhIYG+vr4MCwtj8+bN1eWqh3H16lU6OjqyQ4cONnt2EomkcHkaRGHtQ0Qh4hHnDgJwAMCBoKCgwnliBUC9evVYp04dzp07V10WGjduHE0mE/39/dW9ApmZmSxXrhw1Gg2rV69OvV5PX19fvm6ZAXxtZ8dTJ0/SZDJx3bp1rFatGgEwICCAfn5+OewNZQEeBmgEOA5gsyZNCIAxMTFs3LgxtVqtGm7jl19+IQAOHz6cTk5OBMBvvvnmoX2ZMmVKjlAdEonk6eZpEIVnavkoIyODBoOBw4YNY+vWrdVsaGfOnOGxY8dy7BP44osv1EFdp9MxMDCQPS2CsEqj4dLFi9mtWzd6eHjkEAAAdHd356hRo3j06FFu2LCB06ZNo4tOx8WW69cCjAoPJwAOGzaM4eHhdHJy4sGDB0mSffv2paIo6qBvZ2fHCxcu5OlPeno6y5Yty4oVKzIjI8Omz1IikRQ8T4MotMllaN6XnzqLqygcPHiQALh06VK6uLjQ1dWVdevWJUl++OGHBMDLly/z7t27LFmyJA0Gg5r/oB3ATIC/CsFx2byRnlSCgoI4bdo0njt3jnVq1+Zgiw3iHMDmFkH5+OOPWapUKXp7e/Ps2bNMTExkQEAAK1SowFatWhHAIwPirVmzhgA4c+ZMWz9OiURSwBS5KAD4DsB1AJkArgDoD2AIgCGW3wWA2QDOAjiWH3sCi7EofP755wTAdevWqYP2vHnzSJJNmjRh5cqVSZKTJ09Wf9doNGwEMB3gHoC9LCk2/2mxs7Pj2rVrOXv2bNYEeBFgCsDmzs7U6XT89ttv6e7uzpCQEN66dYsbN25UvZWsy0iffPJJnj6ZTCa2bNmSJUqU4I0bN2z6PCUSScFS5KJQWKW4ikLv3r3p5eXF2bNnEwD1ej3v3LnDlJQU6nQ6jhw5kjdv3qSzszPd3Nzo7OxMFyF4C+BfAJtZ7AYPWy6ybmyrXbs2K1WqRIPB8NBze/XqxSNHjjDE0ZGxAG8DrGWZkaxcuZJ6vZ7R0dHMysrioEGDKITg22+/rS5jPczb6NSpU9RqtRwwYEARPFWJRFJQSFGwMRUqVGDbtm3Zo0cParVatmjRgiS5evVqAuBvv/3G4cOHUwihGqFHWOwArby88gzwoaGh1Ol06rnW4w4ODuzXrx83bNjAt99+m25ubjmu8/T05ObNmxnm7MxLAP8GWEmrZXh4uLqr+b333mNycjJLlSrFsmXLsm7dugTAWrVqPXQZacSIERRCPDH3s0QiKb5IUbAhd+7cIQBOnDiRAQEBBHJGRHV0dGRsbCz1ej2DgoKo1+tpEIJXAf4KUFGUHAO7VQQURVGFIXuxnh8cHMzx48eru5qzl/79+zNMr+ffAC8BDALYrl079urVi0IIrl+/nlu2bCEA9u7dW73Pw5aRkpKS6OPjwzp16si4SBLJU4oUBRuyadMm1b3TOijv3LmTJpOJwcHBbNeunertYy19LbOE5rkG8+yCEBYWxp49e3LGjBn87bffePz4cU6aNInlypXLca6npyfHjBmjhui2liZNmrC6RsPbAM8ALAnw3XffZXh4ON3d3Xn+/Hm+8sorFELwP//5j7rsFRsbm6ePCxcuJPBg34NEInm6kKJgQyZOnEgA/PLLL9WBNT09nadPnyYAfv7556xZsyY9PT2pKAo1QvC0EDz4CEFwdHTM4yZ6/fp1zpkzh8uXL+fevXu5ceNGvv7663R1dVWvHTlyJH19fXPUWb16dda2GJ6PAnSFOUVniRIlGBERwRs3bjAwMJAVKlRghQoVKIRgrVq18gTNMxqNrFWrFkuWLMmkpCRbPl6JRFIASFGwIW3btmVoaChfeeUVKorChg0bkqS6hr9r164cyz4dLLOErg8xFuv1ep49e1at+9q1a3zjjTcealy2t7dn+fLlWbZsWfVYq1at1B3L1hISEsLGMLurbgDo5e7OJUuWEDAHx7N6TPXt21cVpo8++ihPP/ft22e2hYwYYbNnK5FICgYpCjbCZDLRy8uLvXr1YqVKlQjLEg1JtmrViuXLl1c9kgBQEYJ7NRqeBah5yExh27ZtJM2hJl5//XUaDAZqNBr26dOHx44d45EjR/jzzz/zk08+4X//+19GR0cTAJ2cnFTR8ff3Z8OGDXPU7ePjwwEWMZoKMCoqSo3QunjxYvbs2ZNarZbdu3d/7DLSgAEDqNVqeeLECVs+ZolE8i+RomAjrJFRP/zwQ/Ut+7fffmNaWhrt7e35+uuvs1mzZuqbfgPLwPyfh8wSZsyYQdK858HOzo4ajYb9+vV7YmC6P/74g7Vq1SIA1X1Vr9ezZs2aOep3dHTkLMv9ewAcN24cGzZsSEdHR+7bt49eXl6sWrUqfX19qSgK69evn8cbKT4+nq6urmzWrJk0OkskTxFSFGyENTKqdalIURSmpqZyw4YNBMwJazQajfoWv0mn4w2AhlwzhPbt29NkMnHHjh3UaDRs0aKFuox07tw5Dh8+nKVLl2atWrXYo0cPjh07losWLeK2bduYlpZGk8nEpUuXquk8rcJQLdf+BwedjlsBpgGMALh8+XK6ubkxMjKSS5cuVb2RrOfPnj07T59nzZpFAFyxYoWtH7dEIvkfkaJgI4YPH06DwcBRo0apYbBJctiwYTQYDJw/f746wIZZ3tLH5poh+Pv78/bt24yPj6efnx9DQkKYmJjIHTt2sGPHjhRCqLMQnU5Hg8GQY++Cn58fFy1axKysLKampnLcuHHqb3Z2dgwLC8txP28heAFmV9VAvV71LHrzzTfZsWNHGgwGtmjRgoqi0MHBIY/ROzMzk+Hh4QwKCmJqampRPHaJRPIPkaJgI+rWrcu6deuyfv36BMDRo0eTNG9mi46OZkxMjJrY5mshmALQLZcorFixgkajkdHR0bSzs+PChQsZbglqZx38/fz8OGrUKA4aNIj169fPsWnNuscgPDycGzZsoMlk4o8//qjOTuzs7BgaGprjnlUBpgLcATA0OJgDBw6kEILff/89XVxcWKtWLTo5OVGj0bB58+Z5loq2b99OWJagJBJJ8UeKgg3IysqiwWDg0KFDqdfrCZhjH1ntDFOnTqWDgwMNBgMDYc598FEuQfD19WVmZiYnTZpEwLzb2M7OTv29RYsWXLNmTR4XUZPJxBs3bnDt2rWMiIhQB38AbNasGQ8fPsw//vhDbZfBYKC3t3eOe79kmbnMhjlERoUKFejn58eZM2cSALt166aeu3Dhwjz979GjB+3s7GQyHonkKUCKgg04c+YMAfCtt95S3+oTExM5Z84cAlDX3gFwuGUADs4lCuPHj+eWLVuoKApffPFFellCXrRv357nzp3LVzuMRiMXL17MkiVLqrYEnU7HRYsW8cyZM3R0dFRFI7vgAOA0S7tegtkNVa/Xs3379mzcuDGdnJxYvXp1arValihRglevXs1x36tXr9LJyYlt2rSRRmeJpJiTX1FQIPmfOXHiBAAgOTkZABAWFgYXFxds2bIFgYGB2L9/PwwGAwBznPC/AJy3XCuEgKIo6NChA3r06IFy5crhzJkzuHnzJmrWrImVK1ciODhYvdft27cxffp0jB49Gq+++ip69+6NmJgYREdHY9y4cYiKisKZM2cwduxYCCFAEn379sW8efMQFxcHd3d33L9/H1qtFkIItd4xAHYDmA/gy5EjMWbMGKxevRqNGzeG0WiEk5MTAODu3bt45ZVXzG8SFvz8/PDee+9h7dq1WL16dSE9ZYlEYlPyoxzFqRSnmYI1dEV0dDSFEBw+fDhJskyZMoyJiaG7u7s5IirMG8c+yPaGbmdnxxdffJEtWrSgwWBgy5YtVdtBYmJijvusX79e3ams1+vp7u7OUqVKsVKlSqxevToVRaGiKGzfvj03bNjAM2fOMCwsTLVHtG3bltevX1dzRlsTAFlLAMwpPf8EGOLvz+joaBoMBnUfQ4cOHdRzv/vuuxxty8jIYOXKlRkUFMS7d+/a7NlLJJJ/BuTyUeHTq1cv+vr6qoPtTz/9xNu3bxMwB6SDxUW1o2WJpmGupaPPPvtMHbQB8w7l7EtGd+/e5ZAhQwiAFSpU4IwZM/j7778zNjY2xwB84cIFjhkzRl16CgkJ4YwZM9SNbUIIVqpUiXv37lVtDNldVwGwlaWNc2EOnOfp6clq1aoxMjKSnp6eLF26NO3s7Ojh4cH4+Pgcz2HHjh0EwFGjRtns2Uskkn+GFAUbEBkZyTp16qgD682bN/nbb7+pb9dWr6AvASYC1OJBjoTQ0FD26NFD9UxSFIU7duxQ6969ezdDQkIIgHXq1KG7u3uOQRwAXVxcWLFiRTVXQ3p6Or/99lu1TdHR0ezXr5/qoeTh4cEVK1aoXklWG4S1TLEIQ3eYw19YxU2r1bJp06ZqO7t27ZrnWfTt25darZbHjx+35Z9AIpHkEykKhYzJZKKjoyObNGmivp2T5LRp09Q3ceugexXg8lwD+oQJE6jVatXIpgsWLFDrnjt3LoUQdHZ2VndCt2jRgiNGjOD777/P6dOnc/LkyXzttdfYunVrCiHo5OTEsWPH8s6dOzSZTOquaF9fXw4dOpRCCOr1erq6uubIEe3s7Kx+1sLsopoCMBRQ90hYN7M1a9ZMFZQff/wxx/O4efMm3d3d2ahRI2l0lkiKIVIUCpmLFy8SAGvWrKmGnibJbt26qWJgMBhYxfL23RsPYhDZ29tzzJgx6mD88ssvq/WePn2aGo2GQgh1Z3OLFi3U8BXWotfrWbFiRXbo0IFTp05l586dCYCurq6cNGkSk5OTeeTIETXyabdu3WhnZ0edTkd3d3e+8847al3Z8zn4A4wHeARgoKcng4ODGRQUxNDQUPr5+dHLy4sODg709vZmQkJCjmdiFZuvv/7apn8LiUTyZKQoFDLr16/PYbT9/vvvSZLlypVj+fLl1UH2LYso+GQzMPfp04d+fn7mMNoaDe/cuUPS7FpasWJFAuALL7yg5k1wc3PjiBEjeOTIEe7YsYMLFixQdx8HBwcTMIfInjdvHtu3b6/OVPbs2cO7d++yT58+BMCwsDDVXdXT01M9nru0tLR5NsC6detSURS2atWKQgg1AquiKOzVq1eOZ2INr+3t7c3bt2/b/G8ikUgejRSFQuajjz7KMZBevXqVSUlJBMCgoCDVW2gnwP25Bt0PPvhA/RwTE6PWac2XbN1LUK1aNS5YsOCxoSSMRiO/+eYbBgUFEQBbt27Nb775hsHBwdTr9Zw/fz5J8ptvvqG9vT39/f2p1Wqp0+no7e2t2h9yJ+j50CIMnSxLV9a6YREKq73kl19+ydGeQ4cOUVEUdeYkkUiKB8VCFAC0BHAaQByA0Q/5PQjA7wD+BHAUQOsn1VlcRGHAgAGq15Gvry9JcuvWrerSjrOzM90BZgF8L9ssoUaNGmzYsKG6ZHPp0iWS5JEjRyiEUI9//fXXedbmTSYTY2NjuXTpUg4bNoz169eng4MDg4OD2bNnT3bt2pUlSpSgoijs16+fahweMmQI79+/z+3bt9PZ2ZleXl5qqs+SJUuqXktWzyQA1AHcB/A2zKk8K1asSDc3N5YqVYpBQUF0cnKik5PTQ11o33jjDQoh+Mcff9jmjyGRSJ5IkYsCAA2AswDKANADOALghVznfAHgP5bPLwC48KR6i4so1K1bV/UOeumll0jmnT10s7xt18xlYLZ+btCgAUlzgDkPDw/1+KQfVyMTAAAgAElEQVRJk/Lcb9myZTn2FxgMBtatW5evvvoqO3XqlCMWklWsIiIi+Oqrr6pv99euXeP+/fvp4eFBV1dXCiGo0+kYGBiovvlnnzGUBZhkMT57u7vTwcGB1atXJwA2btyYgNndtX///jnampyczMDAQIaFhfH+/fuF/8eQSCRPpDiIQh0AG7N9fwvAW7nOmQdgVLbzdz+p3uIgCiaTiW5ubuqa/5w5c0iS//d//6cOyAC4xGK0VSzfNRoN+/btq/5+5MgRklSXZQCwZ8+eOWYIRqNRjXpau3Ztfvnllzx8+DAzMzNztMloNPLPP//kzJkz2a5dO3XG4eLiwkmTJtHBwYG+vr48fPgw//rrL5YsWZJOTk5qu3IHzLOW7hZhm2CxcwBgvXr1CICVK1dWvaPWrVuXoz1r1qwhAE6cOLGQ/xoSiSQ/FAdR6AxgfrbvPQHMynWOL4BjAK4AuAMg4kn1FgdRuH79uupJBIC7d+8mSVasWJEBAQHUarVUAN60CANgTn5Tu3Zt1V5QqVIlkuaEOtYBuF69ekxPT1fvc/fuXcbExBAwp8rM/tuT+Ouvv9TEO0IIDhgwgAEBAXR3d+fBgwd55swZlipVSh3UhRCPFIaFAI0Ao2DeRKfX6+nr68tSpUpRr9ezRIkS9PPzy2Ncfumll6jX63nq1KkCeOoSieTfUBxEoctDROGzXOf8F8AIPpgpnACgPKSuQQAOADgQFBRUaA8tv1g3qFnX4FNSUpiSkkIhBD09PakoCmtZ3rCz52Fu1aqV+nnz5s1MSEhQl2sCAgJy7BS+cOECq1SpQkVR+NFHH/1Pvv8mk4lfffWVOvCXLl2a/v7+dHV15b59+3jp0iWWK1cuhy0hMDAwjyg4AjwJ834LT5i9oaxeTw0aNCBg9kbq2bNnjvtfv36drq6ubNiwYZ4MbhKJxLYUB1HIz/LRcQCB2b6fA+D9uHqLw0whe/RTLy8vkuTOnTtzDKTjLUbm7LkTrDOLgIAAkmaDLGDebZw95/G+ffvo7e3NEiVKcP369f+6vUlJSaxRo4Zqi/Dz82OJEiW4e/duXrx4kQEBATmip5YoUUKdPViPVQGYDnANQDdXV3X5yDrDsEZiXbVqVY57W5MMffHFF/+6HxKJ5H+nQEUBgD2A0Pycm+0arWWQD8YDQ3OlXOesB9DH8rkigGsAxOPqLQ6i8Morr9DBwYEA2LBhQ5JU03Fay36LgVZ927YMmgC4ePFi3rt3Tx2IZ82apdadmprK0qVLs1SpUjx58mSBtdlkMnHAgAHqYO/t7U0nJydu376dp06dooeHB/V6vWqLyO2iCoBDLbOfYRZhA8z7NIKDg6nT6ejm5kYfHx/eunUrx32joqLo4uLCa9euFVh/JBLJP6PARAFAO5jdSs9bvlcFsDpflQOtAcTC7IX0tuXYBADtLZ9fALDLIhiHAbR4Up3FQRQaN26sDorWIHC9e/dWB/6SlsHzLTzIjFaqVCkC5h3HRqORc+fOJWB2U82eQMeam+H3338v8HabTCaOGDFCHeRdXV3p4ODArVu38uDBg3RycqJOpzPbRLLtcs5eVsEc8bW6RRCsHlF169ZVxSR3bKTY2Fja2dmxc+fOBd4niUSSPwpSFA4CcAHwZ7ZjR/NTeWGU4iAKPj4+apRR63JJWFgYvby8qNFo2MciCuHZBlPrzGLy5Mk0Go3q3oCBAweq9R4/fpw6nS7PTuGCJnsOZwcHBzo5OfHgwYPcunUr9Xq9urktu63BWtxhzu0cC9DJIgJlypQhYI7/ZPVoWr58eY57Tp48mYA5kqxEIrE9BSkKey3/SlEgmZCQQODBXoDz588zNTWVQgg6ODhQCMFvLUbZh71pp6WlcfXq1er3pKQkkua3+EaNGtHNzY03btwo9H6899576lKSVqulh4cHT58+zTVr1lBRFPV47kxtANjAYi9ZYgnEZ11G8vPzo1arpbu7Oz08PHj9+nX1fhkZGaxWrRq9vb158+bNQu+fRCLJSUGKwgIAPWDecVwOwGcA5uan8sIoRS0KVoOyoijU6/U0mUz8448/cgyacQB/sHwWQqhhrwMDA0mS4eHh6pKLlcWLFxMA582bZ7O+jB8/ngDU5SJ/f39euXKFX3/9tdoXa3ym3MLwjmU21AtmbySrh5M1X7ROp2OrVq1yeE0dOXKEOp2O3bp1s1kfJRKJmYIUBQcAkwHst5RJAAz5qbwwSlGLQvaw0+XKlSNJzp49O8fyCgG+iQcb1qyiMHLkSO7Zs0c917q/ISEhgZ6enqxdu7bNXTetS0l2dnYUQrBMmTJMSEhQl3seVRSAW2EOs10O5gRB3t7eBMDg4GDVg+nzzz/Pcb+JEycSAFesWGHTfkokzzsFIgowh6r4MD8V2aoUtSgMGzZMXWu3hrfo168fDQYDhRBsYRGFxtkGUKtr56VLl9Tdyz4+PmqdAwcOpEaj4eHDh23eH5PJpGZ3s9o9wsLCmJKSomaPe1TxhzmN5yGAdlZRdHent7c3dTodfXx8aDAYcmxey8jIYEREBL28vPJkcJNIJIVHQc4UtuSnIluVohaF6Oho1Uj82WefkSSrVq1KFxcXKorCty2iUCKXIDg7OzMuLk79PnXqVJLkrl27CID//e9/i6xPWVlZaj4G6xt+zZo1mZqaqkZIfZQ3UhtLf+doNOoSkqIorFatGgGzK25kZCQzMjLU+x07dox6vZ5dunQpsj5LJM8bBSkKMwCshnlHcoy15KfywihFLQqBgYGqG+aePXt47949ajQaddBcBfPuX+u6unXpKCYmRk2NqdFomJycTJPJxKpVq9Lf35/JyclF2q/09HQ2bdqUGo1GDa7XqlUrJiUlMSwsTLU9ZN/QZi3TLMLQTVFob2+vGp9DQ0PVdKPjxo3Lcb/333+fALhs2bIi6rFE8nxRkKKw6CFlYX4qL4xSlKKQnJxMWNbPAfDu3bvcv39/jgHyKh7EOwIehMJYu3atmj2tQ4cOJB8kvLfmPChqkpOTGRkZSYPBoHpX9evXj1evXlWzyeXOAAeY03juhjmiahmL6Lm6utLJyYmOjo4sWbIkhRCqDYU0R4atUaMGPTw8+PfffxdhryWS54MCE4XiVopSFPbt2/fAoOzuTpLqJjQA9LO8Mb+We9DUajlz5kz1+86dO0mSffr0oZOTE1NSUoqsT7mJj49n+fLlWaJECXUz3rhx43js2DFVDB+2lBQEMAHgAYD6bIJozULn5ubGMmXK5JgRHT9+nHq9nu3bt5d5nSWSQia/oqDgCQghFgkhFuYuT7ruWeTEiRPq5/LlywMADh06BL1eD0VRUMPy237LvzqdDgBQuXJlfPXVVwCAgIAA1K1bFykpKVi+fDm6desGJycnG/XgyXh5eWHjxo1wcHCAi4sLDAYDJk6ciG3btmHVqlUQQsBkMuW57hKAPgAiAMzUaqHVapGVlYXY2FhUrFgRqampOH/+PIYOHape88ILL+CDDz7A6tWr8eWXX9qqixKJ5DE8URQA/AJgraX8BqAEgLuF2ajiyokTJ6DRaAAANWvWBAAcPHgQGo0GJFEDQCbM8Tq0Wi30ej0AYPDgwTh27BgA4NVXX4UQAsuWLUNaWhr69+9fBD15PKVLl8batWuRlJSEUqVKQafTYejQoUhKSsKcOXMeed0aAB8BeCUrC50VBVqtFg4ODjh9+jQcHR3h7e2NJUuWYMmSJeo1b7zxBpo3b45hw4bh1KlThd85iUTyePIzncheYBaSIvNIKsrlo3bt2qneOT/99BMzMzPVjGUAuNHinmn9bjXKLl26VD1mTb9Zu3ZtVqxYsVgvm6xbt44ajYa1atWioihUFIVbtmzhsGHDcvQve9EB3AMwEebMbYrF+Gy1SQQEBNDR0ZGnT59W73P16lV6eHiwevXqMlObRFJIoLBsCgBCAcT90+sKqhSlKJQtW5bOzs4EwIsXLzI2NjbHgJgAcF6uQdLX15fR0dHqZ9K8lg6A06dPL7K+5Jcvv/ySANi8eXM19MXhw4fV/RYP2+0cBPP+hSMA7bMJZMWKFVX33KpVq+ZIGvTTTz8ReBBgUCKRFCwFJgoAUgAkZyuxAF7MT+WFUYpKFNLS0tRk9zqdjiaTKUcMo7IWI/OAbO6oADh8+HDVQGsNfjdixAhqtVqbxDgqCMaOHUsAbNeunep9dfLkSVaqVOmRwhANc7a2pVotRTbjdFBQkDrbev3113PcZ9CgQRRCFEqEWInkeafQZgpFXYpKFA4fPqwOeMHBwSTJadOmqd443SyiUAUPYgYByCEcK1eu5P379+nl5cVOnToVST/+F0wmE3v16kUAbNOmjepNdPz4cXUj38OWkqzxkd4wGKjRaKjX6+ng4EB7e3s19Hj2pDx3795l+fLlGRAQkCe1p0Qi+XfkVxTy4330W36OPetkN4JWrlwZAHDy5EnVG6cmgHswp5IDAJPJBAcHB/z888/qdY0bN8batWtx8+bNYmlgfhRCCHz55Zdo1qwZNm7ciCZNmuDOnTuIjo7G6tWrYWdnZ51V5mAizFmUpqanI5JEZmYmMjMz4eHhgStXrsDX1xd9+/bF5cuXAQCOjo749ttv8ffff2Pw4MEPrVMikRQujxQFIYRBCOEOwFMI4SaEcLeU0gD8bNXA4kJcXJz6uX79+gCAo0ePqgNXDQB/AsjKdk2DBg2wZs0aAEB4eDjc3NywYMEC+Pn5ITo62kYtLxj0ej1WrFiBChUqYN++fahevTquXLmCAQMG4IcffoAQIs81BPAygOsAlplMcLcIw5UrVxASEoL4+Hikp6eja9euyMjIAABERERg8uTJ+OGHHzB37lyb9lEikeDRy0cA3gBwHsB9mNNqnreUIwCG5mcaUhilqJaPevfureYW2Lt3L00mkxrOQQMwFeDHuZZPVqxYoX4eNWoUr1y5QkVR+NZbbxVJHwqCS5cu0c/Pj35+fixbtiwBsF69ejnyVucuETDnd96kKFSyHff29qaHh0ce+4LRaGTr1q2p1+t54MCBIuytRPLsgH+7fETyE5LBAP4fyTIkgy2lCslZBS1OxZ3sM4XKlSvj5s2buHvXvF3jBZjji+/Pdr6iKDh+/Lj6vVmzZliyZAlMJhP69etnm0YXAoGBgVi7di2Sk5NhMBjg5eWFXbt2YfPmzXjttdcAIM+s4SCA1wA0N5nwocEAIQQ0Gg2SkpKQnJyMMmXK4NNPP8V3330HwPzslixZAh8fH3Tp0gWJiYk27qVE8hyTH+UAEAbgJQC9rCU/1xVGKaqZgre3N4UQdHV1JUlu27ZNfePtZzGols/1Fmx1wdTpdExNTWVISAgbNWpUJO0vaDZs2ECNRsP69eurIbf79eunRlV9mOH5C8tz6m1vTyEENRqNun+hVKlSdHBw4F9//aXeY/fu3dRqtezYsWOx3s8hkTwNoAANze/CnG3tMwCNAUwD0L4gham4k5KSgvj4eJBEmTJlAOQ0PNcAkAQgLts1VapUwalTp6AoCurVq4dDhw4hLi7uqZ4lZCc6Ohrz5s3Dzp070bhxY2i1WixcuBDh4eGoWLHiQ43EQwHsAjD73j1UVxQYjUb8/fffCA4OxuXLl6HX6/Hiiy8iOTkZAFCnTh1MmzYNq1atwscff2zbDkokzyn5CXPRGUBTAH+T7AugCgC7/FQuhGgphDgthIgTQox+xDkvCSFOCCGOCyG+zXfLbcjZs2fVzxEREQDMnkdWagA4AIDZlk0qVKgAkjCZTGjWrBnWrVsHjUaDjh072qrZhU7//v0xbtw4rF27FjExMRBCYPr06ejWrRu8vLzynJ8B4EUAdwCsMBrhaTl+/vx5+Pj4gCTOnDmDfv36qaIybNgwdOzYEW+++Sb++OMPW3VNInl+edJUAsA+y78HYY57JAAcz8d1GgBnAZQBoIfZQP1CrnPKwey042b57v2keoti+eiHH35Ql0CWLl1KkoyKiiJgzjiWAXBKrqWShg0bqp/37NnDyMhI1q9f3+ZtL2xMJhP79u1LAOzatava5ylTpqib9nKXCID3AG4Vgtps+zocHBwYGBiYZ7f3nTt3GBwczICAAN68ebMIeyuRPL2gAHc0fw7AFcAQAGcsg/iifFxXB8DGbN/fAvBWrnOmARiQn4ZaS1GIgjUhDACeP3+eJNV8xDUt6+Qx2QY9Ozs76vV6CiHo7OzM69evUwjBCRMm2LzttiAjI4MtW7akRqNhx44dVZvC5MmTH7rbGQBftjy3eZZ8E4qiqCFEgoODqdFouHnzZvUeBw4coJ2dHZs0acLMzMwi7K1E8nRSIKJgmRUEZvteGkB4vio2LzvNz/a9J4BZuc5ZZRGGXQD2AGj5pHqLQhT69eunBoQzGo1MTU1VB7dXLYNbYLYBz/q2a80V8P3336szhmeVlJQUVq9enfb29mzcuLHZVVej4ZgxYx7pqvqh5dm9bjFUCyHo4+NDAPT396erqytjY2PVeyxevJgA+MYbbxRhTyWSp5P8isJjbQqWilZl+36B5NHHXZONvLuZzINBdrQwLyFFAegOYL4QwjVPRUIMEkIcEEIcuHnzZj5vX3DExcWBJFxcXKAoCs6cOaP+VgPA3wCuW0JqA4CDgwMAICMjA02bNsWmTZvg6uqKyMhIG7fcdjg5OWHt2rUoWbIkjh49isjISBiNRkyfPh1Dhgx56DWjAWwAMD0tDa0tu6Jv3LgBX19f3Lp1CyTRrl071SW1V69eGDZsGD755BMsWrTIdp2TSJ4j8mNo3iOEqPHk0/JwBUBgtu8BAK495JyfSWaSPA/gNMwikQOSX5CMJBn5MANmYRMbGwuSCAw0dye359F+AMhmZL59+7b62SoKzZo1U3MxPKuULFkSGzZsAADcunULoaGhyMzMxJIlS9CpU6c85xsBdIX5j770/n1UUsz/OV6/fh2Ojo7QarWIi4tD9+7dYTQaAQAffvghmjZtiiFDhmDv3r026plE8vyQH1FoDLMwnBVCHBVCHBNC5Ge2sB9AOSFEsBBCD6AbgNW5zlllqR9CCE8A5WHePV1sSE1Nxd9//w3A7FEEPMjAZg+gAswW+Kwsc4ALIQRu3rwJg8EAX19fCCFw5coVtGjRoghab3vKly+PX375BfHx8dBqtfDz80N6ejp+//13NGjQIM/5yQDaAEgHsNpkglXyb9++jbS0NAQEBGDDhg148803AZiTFy1btgz+/v7o1KkTrl3L/Z4hkUj+DfkRhVYwexA1AdAOQFvLv4+FZBbMrukbAZwEsJzkcSHEBCGEdZ/DRgAJQogTAH4HMJJkwj/vRuFx7twDjbIu/xw6dAgAUBHmB/hXtvMdHR0BmAPiNWnSBJs3bwYANG/e3BbNLRbUrl0bK1euxOnTp+Hn5wcXFxekpKTg9OnTCA8Pz3P+JZg3vpQE8DMAg+U4SVy6dAkhISH46KOP1CUjDw8P/Pzzz0hOTkZMTAzS09Nt1DOJ5DkgP4YHAPUB9LV89gIQnJ/rCqPY2tC8cuVK1TC6d+9ekmRwcDABsKfFUFoxm/HUz89P/bxo0SK2bt2a5cuXt2mbiwvWjHONGjWivb09NRoN/fz8WLp06YcanmMsz/MHjYbCcszFxYUAWKZMGep0Om7fvl2t3/q36dmzp9zxLJE8ARTwjuZRMLuUAoAOwDcFokhPAdljHoWHh8NkMuHSpUsAgEowRwu8bPdgL5/52ZupX78+tm7d+twsHeWmR48e+PTTT7Ft2zY0bNgQQgjEx8cjKyvroZvbfgQwEkBnoxEfWvJbJyUlwcPDA+fOnYOXlxc6dOigbhyMiYnB+PHj8fXXX2PChAk27JlE8uySn+WjTjDP7lMBgOQ1AM6F2ajiRFxcHIQQsLe3h8FgwKVLl1SjZxjMRtL7JpN6/q1bt6DRaFCmTBlcvnwZaWlpz60oAMBrr72GcePGYePGjWjdujVMJhOuX78OBwcHODk55Tl/OoB5AEZkZOC/BvNCUkJCAjw8PHD79m0IIdCqVStcv34dADBu3Dj06dMH7733HhYvXmzDnkkkzyb5EYUMy9TDvHFBCMfCbVLx4syZMyAJb29vADk9jyrBnFQnMzMTgNkImpmZCUVRUKtWLWzatAlarRZRUVG2b3gxYvz48Rg8eDBWr16NTp06wWg04sqVK/D09ISdXd6IKa/CbFv4MD0dPS2/JyQkwGAwgCTi4+PRpk0bpKSkQAiBefPmoWnTphgwYAB+++25y/8kkRQo+RGF5UKIeQBchRADAfwK4MvCbVbxITY2FgAQHBwM4IHnkSPMO/myG5lLlCgBwCwSERER2LRpE+rWrQtn5+dmYvVQhBCYPXs2unfvjpUrV6Jjx44wGo24fPkyfH19odVqc5xvhNlVbReA+ffvo5VlKSkxMRGZmZlwdXXFkSNH8NJLLyEzMxN6vR4rV65EhQoVEBMTg7/++itPGyQSSf54oiiQnA5gBYCVMLuMvkPys8JuWHHg3r17qsuj1Wtm/35z1oQXLOccz3a+3jJ4AUBISAgOHTr0XC8dZUej0WDx4sXo0KEDVq1ahXbt2qnCEBgYCEXJ+Z9iOsxrlqcBLMvIQE3LHo+7d+8iISEBQUFB2LBhA4YMGaJuLFy7di0cHR3Rpk0bdXlJIpH8M/IzUwCAYwB2ANhu+fxckN0dtW7dugCAI0eOADDbEwDgVLbBzJp0BzAvdwDPlyvqk9DpdFi2bBmio6OxZs0atGzZEkajEZcuXUJgYGCe5DyJAFoCSACwxmhEiOV4VlaW6qq6cOFCvPvuuwCAoKAgrF27FgkJCWjTpo0aglsikeSf/HgfDQCwD0AMzPGM9gghno2kAE8gu+dRo0aNAAAXLlwAYLYn3AMQl83IfPfuXdjb2yM0NBQ7d+6Em5ubGmpbYsbOzg4//vgjGjVqhE2bNqFp06YwGo24evUq/Pzypv6+BiAa5pgpGwH4wrwHRFEUxMXFoXz58pg4cSI++ugjAEC1atXwww8/4NixY2jXrh3u3btnw95JJM8AT/JZhXkG75HtuweA0/nxdy2MYst9CtOnT1cjeJJkQkKC6lO/HuDBbD72BoOBAGhvb8/u3bvT39+fXbp0sVlbnzaSk5NZu3ZtarVaNmvWTM1QZ83ElrtEAkwGeAKgV7ZotADUDHdffPGFWv93331HIQTbtGnDjIyMIuypRFI8QEHtU4A5PlFKtu8pAC7/WzF6GrDOFFxdzTH6Tp8+rf5m9TyyYt3JfO/ePQQGBuLq1avSnvAYnJ2dsX79elSuXBnbt29HkyZNkJmZidu3b6ueXtk5AHM4jCCYPR3cAdy/fx/29vY4efIkypcvj8GDB+P7778HAHTr1g1z587F2rVr0atXL9WNWCKRPJ78iMJVAHuFEO9ZNrLtARAnhPivEOK/hdu8osUaDTUgIADAA88jF5gj/WX3cTFlW0ZKS0sDIO0JT8LV1RW//vorwsLCsGPHDjRq1AgZGRlITEx86Oa2HTAbn8sD2Azz3+HevXtwdHREbGwsypQpg549e2LNmjUAgEGDBuGDDz7A999/j1dffTXHxkKJRPJw8iMKZ2EOXGf9P+pnANdh3sD2TPtaWnfOhoaGAgB2794N4IHn0Yls5yYmJkKv10MIgXPnzqFcuXIoVaqUDVv7dOLu7o5ff/0V4eHh2LVrFxo0aICMjAykpKSoM7TsbIF5N2UlmG0MzjAHLXR0dMTZs2cREBCALl26YMuWLQCAUaNGYfTo0Zg3bx7GjBljw55JJE8p+VljKk7FVjaF9PR0CiEIgFOnTiVJvvDCCwTAAZYYPWUVRU0OA4Bubm4MDQ2lt7c3e/XqZZN2PivcuXOHNWrUoFarZb169QiAjo6OdHJyeqiNoT3MaVB3AHS0HHNwcKAQggEBAXR0dOTWrVtJmlOGDhkyhAA4fvz4Iu6pRFI0oABjH0UKIX4SQhyyhM4+ms/Q2U8158+fV5cbrJ5HFy9eBGB2R00FcM6yZGRvbw/AvGmtYsWKiI+PR61atWze5qcZV1dXbN68GREREdizZw8aNWqE1NRUKIry0F3PqwH0gDnn6zoATjAv2xkMBly/fh1ubm5o1aoVtmzZAiEEZs2ahd69e+Pdd9/Fe++9Z9O+SSRPE/lZPloKYBGAF2EOmW0tzzTZ3VEjIyNx//59pKamAnhgZLaup+l0OgBml1TrrubatWvbsLXPBi4uLti4cSNq1KiBHTt2oEmTJkhOToa9vX2eXc+AeUdlDwB1kdPGoNfrcePGDXh5eaFNmzbYvHkzNBoNFixYgL59+2L8+PF45513pI1BInkI+RGFmyRXkzxP8qK1FHrLihirKBgMBmg0Gpw9e1b9LbfnUVpamrojNz09HQaDAZUrV7Zha58dXFxcsGnTJjRq1AhbtmxBkyZNkJiYCBcXlzyb2wBgOcxvK9Vgtjd4wCwMWq0WV69ehbe3N9q1a4eNGzdCo9Fg/vz5GDBgACZOnIhx48ZJYZBIcpEfUXhXCDFfCNFdCBFjLYXesiLGKgpWLxjrTmZ3mDdQZReFzMxMddC6cOECIiIi1NmD5J/j7OyMdevWoUOHDtiyZQuioqKQkJAAT09PAMgjDqsBdIA56dE2mJP13Lt3DzqdDleuXIGPjw86dOiAdevWQVEUzJs3DwMHDsTkyZPx9ttvS2GQSLKRH1HoC6AqzBEHrEtHbQuzUcUBq+eRNRDe+vXrAZhnCQBwItfA5OjoiNDQUBw5ckQuHRUABoMBK1asQM+ePbF161Y0aNAAN2/eRMmSJc3GsFyxkjYCaA2gFMzCEADzrE2r1eLy5cvw8fFBp06d8OOPP0JRFMydOxeDBw/GlClTMHLkSCkMEomVJ1miARzLj8XaVsVW3kfWnbVDhw4lSVaoUIEAOMTieVTK4nmk0WgIgK6urmzdujUBcPny5TZp4/OA0Wjka6+9RgCsV68eFUWhj49PjmefvdQBmAjwPHqKqWkAACAASURBVMDylmM6nY4AWLp0aSqKou58NhqNHDp0KAGwb9++zMzMLOLeSiSFBwpwR/MeIcQLTz7t2SEzMxM3btwAANSpUwcAcPXqVQDmmUISgIsWzyPrUkZiYqLqhSQ9jwoORVHwySef4J133sGuXbsQERGBpKQkeHp6wmg0QmOJnmrlD5iTiTvAHHq7Jsx/T41GgwsXLiAkJASDBg3C+++/DyEEPv30U7z77rtYtGgROnfuLPM9SyRPUg0AJwFkwBwD6SjMUVKP5kdxCqPYYqYQGxurvnleu3aNJNW9CFsA7sr1duro6EgAjI6OZsmSJWW+4ELi888/p6IorFChAt3c3Oji4kKtVkutVptnxlAW4BmAqQDbWo4pltld+fLlCYDDhg2j0WgkSX722WcUQrBRo0ZMTEws4p5KJAUPCnCm0BJAOQAt8MCekC+XVCFESyHEaSFEnBBi9GPO6yyEoBAiMj/1FjZWI7OiKPD19UVCQoK65hyGnEZmAPDw8ICiKDh79ixq1679UC8Zyb/nP//5D1atWoVLly7B3t4eDg4O0Ol0cHBwyDNjOAuzq+oJAD8B6IcHoUhiY2MREhKCjz/+GL169UJmZiaGDh2KpUuXYteuXYiKilJnihLJ80Z+kuxchDnUTxPL57T8XCeE0ACYDaAVzJEhuj9sGUoI4QzgdQB7/1nTCw+rKFgzpm3duhUA4GUpuUXBzs4O5cqVQ1xcnFw6KmTatWuHrVu3IisrC2lpafDz80NqaupDg+jdBBAFcwC9BQDGZvstLi4OwcHBWLp0Kdq0aYP/396dx1Vdpo0f/9zsIIiCCoIoCiruhJiik1ZmbpXpOI9WtoyZbfZMv5qnfdKpZn7NNM3U2PLY1K/Fdi1LLbXUct8VE1Nc2AUVEUEQ2c79++M+5+sBQcmAA3K9X6/vi3O+5+vhOnUO17m36y4oKOCWW25hyZIlJCcnM2zYMGvXPSFakrr8cZ8NPA48aT/lCXxYh+e+EjiktU7RWpcBn2JmDlb3PPB3zGZbTUL1QnjffPMNcG7m0b5qM19yc3OtayUpNLxBgwaxefNmQkNDyczMpH///uTk5NRYa6oY06x9H/NGewdw7I+Xnp5OeHg4q1evZujQoaSmpjJmzBhWr15NQUEBQ4YMsb4QCNFS1KX7aCKmOGUxgNY6m7oVwgunaontLPs5i1LqCiBCa720TtE2kp07dwLnCuFt3boVOLfbmmM6qmNa5KlTp/D09EQpRXx8k+gBu+x17dqVDRs2MHjwYHbt2sXgwYNJT0+vMTFUAHcBf8Z0I60E2mG6k44ePUrbtm3JyMhg8ODBbNq0iSFDhrBlyxZCQ0O5/vrreffddxvzpQnhUnVJCmX2QQoNoJRqVcfnrqlj3ZoMrpRyA/4FPHrRJ1JqplJqu1Jqe25ubh1//aXbv38/gLVrWmamyW19gJNAlr02v81ms/ZROHXqFH379rW6nETDCw4OZuXKlcyYMYMtW7YwYMAAjhw5QseOHXF3dz9vbGcOMAWIx2wl2AeorKzk1KlTlJeX4+bmxjXXXMMnn3xCt27d2LhxIyNGjGD69Ok89dRTVcqjC3HZuthINPBHYB6QAtyDmfX3UB3+XQKwwun+k8CTTvcDgRNAmv04i9l9Mf5Cz9vQs4/Onj1rzWBZt26d1vrczKO1oNdUm+USHh6ulVK6bdu2esaMGQ0am6iZzWbTc+fO1e7u7rpLly46MDBQBwYG6uDg4Fp3cTti38ltfLWZSV26dNGAnjNnjq6srNRlZWX63nvv1YCePHmyLi4udvXLFeKSUI+zj9pjao99AfQEnsUsGL2YbUB3pVRXpZQXMBVTkcCRjAq01u201pFa60jM5j03aa231+G5G0xKSop1Oz4+npycHGvmUfWaR2CK4UVHR5Ofny/jCS6ilGLWrFl89913nD59Gq01bdu2JT8/31qR7mw7MAgzx3oxZsDM0QpIT08nMjKSOXPmMGnSJM6cOcObb77Jyy+/zBdffMHQoUOrvEeEuNzUJSmM0lp/r7X+H631H7XW32NmFF2Q1roCmIWpQLAP+FxrvVcp9ZxS6qZfF3bDccw48fLywsfHh1WrVgGm3lEQVTfWATh58iQhISGADDK72rXXXsu2bduIiIggPT2dvn37kpqaSnR09HnXZgPDMQX1XsRMWw20P5aWlkZERARLly5l0KBB7N27l0ceeYRvv/2WjIwMBg4caJU9EeJyU2tSUErdr5TaA/R03kdBKZWKWcR2UVrrb7XWPbTWUVrrv9jPPau1XlzDtVe7upUA5/ZhdhTC+/7774FzM4/2V5sPX1hYCIC/vz+9e7eohd9NUrdu3di8eTNTpkzhp59+olevXqSlpdGxY0d8fHyqXFsC3AI8jNn/eTvQ3/7YkSNH8PPz48SJEwwZMoTPP/+cMWPGsH37diIjIxk/fjzPPfecjDOIy86FWgofY2bzLabqPgoDtdbTGiE2l0hMTASge/fuAGzZYpZPOJLCT05/BBybv+Tm5jJo0KDzFlAJ1/D39+fjjz/m9ddf5/DhwwQFBVFaWopSis6dO593/auY9Qy+mD7MOzHdSUVFRZw6dYq2bdsyZcoUHn30USIiItiwYQPTpk1j9uzZTJgwgVOnTjXq6xOiIdWaFOx9/mla61u00z4KWuuTjRlgY3NMR42NjQUgIyMDMEkhFziuz1XTDAoKAuDw4cPSddTEKKV44IEH2LBhA35+fhQWFhIWFkZGRgYxMTHnXb8RsyfDJuA94C3Axz7wlpWVRUREBP/85z8ZMWIEx48f5/3332fu3LksX76cuLg468uDEM1dXcYUWhRHEhg+fDgFBQWUlJQANQ8yBwQEEBQUREVFhSSFJio+Pp6dO3cyduxYDh8+TM+ePTlw4AAdO3bE39+/yrW5mFouf8VMs9sODLA/lpWVRdu2bdm9ezexsbEsWLCAWbNmsXbtWmw2G7/5zW948cUXpTtJNHuSFJw4J4GhQ4daXUlg6nRUTwqlpaXWxi+SFJqutm3b8tVXX/Gvf/2L1NRUAgMDqaio4OzZs3Tr1q3KtZXA08B1QBtM7ZWHAbQmPz+fkpIS/P39mTJlCvfccw/9+/cnMTGRiRMn8uSTT3L99deTk5PT2C9RiHojScGJo7yFm5sbHTp0YPt2M+4dhvkDUX1jnezsbKufumPHjo0crfgl3NzcePjhh9m2bRvh4eHk5ubSpUsXUlJS6NGjx3mb9qzCDDovx6ywXIbZ0U1rbS2Qe/vtt4mPjyctLY3PPvuMt99+m40bN9K/f3+rNIoQzY0kBSeOmUdt27ZFKcW6deuAmmseeXh4UF5ezokTJ6SV0Iz079+fbdu28Yc//IHDhw/TsWNHMjMz8fPzIyIiosq1ecDNwH3AVZgpdxPtj+Xk5ODt7U1OTg7x8fG88MIL3HHHHezYsYOwsDBuuOEGZs6cac1OE6K5kKTgxDHI7OhScLQUHElhj728BZybspqXl8eVV17ZeEGKX83Hx4dXXnmF5cuXo7WmrKyM4OBgMjMz6dGjx3nXzwMGYgp5fQl8gqmdVFpaSkFBAR07duTZZ58lISGByspKtmzZwmOPPcY777xDv379WLlyZaO+PiF+DUkKThyF7wYMGMDZs2er7LZ2HFOTwyEoKMiqreOokSSal9GjR5OUlMStt95Keno6oaGhpKen07p1a8LCwqpcux8YjCm/PQkzvjTZ/lhWVha+vr4kJycTFxfHK6+8wl/+8hc2bNiAr68vo0aN4v777+f06dON+vqEuBSSFJw4VjMPGzaMpKQk63xNM488PDysQeYrrriikSIU9S04OJgPPviApUuX4u7uTllZGf7+/mRnZ9OzZ88qa08qgL8AcUAGsMB+dABKSkooKioiKCiIJ598kmHDhuHr68uuXbt49NFHmTdvHv379+e7775zxcsUos4kKdhprcnLywPMvsy7du2yHqtp5lFubi7e3t5ERUXRpk2bxgtUNIjx48ezd+9eZsyYQXZ2Nu3btycjIwNPT8/zynHvBYYAT2BWc+4DZmLKAh87dgwvLy9+/vln4uLiePbZZ/nzn//MunXr8PLyYvTo0dxyyy0cPXq0sV+iEHUiScHu6NGjVNrHDKKioqykEI6pieNc88jNzY3s7GyKi4uJi4tr9FhFwwgMDOStt97i+++/p3Xr1pSUlBAcHEx6ejoRERHWCnYwU1f/hlnHsBsz7rAJiAXKysqsVsM//vEPevfuzcmTJ9m9ezdz5szhyy+/JCYmhjfffFPWNYgmR5KCnWPmkb+/Px4eHuzYsQM4N8jsPB3V0W2Un58v4wmXoeuuu46kpCTmzJnDiRMn8PHxobCwkLKyMmvjJYdk4FpgGhCJWfD2CmYXqhMnTqCUoqioiJtuuolbb72V3//+9+zZs4eBAwfywAMPnLceRghXk6Rg50gCnTp1orKykt27dwNOM4+cyluEhoZat6WlcHny8fFh9uzZJCUlMWLECAoKCmjXrh2HDh0iICCgynsA4CMgBtNieAiTLO4C0JqTJ0/i6+vL0qVL6dGjBx9++CFff/018+fPJyUlhbi4OGbOnMmxY8ca+VUKcT5JCnabN28GzMyjAwcOUFpaCpikcAyz45pDQEAAXl5mp19JCpe36Oholi1bxsKFC/H29qayshI/Pz+OHj1KWFhYlcqrp4AHMeMN6cC7mJbDcMxAdHl5Oa1ateL5558nJiYGpRTJyck8/PDDvPvuu3Tv3p2XXnrJeu8J4QqSFOwcLYPqg8w1zTwqKioiICCALl26EBwc3HhBCpdQSvHb3/6W5ORknn/+eU6fPo27uztnzpyxSmU4b/25DbPt4C2Y9QxrMDtUdcPsvwFQXFzMtGnTGDduHP/1X/9FUlISw4cP57HHHqNPnz4sWrTI2txJiMYkScHOsSahLjOPMjIyKC8vl1ZCC+Pn58czzzzDoUOHuOuuuzh16hS+vr5WraPIyMgq13+K2arwaUyhvX2Y8YYOmD29lVL89NNPJCQk8MQTT/C3v/2N5cuX4+XlxaRJkxg2bBhr1qxpzJcohCQFgPLycs6cOQNATEyMtYgtAmhN1aQQFBREfn4+hYWFMsjcQjnqHiUmJjJ06FBKSkpo1aoVWVlZeHp6Wjvxgdl4/K9Ad0xJ7gcxm52/ALTWmjNnzuDh4cGyZcvo168fn376KYsXL2bevHmkp6dz9dVXM3bs2CpfVIRoSJIUMNsvgtmC09vb20oKjn3UnJOC8yYt0lJo2QYMGMDKlStZuXIlffr0oaKiAh8fH44fP46fn1+VrsWjwL1AL8yuVU8DqZi1Dl4VFZSWluLh4cH8+fPp3bs3e/fuZe3atfz9739ny5YtxMXFMXXqVPbt2+eKlypaEEkKnKtx1KFDB3bu3MnZs2cBp+moTtd26NDBui1JQQCMHDmSTZs2sXjxYiIjI9Fa4+HhQV5eHn5+flUWNx4CbsWsb1gP/F8gDZMcfMrLrbUyc+fOpXfv3mRkZLBu3TqeeeYZli5dSp8+fZgyZQp79uxp7JcpWghJCsCGDRsA6N27N+vXr7fO98F8w3OeeVReXo6vry/h4eFVuglEy6aU4sYbbyQxMZFPPvmE8PBwwCx0LCgowNfXl4CAAOv6n4CbMAPSWzHJIR2YAwSUl6O1xmaz8cYbbxAbG0tubi5r1qzhiSeeYNmyZfTv35+JEydaRRyFqC8NmhSUUmOUUslKqUNKqSdqePwRpdTPSqmflFKrlFJdanqehuZoKQwbNowNGzbg6ekJ1Dzz6NixY7i5ucl4gqiRm5sbU6dOJSkpiYULFxIdHY3W2lrE5uXlRWBgoHX9ZuAGTBXW1cBsTMvhb0BIRQU2mw2bzcZ//vMfBg8eTGZmJkuXLmX27Nn8+OOPDBw4kLFjx7Jq1SqZrSTqh7bvQ1vfB+AOHMbMxPPCVAPoXe2aawA/++37gc8u9rwDBw7U9a1du3Ya0CtWrNDBwcHa3d1dA7oQ9KugsR9t27bVnp6eGtBz5syp9zjE5cdms+lvvvlGDx06VAPa29tbu7u7a6WUbtOmjfXechx9QH8EugJ0Gej5oGPtjymlrPfmqFGj9JdffqlfeOEFHRISogEdGxur58+fr8vKylz9skUTBGzXdfjb3ZAthSuBQ1rrFK11GWaG3gTnC7TWP2itz9jvbgY6NWA8tTp16hRgBprz8vKorKykM6ZUgXNLISYmhvLyckDKZYu6UUoxbtw41q9fz5o1axgzZgyVlZUopSguLgbMjDaHvcBtQBTwGuYDswvTihinNTb7mMOPP/7IpEmT+PTTT3nuued44403KC0t5fbbb6dbt2689NJL1poIIX6JhkwK4Zh9SRyy7Odqczdm18NGVVRUREVFBW5ubtZ2nHBukNk5KTj6iUEGmcUvo5Ri+PDhfPXVVxw4cID7778fDw8PAGs6tL+/v7UILh14BDMt+o9ANLAUOGi/39r+5WT//v3ce++9PP3004wbN463336bHj168NhjjxEeHs706dOt7lEh6qIhk4Kq4VyNnZ5KqWlAPPBSLY/PVEptV0ptz83NrccQz+221qZNGzZu3Iivry9wbjqq88wjd3d3lFKEhITInsziknXv3p3XXnuNrKwsXnzxRWtGW1lZmWm+u7lZFVkLgJcxfbBTMd+sXrL/fA+Iq6gATGv35ZdfZsaMGfj4+PDaa69x55138vnnnzNo0CAGDx7MBx98QElJSWO/XNHMNGRSyMJ80XHoBGRXv0gpdR1m2vZNWusai75ord/SWsdrreMd22DWF8eK0aioqPNmHuUA+U7Xnjx5Ei8vLwYOHFilrIEQlyIoKIjHH3+clJQUlixZwnXXXQeAzWazpqY6vqRUAJ8BVwN9gXcwO8BtwXQvzdIaRyfU8uXLmTVrFkuWLOG+++7j2WefpbCwkDvvvJOwsDAefPBBduzYIQPTomZ1GXi4lAPwwCze7Mq5geY+1a65AjMY3b2uz1vfA83jxo3TgJ4xY0aVAb8toL93ut++fXsdERGhlVL6mWeeqdcYhHBISUnRjz/+uG7fvr0GrIkNzoPMjsMf9H2gt4HWoM+C/gz0GNBu1QawR4wYoZ966ik9ZcoU7ePjowHdv39//eqrr+rc3FxXv2zRCKjjQHODJQUTA+OAA/Y//E/bzz2HaRUArMQUIU20H4sv9pz1nRQiIyM1oB955BHrA6RAnwb9itOHavTo0dbtRYsW1WsMQlRXVlamv/76az1hwgQrGXh4eFRJFM5HP9D/BH3cniCyQL8E+gqnmUuA9vPz01OnTtUPPfSQHjhwoPW8N9xwg/700091cXGxq1+6aCBNIik0xFHfScHX11cD+vbbb7c+dF3sH6x7nD509957r3U7PT29XmMQ4kKOHTumX375Zd23b1/rPejm5qaB81oPnqAngv4KdKn9fbwP9J9AR1dLEO3bt9e33XabnjZtmg4PDzetD39/fccdd+jly5fL1NbLTF2TgjLXNh/x8fG6vmZTaG0G9cBMMU1OTqaoqIhxwDfAMGCj/doHH3yQ119/naCgIGtHLSEaW1JSEp988gnz588nMzMTpZSjVV7lNkBb4LeYshojMAOIO4GFmFLeB5yeNzw8nKFDh3L27FnWrl1LQUEBQUFB3HzzzUyePJmRI0dae4iI5kkptUNrHX/R61pyUtiyZQtDhgyhVatWlJSUWJnyf4C/Yz5Up+zXTp06lQULFjBy5EhWrFhRL79fiEultWbz5s18/PHHfP755xw/fvy8pOCsEzAFkyQS7OeSMAliEabshkNoaCjx8fGcPXuWrVu3UlhYSGBgIBMmTODmm29m1KhR+Pv7N+CrEw2hrkmhRdc++uijjwBTCtlms1kfqL6YaVKOhNCxY0f27NmDzWaTRWuiSVBKkZCQwNy5c8nJyWH9+vU8/PDDVdbSOMvCTG0dikkQ/w3kAc9iZoCkYRbLjQZOHj3K0qVLWblyJZWVlQwfPpzY2Fi+/vprJk2aRLt27Rg/fjzz5s2z9iERl48W3VLo06cPP//8MwkJCWzatMk6n4JpZk+2358wYQJLlizBZrOxYMECJk+eXNPTCeFyWmt27tzJokWLWLRoET//bFba1NaK6ICpvXQjMApoBRQB3wHLgRVAhtP1ffv2JTAwkPT0dLKysgCzkHP06NGMGTOGhIQEq3aYaFqk+6gOvL29KSsrIzY2loMHD1JcXEwXzLemhzDfnODceAKYvRe6dHFJ3T4hfrGMjAy++eYbvvrqK1avXk2FfbFbTXwwxchuBMYDjp1D9mOSwwrM1qKOujSBgYFERUVx5swZDh48SGVlJQEBAVx33XWMHj2akSNHEhUVJeNvTYQkhYvIzMy0Nszx8/Pj7Nmz2Gw27gDeB/ph+lwBHnroIebOnUtoaCjZ2dnyJhfN0pkzZ1izZg3Lly9nyZIlpKamXvD6GEx30hjMQLUvUIZZMLca+AFTsMyx4jQ0NJTg4GCOHz+Oo/JA586dufbaa62jtu4t0fDqmhQ8GiOYpmj+/PkAhISEcOzYMev81cAJqtY8OnHiBG5ubowYMUISgmi2/Pz8GDt2LGPHjuXVV18lIyOD7777jhUrVrBy5UqrMKTDfvvxKqYV8RvgWvvxDKbMdwmwCVgHrDt6lM1Hj1Js//cdOnTA09OThQsX8t577wEQHR3NVVddxfDhw7nqqqvo1q2bfKaamBbbUkhISGDz5s0MHz6ctWvXWudTge3A7+z3O3fujJ+fH/v372fu3LnMmjXrV/9uIZoarTX79u3jhx9+YNmyZaxZs4aioqJar28NXAWMBIYDsZha+RWY8bj1mGSxGTPIDRAQEEBgYCD5+flWhdiwsDCGDRtGQkICCQkJXHHFFVbdJ1G/pPvoIvz9/SkuLubqq69m06ZNlJaWWuMJs4DX7ddNnjyZL774Aq01u3btIjY29lf/biGaOq01ycnJrFu3jlWrVvHjjz9WaVFXF4CZ6nqV/bgS090EcASTHDZjdpnbBZzGFJgMDAykoqKCwsJCwJSvj4uLIyEhgUGDBhEfH090dLS0JuqBdB9dQHFxMcXFxSilSE5OtvZIuNr++I9O13bs2BGtNb6+vvTr16+RIxXCNZRSxMTEEBMTwz333ANAdnY2W7ZsYf369axevZqkpCRr4Po0ZsbSd/Z/7wn0B4bYjwTMGgkAG5AMbK+sZPvJk+zETIs9jSkGuHfvXrZt22YVBWzTpg0DBw5k0KBBxMXFERsbS1RUlLXwVNSvFpkUFixYAED79u3Jycmxzl8D5HJ+uWyA+Ph467YQLVFYWBgTJ05k4sSJgNmvfM+ePWzfvp21a9eyadMmUlNT0VpTDuywH45Wd3vMtqPx9uNa4Han5z8EJFZUkHj6NInAHsx02IKCAjZu3Mjq1autabUBAQEMGDCA2NhYBgwYQL9+/ejTp48sqqsHLTIpfPbZZwAEBweTn59fpaXwI1U3fcjONtW+HWWNhRCGp6cncXFxxMXFMXPmTABKS0utRLFhwwa2bdtGSkoK5eXl5GLWPix3eo6OmPGIK+w/Yzm3PgigEEjSmqSSEvZgvrDtA3JOn2bTpk1s3LgRm81mXd+tWzf69etHv3796N27N7169aJnz55WCXJxcS1yTCEoKIj8/Hx8fX1p06YNOTk5RGIGmR8E3rBfFxUVhc1mIzU1ldWrV3PNNdf8uuCFaIFsNhuHDx8mMTGRbdu2sXnzZvbv309tG2YFYKaE93U6+gHtnK4pwCSHfZgZUsmcK8dcXm2hnlKKyMhI+vTpQ0xMDD169KBnz5706NGDkJCQFjNeIQPNtbDZbNYOalprPDw8qKio4C7gXczmOo7uo6lTp/LZZ5+hlKKwsJBWrVr9+hcghACgpKSEffv28fPPP7N161Z27NjBoUOHyM3NrXH1dQhmR8Re9sNx23kPxErMVqYHMd1Rh+3HIUylgrPVnrNVq1Z0796dnj17EhUVRXR0NNHR0URFRREaGnpZjVvIQHMtfvjhB8A0ff39/a3Nza8BjlN1PCE0NBStNVFRUZIQhKhnvr6+VvfTtGnTrPPl5eWkpKRw8OBBdu3axY4dO0hOTiYrK4sfior4odrztAa6Az2qHVdiilo6y8H0CKRiZhqmFheTlpjIzsREvqZq0vDw8KBTp05ERUXRo0cPunbtSpcuXejSpQudO3cmJCTkskoaDi0uKbz11luA2Q+3qKgILy8vysrKrPEEZ46+yhEjRjRmiEK0aJ6envTs2ZOePXtyww03VHns7NmzpKamcujQIRITE9m9ezeHDh0iIzubnSdOnNfCaAtEA1H2n5GYrSATMFVjq/8BPI5paWQCmRUVZKWlkZWWxp5Vq1iOKZTpWMHt5uZGSEgInTt3JjIykq5du9K5c2c6depEp06dCA8Pp127ds0ucbS4pODYkxnMNxKtNV0xdV5etJ/39PSkoqKCAwdMxfmxY8c2epxCiPP5+PjQq1cvevXqxY033ljlMa01eXl5pKWlceDAAfbs2UNycjJpaWmsysri87y8KoPS7piKsZ2BLvafjqMnZmFeYA0x5GGSQ7bNRk5ODtk5ORzbsoUUzP4rRzHbSRbYrw8MDCQkJISIiAgrgYSFhREaGmod7du3bzKL9lrcmIKbm5v1bcLPz48zZ84wHbMRem/MwFVoaCgdOnTgxIkTZGdnc/ToUUJCQuojfCGEC5WWlpKTk0NmZib79u1j3759pKamkpWVRVZWFnl5eVWKBgYA4ZjkEYEZvwhzOjraj5rqwpZiWh7HMUniOGbKey6mlM4Jp9t5QKFS+LZqRXBwMBEREYSFhREeHk5ERISVOLp3737JBTlloLkG+/fvp1evXued/wC4Hgh1OvfCCy/wpz/9iTZt2ljjDkKIluH06dMcO3aMrKwsDh48yKFDh0hLS+PIkSOkp6eTm5tLaanpSFKYbqrQakeHt+YYcgAAClhJREFUakcIZq1GbZNjK4F8TILIA07aj3ynn1uALZf4N1sGmmvw8ssvV7nv6+tLSUlJlfEEb29vKioqiI2NRWstZS2EaIECAgIICAggOjqaq6+++oLXlpaWkpeXx/Hjx0lLSyM1NZXMzEwSc3JISUkhMzOTvLw8ysrKAPDDTK9tbz/aAcE1HB0xsyHbAm3sv+uv9f5Kz9egSUEpNQZTZNEdeFtr/WK1x70xX9QHYpLjFK11WkPF8+2331a5X1JSQjdMs/BH+zlPT0+uv/56tm7dCsCYMWMaKhwhxGXA29ubsLAwwsLC6vwl0mazcfr0afLz88nOzraOI8eOsTE9nYMHD5KVlUV+fj4lJSW4YxKDTSmeatBX04BJQSnljlnhPgpTKHGbUmqx1tp51ufdQL7WOlopNRX4G2ZSQINwrE525liO9qP9Z1FREdOnT7daFTfddFNDhSOEaKHc3NwIDAwkMDCQyMhIV4dTRUPOlboSOKS1TtFalwGfAhOqXTMBs6cNmD3ER6oGWl5YW4XHqzGzBfZjKjR26NCB8ePHs3fvXmtqnBBCtBQN2X0Ujpnu65AFDK7tGq11hVKqANOddqK+g5kSGsoqp/tubm506NCBzseOsdQ+cFNeXs60adOorKwkPz+f7t27t5gl8EIIAQ2bFGr6a1p92Lwu16CUmgnMBKwtNC8lmCov1mbjxNGjHAf+1/GLteb3v/+9teo5ISHhkn6XEEI0Vw2ZFLIwY7gOnTBrPmq6Jksp5YFZK3Le/E+t9VvAW2CmpF5KMKtttvNWFjo22tFa4+bmRnx8PH379uXf//43ADfffPOl/CohhGi2GnJMYRvQXSnVVSnlBUwFFle7ZjFwp/32ZGC1bqCFE0opHnrooSrnioqKrIVsNpuN6dOnA1jbc44bN64hQhFCiCarwZKC1roCs7PlCsxC4c+11nuVUs8ppRxTet4BgpVSh4BHgCcaKh6A2bNn1/qYt7c3U6dOZf369SQnJxMVFdVklp0LIURjadB1Clrrb4Fvq5171un2WeB3DRmDs+DgYH7zm9+wfv16lFIopaxaKL/73e8ICAhg6tSpAHz44YeNFZYQQjQZzat8Xz14+umnAVNH3Waz4elpqpY41iYcOXKEwYMHM2TIEFeGKYQQLtGiylwAjBo1iqCgIE6ePElMTAz79+8nMjKSqKgoa/Xyf/7zHxdHKYQQrtHiWgru7u7cd999gJmCGhERwSOPPMLMmTMpKytj1KhR9OvXz8VRCiGEa7S4pABw9913A5CcnMy8efNo164dK1asAOCvf22MklNCCNE0tbjuI4Bu3boxYsQI1q9fz5w5c0hJScHDw4Nrr72W+PiLVpYVQojLVotsKQDMmDGDyspKtm7dysmTJ6moqOCZZ55xdVhCCOFSLTYpTJo0iYCAAFq3bo2/vz/Dhw/nqquucnVYQgjhUi02Kfj5+XHbbbdRWFhIYWGhNVVVCCFashabFACrrMWgQYMYNWqUi6MRQgjXa5EDzQ7x8fHMmTOH8ePHS4lsIYSghScFpdQF6yEJIURL06K7j4QQQlQlSUEIIYRFkoIQQgiLJAUhhBAWSQpCCCEskhSEEEJYJCkIIYSwSFIQQghhUVprV8fwiyilcoH0aqfbASdcEM4v1VziBIm1oTSXWJtLnCCx1lUXrXX7i13U7JJCTZRS27XWTX4jhOYSJ0isDaW5xNpc4gSJtb5J95EQQgiLJAUhhBCWyyUpvOXqAOqoucQJEmtDaS6xNpc4QWKtV5fFmIIQQoj6cbm0FIQQQtSDZp0UlFJjlFLJSqlDSqknXB1PbZRSEUqpH5RS+5RSe5VSf3B1TBeilHJXSu1SSi11dSwXopRqo5RaqJTab/9vm+DqmGqjlPo/9v/3SUqpT5RSPq6OyUEp9f+UUseVUklO54KUUt8rpQ7af7Z1ZYwOtcT6kv098JNSapFSqo0rY3SoKVanx/6olNJKqXauiO1Cmm1SUEq5A68DY4HewC1Kqd6ujapWFcCjWutewBDgwSYcK8AfgH2uDqIOXgWWa61jgAE00ZiVUuHAfwPxWuu+gDsw1bVRVfEeMKbauSeAVVrr7sAq+/2m4D3Oj/V7oK/Wuj9wAHiysYOqxXucHytKqQhgFJDR2AHVRbNNCsCVwCGtdYrWugz4FJjg4phqpLXO0VrvtN8+jfnjFe7aqGqmlOoEjAfednUsF6KUag0MB94B0FqXaa1PuTaqC/IAfJVSHoAfkO3ieCxa67XAyWqnJwDv22+/D9zcqEHVoqZYtdbfaa0r7Hc3A50aPbAa1PLfFeBfwGNAkxzQbc5JIRzIdLqfRRP9Q+tMKRUJXAFscW0ktXoF84a1uTqQi+gG5ALv2ru63lZKtXJ1UDXRWh8B/oH5ZpgDFGitv3NtVBcVorXOAfOlBujg4njqajqwzNVB1EYpdRNwRGu929Wx1KY5JwVVw7kmmXkdlFL+wBfAw1rrQlfHU51S6gbguNZ6h6tjqQMPIA54U2t9BVBM0+niqMLeHz8B6AqEAa2UUtNcG9XlRyn1NKar9iNXx1ITpZQf8DTwrKtjuZDmnBSygAin+51oQk3y6pRSnpiE8JHW+ktXx1OLYcBNSqk0THfctUqpD10bUq2ygCyttaPFtRCTJJqi64BUrXWu1roc+BIY6uKYLuaYUqojgP3ncRfHc0FKqTuBG4DbdNOdZx+F+WKw2/4Z6wTsVEqFujSqappzUtgGdFdKdVVKeWEG7ha7OKYaKaUUpu97n9b6n66OpzZa6ye11p201pGY/56rtdZN8hut1vookKmU6mk/NRL42YUhXUgGMEQp5Wd/L4ykiQ6KO1kM3Gm/fSfwtQtjuSCl1BjgceAmrfUZV8dTG631Hq11B611pP0zlgXE2d/LTUazTQr2gaVZwArMB+xzrfVe10ZVq2HA7Zhv3on2Y5yrg7oMPAR8pJT6CYgF/urieGpkb80sBHYCezCfuyazslUp9QmwCeiplMpSSt0NvAiMUkodxMyUedGVMTrUEutrQADwvf2z9b8uDdKullibPFnRLIQQwtJsWwpCCCHqnyQFIYQQFkkKQgghLJIUhBBCWCQpCCGEsEhSEOIi7NVYH7DfDlNKLXR1TEI0FJmSKsRF2OtVLbVXOBXisubh6gCEaAZeBKKUUonAQaCX1rqvUuouTPVQd6Av8DLghVmoWAqM01qfVEpFYcq8twfOAPdorfc3/ssQ4uKk+0iIi3sCOKy1jgX+p9pjfYFbMaXc/wKcsRfo2wTcYb/mLeAhrfVA4I/AG40StRCXQFoKQvw6P9j3yDitlCoAltjP7wH62yvjDgUWmLJHAHg3fphC1I0kBSF+nVKn2zan+zbM58sNOGVvZQjR5En3kRAXdxpTcO0Xs++bkaqU+h2YirlKqQH1GZwQ9UmSghAXobXOAzbYN2B/6RKe4jbgbqXUbmAvTXTbWCFApqQKIYRwIi0FIYQQFkkKQgghLJIUhBBCWCQpCCGEsEhSEEIIYZGkIIQQwiJJQQghhEWSghBCCMv/B3otj6nKSHq9AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(output.loc[output.index < 15,output.columns.str.startswith('ensemble')],'k')\n", + "plt.plot(output.loc[output.index < 15,'x'],'r')\n", + "plt.xlabel('time')\n", + "plt.ylabel('temperature')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Remember to remove the models from memory" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "x.finalize()\n", + "for ensembleMember in range(nEnsemble):\n", + " ensemble[ensembleMember].finalize()" + ] + } + ], + "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.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/run-model-from-bmi.-wth-different-states.ipynb b/examples/run-model-from-bmi.-wth-different-states.ipynb new file mode 100644 index 0000000..28b2b7e --- /dev/null +++ b/examples/run-model-from-bmi.-wth-different-states.ipynb @@ -0,0 +1,713 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Run the `Heat` model through its BMI" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`Heat` models the diffusion of temperature on a uniform rectangular plate with Dirichlet boundary conditions. View the source code for the [model](https://github.com/csdms/bmi-example-python/blob/master/heat/heat.py) and its [BMI](https://github.com/csdms/bmi-example-python/blob/master/heat/bmi_heat.py) on GitHub." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Start by importing `os`, `numpy` and the `Heat` BMI:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from heat import BmiHeat" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create an instance of the model's BMI." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "x = BmiHeat()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What's the name of this model?" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The 2D Heat Equation\n" + ] + } + ], + "source": [ + "print(x.get_component_name())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Start the `Heat` model through its BMI using a configuration file:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "# Heat model configuration\n", + "shape:\n", + " - 6\n", + " - 8\n", + "spacing:\n", + " - 1.0\n", + " - 1.0\n", + "origin:\n", + " - 0.0\n", + " - 0.0\n", + "alpha: 1.0\n" + ] + } + ], + "source": [ + "cat heat.yaml" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "x.initialize(\"heat.yaml\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Check the time information for the model." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{\n", + " \"time\": 0.0,\n", + " \"plate_surface__temperature\": {\n", + " \"value\": [\n", + " 0.7444227320828752,\n", + " 0.34022663250676577,\n", + " 0.8432555496724361,\n", + " 0.060846504441758764,\n", + " 0.8446697787107923,\n", + " 0.27627788779009277,\n", + " 0.5877718799413219,\n", + " 0.5737654859660629,\n", + " 0.34645112782286935,\n", + " 0.37414828806554334,\n", + " 0.7846696450611258,\n", + " 0.5304173216068899,\n", + " 0.0792535018455699,\n", + " 0.8793127450388358,\n", + " 0.15703439775879213,\n", + " 0.3912156458741114,\n", + " 0.9707058465142714,\n", + " 0.5511025163784199,\n", + " 0.1470283563347624,\n", + " 0.05645959012302082,\n", + " 0.9089090697986834,\n", + " 0.4495395294486181,\n", + " 0.011832496810315396,\n", + " 0.4171122020390934,\n", + " 0.33814771386922027,\n", + " 0.8466251708240787,\n", + " 0.6636025903978577,\n", + " 0.342616329035892,\n", + " 0.9154478512615142,\n", + " 0.14537348896585434,\n", + " 0.545772555177156,\n", + " 0.8827706254573504,\n", + " 0.2851195288529734,\n", + " 0.47448421167380395,\n", + " 0.2783120501807711,\n", + " 0.28066701069377975,\n", + " 0.841622847809508,\n", + " 0.7140017886819985,\n", + " 0.9744899247558164,\n", + " 0.45762952067167606,\n", + " 0.22664959884228209,\n", + " 0.7177555043888103,\n", + " 0.9900324558495176,\n", + " 0.7858216954853029,\n", + " 0.11380231880676006,\n", + " 0.16366335988935754,\n", + " 0.7884014952466945,\n", + " 0.14975089428180044\n", + " ],\n", + " \"type\": \"float64\",\n", + " \"itemsize\": 8,\n", + " \"nbytes\": 384\n", + " }\n", + "}\n" + ] + } + ], + "source": [ + "stateOut = x.get_state()\n", + "print(stateOut)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "stateOutDict = json.loads(stateOut)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'time': 0.0, 'plate_surface__temperature': {'value': [0.7444227320828752, 0.34022663250676577, 0.8432555496724361, 0.060846504441758764, 0.8446697787107923, 0.27627788779009277, 0.5877718799413219, 0.5737654859660629, 0.34645112782286935, 0.37414828806554334, 0.7846696450611258, 0.5304173216068899, 0.0792535018455699, 0.8793127450388358, 0.15703439775879213, 0.3912156458741114, 0.9707058465142714, 0.5511025163784199, 0.1470283563347624, 0.05645959012302082, 0.9089090697986834, 0.4495395294486181, 0.011832496810315396, 0.4171122020390934, 0.33814771386922027, 0.8466251708240787, 0.6636025903978577, 0.342616329035892, 0.9154478512615142, 0.14537348896585434, 0.545772555177156, 0.8827706254573504, 0.2851195288529734, 0.47448421167380395, 0.2783120501807711, 0.28066701069377975, 0.841622847809508, 0.7140017886819985, 0.9744899247558164, 0.45762952067167606, 0.22664959884228209, 0.7177555043888103, 0.9900324558495176, 0.7858216954853029, 0.11380231880676006, 0.16366335988935754, 0.7884014952466945, 0.14975089428180044], 'type': 'float64', 'itemsize': 8, 'nbytes': 384, 15: 0}}\n" + ] + } + ], + "source": [ + "stateOutDict['plate_surface__temperature'][15]=0\n", + "print(stateOutDict)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "x.set_state(json.dumps(stateOutDict))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Start time: 0.0\n", + "End time: 1.7976931348623157e+308\n", + "Current time: 0.0\n", + "Time step: 0.25\n", + "Time units: s\n" + ] + } + ], + "source": [ + "print(\"Start time:\", x.get_start_time())\n", + "print(\"End time:\", x.get_end_time())\n", + "print(\"Current time:\", x.get_current_time())\n", + "print(\"Time step:\", x.get_time_step())\n", + "print(\"Time units:\", x.get_time_units())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Show the input and output variables for the component (aside on [Standard Names](https://csdms.colorado.edu/wiki/CSDMS_Standard_Names)):" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "('plate_surface__temperature',)\n", + "('plate_surface__temperature',)\n" + ] + } + ], + "source": [ + "print(x.get_input_var_names())\n", + "print(x.get_output_var_names())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, get the identifier for the grid on which the temperature variable is defined:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Grid id: 0\n" + ] + } + ], + "source": [ + "grid_id = x.get_var_grid(\"plate_surface__temperature\")\n", + "print(\"Grid id:\", grid_id)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then get the grid attributes:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Grid type: uniform_rectilinear\n", + "Grid rank: 2\n", + "Grid shape: [6 8]\n", + "Grid spacing: [1. 1.]\n" + ] + } + ], + "source": [ + "print(\"Grid type:\", x.get_grid_type(grid_id))\n", + "\n", + "rank = x.get_grid_rank(grid_id)\n", + "print(\"Grid rank:\", rank)\n", + "\n", + "shape = np.ndarray(rank, dtype=int)\n", + "x.get_grid_shape(grid_id, shape)\n", + "print(\"Grid shape:\", shape)\n", + "\n", + "spacing = np.ndarray(rank, dtype=float)\n", + "x.get_grid_spacing(grid_id, spacing)\n", + "print(\"Grid spacing:\", spacing)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These commands are made somewhat un-Pythonic by the generic design of the BMI." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Through the model's BMI, zero out the initial temperature field, except for an impulse near the middle.\n", + "Note that *set_value* expects a one-dimensional array for input." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "temperature = np.zeros(shape)\n", + "temperature[3, 4] = 100.0\n", + "x.set_value(\"plate_surface__temperature\", temperature)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Check that the temperature field has been updated. Note that *get_value* expects a one-dimensional array to receive output." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 100. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0.]]\n" + ] + } + ], + "source": [ + "temperature_flat = np.empty_like(temperature).flatten()\n", + "x.get_value(\"plate_surface__temperature\", temperature_flat)\n", + "print(temperature_flat.reshape(shape))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now advance the model by a single time step:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "x.update()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "View the new state of the temperature field:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 12.5 0. 0. 0. ]\n", + " [ 0. 0. 0. 12.5 50. 12.5 0. 0. ]\n", + " [ 0. 0. 0. 0. 12.5 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. ]]\n" + ] + } + ], + "source": [ + "x.get_value(\"plate_surface__temperature\", temperature_flat)\n", + "print(temperature_flat.reshape(shape))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There's diffusion!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Generate an ensemble of models using get_state and set_state\n", + "We will demonstrate (some of) the use of get_state and set_state by generating an ensemble of models starting from the state of the one we just initialized, adding noise to each ensemble member and than running all of them forward in time. While running forward, we save the temperature at a single point of interest and plot a graph of all ensemble members to show the ensemble spread over time.\n", + "\n", + "First, we save the current state of the model." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "state_out = x.get_state()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "State is always returned as a string. In this case the string is formated as JSON, so Let's see what keys are in this state:" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dict_keys(['time', 'plate_surface__temperature'])\n" + ] + } + ], + "source": [ + "dictState = json.loads(state_out)\n", + "print(dictState.keys())" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dict_keys(['value', 'type', 'itemsize', 'nbytes'])\n" + ] + } + ], + "source": [ + "print(dictState['plate_surface__temperature'].keys())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we create an ensemble of BMI model objects. Each model gets initzialized in the same way as ```x``` was above." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "nEnsemble = 25\n", + "ensemble = []\n", + "\n", + "for ensembleMember in range(nEnsemble):\n", + " ensemble.append(BmiHeat())\n", + " ensemble[ensembleMember].initialize(\"heat.yaml\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, for each ensemble member, we take the state of x, add some noise, and set the state of the member" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "plate_surface_temp_from_x = np.array(dictState['plate_surface__temperature']['value']).reshape(shape)\n", + "\n", + "for ensembleMember in range(nEnsemble):\n", + " #for each ensemble, start with the same state as x.\n", + " ensembleMemberDictState = dictState\n", + " \n", + " noise = np.zeros_like(plate_surface_temp_from_x)\n", + " noise[2:5,3:6] = 2 * np.random.randn(3,3)\n", + " plate_surface_temp = plate_surface_temp_from_x + noise\n", + " \n", + " #this array we want to set in the state is a numpy array, which has to be transformed \n", + " #to a list otherwise we can not turn it into a json string.\n", + " ensembleMemberDictState['plate_surface__temperature']['value'] = plate_surface_temp.tolist()\n", + " \n", + " #Turn the state with noise added back into a json string\n", + " ensembleStateJSON = json.dumps(ensembleMemberDictState)\n", + " \n", + " #finally set the state of the ensemble member.\n", + " ensemble[ensembleMember].set_state(ensembleStateJSON)\n", + " \n", + "#Note that the above loop is needlessly verbose for educational purposes.\n", + "#This does use a lot of memmory. The code below is less readable, but uses less memmory.\n", + "#When using bigger models than this example model, this can make a big difference in \n", + "#performance.\n", + "\n", + "\n", + "plate_surface_temp_from_x = dictState['plate_surface__temperature']['value']\n", + "noise = np.zeros(shape)\n", + " \n", + "for ensembleMember in range(nEnsemble):\n", + "\n", + " noise[2:5,3:6] = 2 * np.random.randn(3,3)\n", + " \n", + "\n", + " dictState['plate_surface__temperature']['value'] = (plate_surface_temp_from_x + \n", + " noise).flatten().tolist()\n", + " ensemble[ensembleMember].set_state(json.dumps(dictState))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now run the entire ensemble forward in time. Every timestep, we are saving the temperature of one location of the plate of interest, for each ensemble member.\n", + "\n", + "Note that for bigger models, the for loop in this step can be run in parallel for all the models. " + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "loc_of_interest = [14]\n", + "\n", + "distant_time = 15.0\n", + "output = pd.DataFrame(columns = ['x'])\n", + "\n", + "outputValue = []\n", + "while x.get_current_time() < distant_time:\n", + " x.update()\n", + " x.get_value_at_indices('plate_surface__temperature',outputValue, loc_of_interest)\n", + " output.loc[x.get_current_time(),'x'] = outputValue[0]\n", + " \n", + " for ensembleMember in range(nEnsemble):\n", + " ensemble[ensembleMember].update()\n", + " ensemble[ensembleMember].get_value_at_indices('plate_surface__temperature',outputValue, loc_of_interest)\n", + " output.loc[ensemble[ensembleMember].get_current_time(),'ensemble' + str(ensembleMember)] = outputValue[0]\n", + " \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "View the final state of the temperatures of the entire ensemble at the locaiton of interest:" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'temperature')" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd4VEX7/u85W7JppCekQiAQkBBKQm+hhl4iSPlKr6+iwssPQQSVJoIgFhBQiqCoICiCdEW6dCnSQuhNAoEUEkKS3fv3x+4eUijxNdkEmM91zcXu2XPmzJzo3GfmeeZ5BElIJBKJRAIASlE3QCKRSCTFBykKEolEIlGRoiCRSCQSFSkKEolEIlGRoiCRSCQSFSkKEolEIlGRoiCRSCQSFSkKEolEIlGRoiCRSCQSFW1RN+Cf4unpydKlSxd1MyQSieSp4uDBg7dIej3pvKdOFEqXLo0DBw4UdTMkEonkqUIIcTE/58nlI4lEIpGoSFGQSCQSiYoUBYlEIpGoSFGQSCQSiYoUBYlEIpGoSFGQSCQSiYoUBYlEIpGoPHX7FJ4Hdu3ahf379yMsLAxVq1aFp6dnUTdJIpE8J0hRKEYkJSXhzTffxBdffJHjeEBAAKpWrYoaNWrgjTfegIuLSxG1UCKRPOvI5aNiwurVq/HCCy9g/vz5GDFiBC5duoRff/0V06dPR1RUFC5cuIDx48ejZs2aOHHiRFE3VyKRPKPImUIRc+PGDbz++utYvnw5wsPDsWrVKtSoUQMAEBgYiKZNm6rnbt++HV26dEGtWrXw1Vdf4cUXXyyqZkskkmcUOVMoQpKTkxEZGYlVq1Zh0qRJOHDggCoID6Nhw4Y4dOgQwsLC0LlzZ4waNQpZWVk2bLFEInnWkaJQhEyePBlXrlzBli1b8Pbbb0On0z3xGn9/f2zduhVDhgzBtGnT0LJlS9y6dcsGrZVIJM8DUhSKiLNnz+Ljjz9G7969Ua9evX90rZ2dHebMmYOFCxdi586daNWqFdLS0gqppRKJ5Hmi0ERBCLFQCBEvhPjrEb//nxDiqKXsFkJUKay2FEdGjhwJnU6H999//3+uo2/fvli+fDkOHjyIfv36gWQBtlAikTyPFOZM4SsALR/z+3kAjUiGA5gI4IvHnPtM8fvvv+Onn37CmDFj4Ofn96/qat++PaZMmYJly5Zh4sSJBdRCiUTyvCIK8+1SCFEawC8kw55wnhuAv0j6P6nOyMhIPs1JdoxGIyIiIpCYmIiTJ0/C3t7+X9dJEn369MGSJUuwfPlydOnSpQBaKpFIniWEEAdJRj7pvOLiktofwPqiboQtWLhwIY4cOYLly5cXiCAAgBACX3zxBeLi4tC7d2+UKVMGERERBVK3RCJ5vijymYIQojGAzwHUJ5nwiHMGARgEAEFBQREXL+Yrq1yxIykpCeXKlUOFChWwbds2CCEKtP74+HjUqFEDWVlZ2L9//79empJIJM8O+Z0pFKn3kRAiHMB8AB0eJQgAQPILkpEkI728nph3utgyefJk3Lp1CzNnzixwQQAAb29vrFmzBklJSYiJiUFmZmaB30MikTzbFJkoCCGCAPwIoCfJ2KJqh604f/48Pv74Y/Tp06dQl3bCw8OxaNEi7N27F1OmTCm0+0gkkmeTwnRJ/Q7AHwBChRBXhBD9hRBDhBBDLKe8A8ADwOdCiMNCiKfXepwPvvnmG2RlZWHChAmFfq8uXbqgR48emDhxIg4ePFjo95NIJM8OhWpTKAyeVu+jatWqwdHRETt37nzseQcOHMD777+PuLg4REVFoXnz5oiKioKzs/M/ut+dO3cQFhYGV1dXHDx4EAaD4d80XyKRPOU8FTaF54Vz587h8OHDiImJeeQ5u3btQqtWrVCjRg1s3boV3t7emD9/Ptq3bw93d3c0aNAA77//PlJTU/N1Tzc3NyxYsAAnTpzAuHHjCqorEonkGUeKgg346aefAACdOnXK89vOnTvRuHFj1K9fHwcPHsQHH3yACxcu4Ndff8WdO3ewZcsWjBw5Eunp6Xj77bdRt25dnD9/Pl/3bdmyJQYPHowZM2Zgx44dBdoniUTybCKXj2xA3bp1kZ6ejkOHDuU4fvToUURGRsLT0xNvvvkmBg4cCEdHx0fWs3HjRnTr1g2KomD58uU5wmo/irt376JKFXMEkSNHjsDJyenfdUYikTyVyOWjYsK1a9fwxx9/5Ml9kJmZiT59+sDNzQ1Hjx7FsGHD4OjoiPj4eIwZMwZdu3bFrFmzcPz4cTWmUXR0NPbv3w9fX19ER0fj448/fmK8IycnJ3z11Vc4f/48Ro4cWWj9lEgkzwgkn6oSERHBp4nZs2cTAE+cOJHj+IQJEwiAK1euJEleunSJr732Gg0GAwHQ29ubAAiABoOB3t7e9Pf355QpU3jt2jV26tSJANizZ0+mpaU9sR3/7//9PwLgr7/+Wij9lEgkxRsAB5iPMVYuHxUyzZo1w9WrV3Hy5En1mHXZ6MUXX8TkyZMxZcoULF68GCTRrVs3nDp1Ctn7aDAYYDKZkJGRAcBsRB4xYgTu3buHyZMno2XLllizZg202kdHLUlPT0flypUhhMDRo0elN5JE8pyR3+WjIn/z/6flaZop3Lp1ixqNhmPGjFGPZWRksFq1avT29uaBAwfo7OxMOzs7vvrqqzxx4gQbNmxIRVE4c+ZM/vnnn0xOTiZJmkwmfvvtt7Szs6NerycAurm5sX379gTAIUOG0GQyPbY9GzduJACOHz++UPstkUiKH8jnTKHIB/l/Wp4mUVi4cCEB8MCBA+qx8ePHEwBXrFjBRo0a0dnZmWfOnGFqaiobN25MRVH47bffPrLOv/76i+XLl6eiKKxUqRIB0NPTkwA4Y8aMJ7apW7dutLOzY2xsbIH0USKRPB1IUSgGtG3blkFBQeob/OHDh6nVatm9e3fOnDmTALhgwQKmpqaySZMmFELw66+/fmK9SUlJqk2hYcOGdHBwoKOjI4UQ/Omnnx577bVr11iiRAk2b978iTMLiUTy7CBFoYhJTk6mXq/nsGHDSOZcNtq9ezcNBgPbtm3L1NRUNmvWjEIILl68ON/1m0wmTp06lYqisEqVKnRxcaFOp6PBYOD+/fsfe+2sWbMIgN99992/6qNEInl6kKJQxHz//fcEwB07dpAkFy9eTABctmwZa9SoQXd3d167do1t27alEIKLFi36n+7z3XffEQBbtGhBLy8vKopCDw8PXrx48ZHXZGVlMTIykiVLlmRiYuL/dF+JRPJ0kV9RkPsUCokff/wRPj4+qFOnDgBzQLzSpUvj1KlT2L9/P+bMmYOTJ0/il19+wdSpU9GnT5//6T7dunXDhx9+iE2bNqFDhw7w9fVFQkICmjRpgvT09Ideo9FoMHfuXMTHx2Ps2LH/axclEsmzSH6UoziVp2GmkJaWRkdHRw4ePJgkefXqVSqKwv79+1Or1bJbt24kyRYtWtDHx4f37t37V/czmUwcOnSo6lkUGBhIAOzevftjr3vttdcohOC+ffv+1f0lEknxB3L5qOj4+eefCYAbN24kSc6YMYMAGBISQl9fXyYkJPDQoUMEwClTphTIPbOystixY0cKITh//nw6OjoSwGM9mRITE+nr68saNWrQaDQWSDskEknxRIpCEdK/f3+6uLjw/v37JMlq1aqxdOnSBMA1a9aQNLuGOjs7886dO/mqMzMzk5cuXXqsx1BaWhrr1KlDg8HAOXPmUAhBnU7HS5cuPfKar7/+mgA4f/78f9BDiUTytCFFoQipWLEi27RpQ9K8rwAAy5Qpw9DQUJpMJp49e5aKonDkyJGPrOPGjRtctWoVR48ezUaNGtHBwYEAGB4ezlmzZj3SQHzz5k2WK1eOnp6efOuttwiAgYGBj5wJmEwm1qtXj15eXvkWKIlE8vQhRaGISExMpBCCEyZMIEmOHj2aGo2GAPj++++TJF955RXq9XpevXo1z/X3799nlSpV1LhHsMQ+cnV1pZ+fH4OCggiA9vb27Nu3L/fs2ZNn9nDy5Ek6ODiwcePGbNq0KQGwc+fOj2zzoUOHKITg66+/XoBPQiKRFCekKBQRmzdvVu0JRqORgYGBDAkJoRCCly9f5o0bN2gwGNi/f/8818bHx9Pf3199u2/cuDFjYmLYtWtXdu3aldWrVycAhoaGsm3btqrdoFOnTkxJSclR16JFiwiA77zzjhpc73FLREOGDKFGo+GxY8cK/JlIJJKiR4pCETFp0iQC4J07d7h161YCoJeXF5s1a0aSHDt2LIUQPHXqVI7rDh48SBcXFwJgdHT0Q20H1vhHVu+itm3bcsSIEeoGtux7E0wmE19++WUqisKlS5dSo9FQo9Hkua+VW7du0d3dnVFRUXKns0TyDCJFoYho164dQ0NDSZIDBgygvb09AXDJkiVMSUmhm5sbO3XqlOOaJUuWUKvVEgDbt2/PtLQ0bt++nd988w0nT57MQYMGMTo6mrVr1+bMmTN548YNTpo0iY6OjtTpdBw4cCBLlChBHx8f7t69W603OTmZ5cqVo7+/Pz/55BMCoL+/P7Oysh7a9jlz5qgb7CQSybNFkYsCgIUA4gH89YjfBYBPAcQBOAqgen7qLc6iYDKZ6O3tzV69evHevXt0cXFhSEgInZycePfuXX700UcEwD179qjXvPnmm6rtoG3btrx8+XIem4KXlxcbVq3K/ytXjjqArq6uHD16NP/880927NiRAPjKK6+wbNmy1Ov1OeInHTp0iHq9nm3btmVUVBQBqPsncpOVlcVq1aoxICCAd+/eLfTnJZFIbEdxEIWGAKo/RhRaA1hvEYfaAPbmp97iLArnz58nAM6ePZsrVqxQDcJ9+vTh/fv36e/vz6ioKPX8bdu2EQCFEGzevDmPHz/O4OBgOjg4sEe7dpxQuzaXlCzJQxoNswAS4G0huMLTk1FCUK/Vsk+fPuzcubMqDNaBf+zYseoy0KeffkoAnDRpkjpz2b59+0P7sHPnTgLg22+/bZNnJpFIbEORi4K5DSj9GFGYB6B7tu+nAfg+qc7iLArWeEf79+9nx44dVRvB77//zpUrVxIA165dS9I8q4iIiKAQgjVr1uSuXbvo7e3NEs7O/EirZaZFBO4D3AZwkqJwSIkS/AbgXctv1zQaztBoGO7mpgrDwIED2a9fPwLg1KlT1Xt17NiROp2OU6dOJSyzDWuuhty8/PLL1Ov1jIuLs9mzk0gkhcvTIAq/AKif7ftvACKfVGdxFoXhw4fTzs6O169fp06nY6lSpViqVCkajUZ1Q1tmZiZJ8scff1SXhxYtWkRnZ2d6enpyohAkwG8Uhf+tUoVTxo3j9u3b+e2337J9+/Z88cUXWbNSJb6s0XANwAyAtwC21WrVhDsvv/wyu3btqtoySDIhIYG+vr4MCwtj8+bN1eWqh3H16lU6OjqyQ4cONnt2EomkcHkaRGHtQ0Qh4hHnDgJwAMCBoKCgwnliBUC9evVYp04dzp07V10WGjduHE0mE/39/dW9ApmZmSxXrhw1Gg2rV69OvV5PX19fvm6ZAXxtZ8dTJ0/SZDJx3bp1rFatGgEwICCAfn5+OewNZQEeBmgEOA5gsyZNCIAxMTFs3LgxtVqtGm7jl19+IQAOHz6cTk5OBMBvvvnmoX2ZMmVKjlAdEonk6eZpEIVnavkoIyODBoOBw4YNY+vWrdVsaGfOnOGxY8dy7BP44osv1EFdp9MxMDCQPS2CsEqj4dLFi9mtWzd6eHjkEAAAdHd356hRo3j06FFu2LCB06ZNo4tOx8WW69cCjAoPJwAOGzaM4eHhdHJy4sGDB0mSffv2paIo6qBvZ2fHCxcu5OlPeno6y5Yty4oVKzIjI8Omz1IikRQ8T4MotMllaN6XnzqLqygcPHiQALh06VK6uLjQ1dWVdevWJUl++OGHBMDLly/z7t27LFmyJA0Gg5r/oB3ATIC/CsFx2byRnlSCgoI4bdo0njt3jnVq1+Zgiw3iHMDmFkH5+OOPWapUKXp7e/Ps2bNMTExkQEAAK1SowFatWhHAIwPirVmzhgA4c+ZMWz9OiURSwBS5KAD4DsB1AJkArgDoD2AIgCGW3wWA2QDOAjiWH3sCi7EofP755wTAdevWqYP2vHnzSJJNmjRh5cqVSZKTJ09Wf9doNGwEMB3gHoC9LCk2/2mxs7Pj2rVrOXv2bNYEeBFgCsDmzs7U6XT89ttv6e7uzpCQEN66dYsbN25UvZWsy0iffPJJnj6ZTCa2bNmSJUqU4I0bN2z6PCUSScFS5KJQWKW4ikLv3r3p5eXF2bNnEwD1ej3v3LnDlJQU6nQ6jhw5kjdv3qSzszPd3Nzo7OxMFyF4C+BfAJtZ7AYPWy6ybmyrXbs2K1WqRIPB8NBze/XqxSNHjjDE0ZGxAG8DrGWZkaxcuZJ6vZ7R0dHMysrioEGDKITg22+/rS5jPczb6NSpU9RqtRwwYEARPFWJRFJQSFGwMRUqVGDbtm3Zo0cParVatmjRgiS5evVqAuBvv/3G4cOHUwihGqFHWOwArby88gzwoaGh1Ol06rnW4w4ODuzXrx83bNjAt99+m25ubjmu8/T05ObNmxnm7MxLAP8GWEmrZXh4uLqr+b333mNycjJLlSrFsmXLsm7dugTAWrVqPXQZacSIERRCPDH3s0QiKb5IUbAhd+7cIQBOnDiRAQEBBHJGRHV0dGRsbCz1ej2DgoKo1+tpEIJXAf4KUFGUHAO7VQQURVGFIXuxnh8cHMzx48eru5qzl/79+zNMr+ffAC8BDALYrl079urVi0IIrl+/nlu2bCEA9u7dW73Pw5aRkpKS6OPjwzp16si4SBLJU4oUBRuyadMm1b3TOijv3LmTJpOJwcHBbNeunertYy19LbOE5rkG8+yCEBYWxp49e3LGjBn87bffePz4cU6aNInlypXLca6npyfHjBmjhui2liZNmrC6RsPbAM8ALAnw3XffZXh4ON3d3Xn+/Hm+8sorFELwP//5j7rsFRsbm6ePCxcuJPBg34NEInm6kKJgQyZOnEgA/PLLL9WBNT09nadPnyYAfv7556xZsyY9PT2pKAo1QvC0EDz4CEFwdHTM4yZ6/fp1zpkzh8uXL+fevXu5ceNGvv7663R1dVWvHTlyJH19fXPUWb16dda2GJ6PAnSFOUVniRIlGBERwRs3bjAwMJAVKlRghQoVKIRgrVq18gTNMxqNrFWrFkuWLMmkpCRbPl6JRFIASFGwIW3btmVoaChfeeUVKorChg0bkqS6hr9r164cyz4dLLOErg8xFuv1ep49e1at+9q1a3zjjTcealy2t7dn+fLlWbZsWfVYq1at1B3L1hISEsLGMLurbgDo5e7OJUuWEDAHx7N6TPXt21cVpo8++ihPP/ft22e2hYwYYbNnK5FICgYpCjbCZDLRy8uLvXr1YqVKlQjLEg1JtmrViuXLl1c9kgBQEYJ7NRqeBah5yExh27ZtJM2hJl5//XUaDAZqNBr26dOHx44d45EjR/jzzz/zk08+4X//+19GR0cTAJ2cnFTR8ff3Z8OGDXPU7ePjwwEWMZoKMCoqSo3QunjxYvbs2ZNarZbdu3d/7DLSgAEDqNVqeeLECVs+ZolE8i+RomAjrJFRP/zwQ/Ut+7fffmNaWhrt7e35+uuvs1mzZuqbfgPLwPyfh8wSZsyYQdK858HOzo4ajYb9+vV7YmC6P/74g7Vq1SIA1X1Vr9ezZs2aOep3dHTkLMv9ewAcN24cGzZsSEdHR+7bt49eXl6sWrUqfX19qSgK69evn8cbKT4+nq6urmzWrJk0OkskTxFSFGyENTKqdalIURSmpqZyw4YNBMwJazQajfoWv0mn4w2AhlwzhPbt29NkMnHHjh3UaDRs0aKFuox07tw5Dh8+nKVLl2atWrXYo0cPjh07losWLeK2bduYlpZGk8nEpUuXquk8rcJQLdf+BwedjlsBpgGMALh8+XK6ubkxMjKSS5cuVb2RrOfPnj07T59nzZpFAFyxYoWtH7dEIvkfkaJgI4YPH06DwcBRo0apYbBJctiwYTQYDJw/f746wIZZ3tLH5poh+Pv78/bt24yPj6efnx9DQkKYmJjIHTt2sGPHjhRCqLMQnU5Hg8GQY++Cn58fFy1axKysLKampnLcuHHqb3Z2dgwLC8txP28heAFmV9VAvV71LHrzzTfZsWNHGgwGtmjRgoqi0MHBIY/ROzMzk+Hh4QwKCmJqampRPHaJRPIPkaJgI+rWrcu6deuyfv36BMDRo0eTNG9mi46OZkxMjJrY5mshmALQLZcorFixgkajkdHR0bSzs+PChQsZbglqZx38/fz8OGrUKA4aNIj169fPsWnNuscgPDycGzZsoMlk4o8//qjOTuzs7BgaGprjnlUBpgLcATA0OJgDBw6kEILff/89XVxcWKtWLTo5OVGj0bB58+Z5loq2b99OWJagJBJJ8UeKgg3IysqiwWDg0KFDqdfrCZhjH1ntDFOnTqWDgwMNBgMDYc598FEuQfD19WVmZiYnTZpEwLzb2M7OTv29RYsWXLNmTR4XUZPJxBs3bnDt2rWMiIhQB38AbNasGQ8fPsw//vhDbZfBYKC3t3eOe79kmbnMhjlERoUKFejn58eZM2cSALt166aeu3Dhwjz979GjB+3s7GQyHonkKUCKgg04c+YMAfCtt95S3+oTExM5Z84cAlDX3gFwuGUADs4lCuPHj+eWLVuoKApffPFFellCXrRv357nzp3LVzuMRiMXL17MkiVLqrYEnU7HRYsW8cyZM3R0dFRFI7vgAOA0S7tegtkNVa/Xs3379mzcuDGdnJxYvXp1arValihRglevXs1x36tXr9LJyYlt2rSRRmeJpJiTX1FQIPmfOXHiBAAgOTkZABAWFgYXFxds2bIFgYGB2L9/PwwGAwBznPC/AJy3XCuEgKIo6NChA3r06IFy5crhzJkzuHnzJmrWrImVK1ciODhYvdft27cxffp0jB49Gq+++ip69+6NmJgYREdHY9y4cYiKisKZM2cwduxYCCFAEn379sW8efMQFxcHd3d33L9/H1qtFkIItd4xAHYDmA/gy5EjMWbMGKxevRqNGzeG0WiEk5MTAODu3bt45ZVXzG8SFvz8/PDee+9h7dq1WL16dSE9ZYlEYlPyoxzFqRSnmYI1dEV0dDSFEBw+fDhJskyZMoyJiaG7u7s5IirMG8c+yPaGbmdnxxdffJEtWrSgwWBgy5YtVdtBYmJijvusX79e3ams1+vp7u7OUqVKsVKlSqxevToVRaGiKGzfvj03bNjAM2fOMCwsTLVHtG3bltevX1dzRlsTAFlLAMwpPf8EGOLvz+joaBoMBnUfQ4cOHdRzv/vuuxxty8jIYOXKlRkUFMS7d+/a7NlLJJJ/BuTyUeHTq1cv+vr6qoPtTz/9xNu3bxMwB6SDxUW1o2WJpmGupaPPPvtMHbQB8w7l7EtGd+/e5ZAhQwiAFSpU4IwZM/j7778zNjY2xwB84cIFjhkzRl16CgkJ4YwZM9SNbUIIVqpUiXv37lVtDNldVwGwlaWNc2EOnOfp6clq1aoxMjKSnp6eLF26NO3s7Ojh4cH4+Pgcz2HHjh0EwFGjRtns2Uskkn+GFAUbEBkZyTp16qgD682bN/nbb7+pb9dWr6AvASYC1OJBjoTQ0FD26NFD9UxSFIU7duxQ6969ezdDQkIIgHXq1KG7u3uOQRwAXVxcWLFiRTVXQ3p6Or/99lu1TdHR0ezXr5/qoeTh4cEVK1aoXklWG4S1TLEIQ3eYw19YxU2r1bJp06ZqO7t27ZrnWfTt25darZbHjx+35Z9AIpHkEykKhYzJZKKjoyObNGmivp2T5LRp09Q3ceugexXg8lwD+oQJE6jVatXIpgsWLFDrnjt3LoUQdHZ2VndCt2jRgiNGjOD777/P6dOnc/LkyXzttdfYunVrCiHo5OTEsWPH8s6dOzSZTOquaF9fXw4dOpRCCOr1erq6uubIEe3s7Kx+1sLsopoCMBRQ90hYN7M1a9ZMFZQff/wxx/O4efMm3d3d2ahRI2l0lkiKIVIUCpmLFy8SAGvWrKmGnibJbt26qWJgMBhYxfL23RsPYhDZ29tzzJgx6mD88ssvq/WePn2aGo2GQgh1Z3OLFi3U8BXWotfrWbFiRXbo0IFTp05l586dCYCurq6cNGkSk5OTeeTIETXyabdu3WhnZ0edTkd3d3e+8847al3Z8zn4A4wHeARgoKcng4ODGRQUxNDQUPr5+dHLy4sODg709vZmQkJCjmdiFZuvv/7apn8LiUTyZKQoFDLr16/PYbT9/vvvSZLlypVj+fLl1UH2LYso+GQzMPfp04d+fn7mMNoaDe/cuUPS7FpasWJFAuALL7yg5k1wc3PjiBEjeOTIEe7YsYMLFixQdx8HBwcTMIfInjdvHtu3b6/OVPbs2cO7d++yT58+BMCwsDDVXdXT01M9nru0tLR5NsC6detSURS2atWKQgg1AquiKOzVq1eOZ2INr+3t7c3bt2/b/G8ikUgejRSFQuajjz7KMZBevXqVSUlJBMCgoCDVW2gnwP25Bt0PPvhA/RwTE6PWac2XbN1LUK1aNS5YsOCxoSSMRiO/+eYbBgUFEQBbt27Nb775hsHBwdTr9Zw/fz5J8ptvvqG9vT39/f2p1Wqp0+no7e2t2h9yJ+j50CIMnSxLV9a6YREKq73kl19+ydGeQ4cOUVEUdeYkkUiKB8VCFAC0BHAaQByA0Q/5PQjA7wD+BHAUQOsn1VlcRGHAgAGq15Gvry9JcuvWrerSjrOzM90BZgF8L9ssoUaNGmzYsKG6ZHPp0iWS5JEjRyiEUI9//fXXedbmTSYTY2NjuXTpUg4bNoz169eng4MDg4OD2bNnT3bt2pUlSpSgoijs16+fahweMmQI79+/z+3bt9PZ2ZleXl5qqs+SJUuqXktWzyQA1AHcB/A2zKk8K1asSDc3N5YqVYpBQUF0cnKik5PTQ11o33jjDQoh+Mcff9jmjyGRSJ5IkYsCAA2AswDKANADOALghVznfAHgP5bPLwC48KR6i4so1K1bV/UOeumll0jmnT10s7xt18xlYLZ+btCgAUlzgDkPDw/1+KQfVyMTAAAgAElEQVRJk/Lcb9myZTn2FxgMBtatW5evvvoqO3XqlCMWklWsIiIi+Oqrr6pv99euXeP+/fvp4eFBV1dXCiGo0+kYGBiovvlnnzGUBZhkMT57u7vTwcGB1atXJwA2btyYgNndtX///jnampyczMDAQIaFhfH+/fuF/8eQSCRPpDiIQh0AG7N9fwvAW7nOmQdgVLbzdz+p3uIgCiaTiW5ubuqa/5w5c0iS//d//6cOyAC4xGK0VSzfNRoN+/btq/5+5MgRklSXZQCwZ8+eOWYIRqNRjXpau3Ztfvnllzx8+DAzMzNztMloNPLPP//kzJkz2a5dO3XG4eLiwkmTJtHBwYG+vr48fPgw//rrL5YsWZJOTk5qu3IHzLOW7hZhm2CxcwBgvXr1CICVK1dWvaPWrVuXoz1r1qwhAE6cOLGQ/xoSiSQ/FAdR6AxgfrbvPQHMynWOL4BjAK4AuAMg4kn1FgdRuH79uupJBIC7d+8mSVasWJEBAQHUarVUAN60CANgTn5Tu3Zt1V5QqVIlkuaEOtYBuF69ekxPT1fvc/fuXcbExBAwp8rM/tuT+Ouvv9TEO0IIDhgwgAEBAXR3d+fBgwd55swZlipVSh3UhRCPFIaFAI0Ao2DeRKfX6+nr68tSpUpRr9ezRIkS9PPzy2Ncfumll6jX63nq1KkCeOoSieTfUBxEoctDROGzXOf8F8AIPpgpnACgPKSuQQAOADgQFBRUaA8tv1g3qFnX4FNSUpiSkkIhBD09PakoCmtZ3rCz52Fu1aqV+nnz5s1MSEhQl2sCAgJy7BS+cOECq1SpQkVR+NFHH/1Pvv8mk4lfffWVOvCXLl2a/v7+dHV15b59+3jp0iWWK1cuhy0hMDAwjyg4AjwJ834LT5i9oaxeTw0aNCBg9kbq2bNnjvtfv36drq6ubNiwYZ4MbhKJxLYUB1HIz/LRcQCB2b6fA+D9uHqLw0whe/RTLy8vkuTOnTtzDKTjLUbm7LkTrDOLgIAAkmaDLGDebZw95/G+ffvo7e3NEiVKcP369f+6vUlJSaxRo4Zqi/Dz82OJEiW4e/duXrx4kQEBATmip5YoUUKdPViPVQGYDnANQDdXV3X5yDrDsEZiXbVqVY57W5MMffHFF/+6HxKJ5H+nQEUBgD2A0Pycm+0arWWQD8YDQ3OlXOesB9DH8rkigGsAxOPqLQ6i8Morr9DBwYEA2LBhQ5JU03Fay36LgVZ927YMmgC4ePFi3rt3Tx2IZ82apdadmprK0qVLs1SpUjx58mSBtdlkMnHAgAHqYO/t7U0nJydu376dp06dooeHB/V6vWqLyO2iCoBDLbOfYRZhA8z7NIKDg6nT6ejm5kYfHx/eunUrx32joqLo4uLCa9euFVh/JBLJP6PARAFAO5jdSs9bvlcFsDpflQOtAcTC7IX0tuXYBADtLZ9fALDLIhiHAbR4Up3FQRQaN26sDorWIHC9e/dWB/6SlsHzLTzIjFaqVCkC5h3HRqORc+fOJWB2U82eQMeam+H3338v8HabTCaOGDFCHeRdXV3p4ODArVu38uDBg3RycqJOpzPbRLLtcs5eVsEc8bW6RRCsHlF169ZVxSR3bKTY2Fja2dmxc+fOBd4niUSSPwpSFA4CcAHwZ7ZjR/NTeWGU4iAKPj4+apRR63JJWFgYvby8qNFo2MciCuHZBlPrzGLy5Mk0Go3q3oCBAweq9R4/fpw6nS7PTuGCJnsOZwcHBzo5OfHgwYPcunUr9Xq9urktu63BWtxhzu0cC9DJIgJlypQhYI7/ZPVoWr58eY57Tp48mYA5kqxEIrE9BSkKey3/SlEgmZCQQODBXoDz588zNTWVQgg6ODhQCMFvLUbZh71pp6WlcfXq1er3pKQkkua3+EaNGtHNzY03btwo9H6899576lKSVqulh4cHT58+zTVr1lBRFPV47kxtANjAYi9ZYgnEZ11G8vPzo1arpbu7Oz08PHj9+nX1fhkZGaxWrRq9vb158+bNQu+fRCLJSUGKwgIAPWDecVwOwGcA5uan8sIoRS0KVoOyoijU6/U0mUz8448/cgyacQB/sHwWQqhhrwMDA0mS4eHh6pKLlcWLFxMA582bZ7O+jB8/ngDU5SJ/f39euXKFX3/9tdoXa3ym3MLwjmU21AtmbySrh5M1X7ROp2OrVq1yeE0dOXKEOp2O3bp1s1kfJRKJmYIUBQcAkwHst5RJAAz5qbwwSlGLQvaw0+XKlSNJzp49O8fyCgG+iQcb1qyiMHLkSO7Zs0c917q/ISEhgZ6enqxdu7bNXTetS0l2dnYUQrBMmTJMSEhQl3seVRSAW2EOs10O5gRB3t7eBMDg4GDVg+nzzz/Pcb+JEycSAFesWGHTfkokzzsFIgowh6r4MD8V2aoUtSgMGzZMXWu3hrfo168fDQYDhRBsYRGFxtkGUKtr56VLl9Tdyz4+PmqdAwcOpEaj4eHDh23eH5PJpGZ3s9o9wsLCmJKSomaPe1TxhzmN5yGAdlZRdHent7c3dTodfXx8aDAYcmxey8jIYEREBL28vPJkcJNIJIVHQc4UtuSnIluVohaF6Oho1Uj82WefkSSrVq1KFxcXKorCty2iUCKXIDg7OzMuLk79PnXqVJLkrl27CID//e9/i6xPWVlZaj4G6xt+zZo1mZqaqkZIfZQ3UhtLf+doNOoSkqIorFatGgGzK25kZCQzMjLU+x07dox6vZ5dunQpsj5LJM8bBSkKMwCshnlHcoy15KfywihFLQqBgYGqG+aePXt47949ajQaddBcBfPuX+u6unXpKCYmRk2NqdFomJycTJPJxKpVq9Lf35/JyclF2q/09HQ2bdqUGo1GDa7XqlUrJiUlMSwsTLU9ZN/QZi3TLMLQTVFob2+vGp9DQ0PVdKPjxo3Lcb/333+fALhs2bIi6rFE8nxRkKKw6CFlYX4qL4xSlKKQnJxMWNbPAfDu3bvcv39/jgHyKh7EOwIehMJYu3atmj2tQ4cOJB8kvLfmPChqkpOTGRkZSYPBoHpX9evXj1evXlWzyeXOAAeY03juhjmiahmL6Lm6utLJyYmOjo4sWbIkhRCqDYU0R4atUaMGPTw8+PfffxdhryWS54MCE4XiVopSFPbt2/fAoOzuTpLqJjQA9LO8Mb+We9DUajlz5kz1+86dO0mSffr0oZOTE1NSUoqsT7mJj49n+fLlWaJECXUz3rhx43js2DFVDB+2lBQEMAHgAYD6bIJozULn5ubGMmXK5JgRHT9+nHq9nu3bt5d5nSWSQia/oqDgCQghFgkhFuYuT7ruWeTEiRPq5/LlywMADh06BL1eD0VRUMPy237LvzqdDgBQuXJlfPXVVwCAgIAA1K1bFykpKVi+fDm6desGJycnG/XgyXh5eWHjxo1wcHCAi4sLDAYDJk6ciG3btmHVqlUQQsBkMuW57hKAPgAiAMzUaqHVapGVlYXY2FhUrFgRqampOH/+PIYOHape88ILL+CDDz7A6tWr8eWXX9qqixKJ5DE8URQA/AJgraX8BqAEgLuF2ajiyokTJ6DRaAAANWvWBAAcPHgQGo0GJFEDQCbM8Tq0Wi30ej0AYPDgwTh27BgA4NVXX4UQAsuWLUNaWhr69+9fBD15PKVLl8batWuRlJSEUqVKQafTYejQoUhKSsKcOXMeed0aAB8BeCUrC50VBVqtFg4ODjh9+jQcHR3h7e2NJUuWYMmSJeo1b7zxBpo3b45hw4bh1KlThd85iUTyePIzncheYBaSIvNIKsrlo3bt2qneOT/99BMzMzPVjGUAuNHinmn9bjXKLl26VD1mTb9Zu3ZtVqxYsVgvm6xbt44ajYa1atWioihUFIVbtmzhsGHDcvQve9EB3AMwEebMbYrF+Gy1SQQEBNDR0ZGnT59W73P16lV6eHiwevXqMlObRFJIoLBsCgBCAcT90+sKqhSlKJQtW5bOzs4EwIsXLzI2NjbHgJgAcF6uQdLX15fR0dHqZ9K8lg6A06dPL7K+5Jcvv/ySANi8eXM19MXhw4fV/RYP2+0cBPP+hSMA7bMJZMWKFVX33KpVq+ZIGvTTTz8ReBBgUCKRFCwFJgoAUgAkZyuxAF7MT+WFUYpKFNLS0tRk9zqdjiaTKUcMo7IWI/OAbO6oADh8+HDVQGsNfjdixAhqtVqbxDgqCMaOHUsAbNeunep9dfLkSVaqVOmRwhANc7a2pVotRTbjdFBQkDrbev3113PcZ9CgQRRCFEqEWInkeafQZgpFXYpKFA4fPqwOeMHBwSTJadOmqd443SyiUAUPYgYByCEcK1eu5P379+nl5cVOnToVST/+F0wmE3v16kUAbNOmjepNdPz4cXUj38OWkqzxkd4wGKjRaKjX6+ng4EB7e3s19Hj2pDx3795l+fLlGRAQkCe1p0Qi+XfkVxTy4330W36OPetkN4JWrlwZAHDy5EnVG6cmgHswp5IDAJPJBAcHB/z888/qdY0bN8batWtx8+bNYmlgfhRCCHz55Zdo1qwZNm7ciCZNmuDOnTuIjo7G6tWrYWdnZ51V5mAizFmUpqanI5JEZmYmMjMz4eHhgStXrsDX1xd9+/bF5cuXAQCOjo749ttv8ffff2Pw4MEPrVMikRQujxQFIYRBCOEOwFMI4SaEcLeU0gD8bNXA4kJcXJz6uX79+gCAo0ePqgNXDQB/AsjKdk2DBg2wZs0aAEB4eDjc3NywYMEC+Pn5ITo62kYtLxj0ej1WrFiBChUqYN++fahevTquXLmCAQMG4IcffoAQIs81BPAygOsAlplMcLcIw5UrVxASEoL4+Hikp6eja9euyMjIAABERERg8uTJ+OGHHzB37lyb9lEikeDRy0cA3gBwHsB9mNNqnreUIwCG5mcaUhilqJaPevfureYW2Lt3L00mkxrOQQMwFeDHuZZPVqxYoX4eNWoUr1y5QkVR+NZbbxVJHwqCS5cu0c/Pj35+fixbtiwBsF69ejnyVucuETDnd96kKFSyHff29qaHh0ce+4LRaGTr1q2p1+t54MCBIuytRPLsgH+7fETyE5LBAP4fyTIkgy2lCslZBS1OxZ3sM4XKlSvj5s2buHvXvF3jBZjji+/Pdr6iKDh+/Lj6vVmzZliyZAlMJhP69etnm0YXAoGBgVi7di2Sk5NhMBjg5eWFXbt2YfPmzXjttdcAIM+s4SCA1wA0N5nwocEAIQQ0Gg2SkpKQnJyMMmXK4NNPP8V3330HwPzslixZAh8fH3Tp0gWJiYk27qVE8hyTH+UAEAbgJQC9rCU/1xVGKaqZgre3N4UQdHV1JUlu27ZNfePtZzGols/1Fmx1wdTpdExNTWVISAgbNWpUJO0vaDZs2ECNRsP69eurIbf79eunRlV9mOH5C8tz6m1vTyEENRqNun+hVKlSdHBw4F9//aXeY/fu3dRqtezYsWOx3s8hkTwNoAANze/CnG3tMwCNAUwD0L4gham4k5KSgvj4eJBEmTJlAOQ0PNcAkAQgLts1VapUwalTp6AoCurVq4dDhw4hLi7uqZ4lZCc6Ohrz5s3Dzp070bhxY2i1WixcuBDh4eGoWLHiQ43EQwHsAjD73j1UVxQYjUb8/fffCA4OxuXLl6HX6/Hiiy8iOTkZAFCnTh1MmzYNq1atwscff2zbDkokzyn5CXPRGUBTAH+T7AugCgC7/FQuhGgphDgthIgTQox+xDkvCSFOCCGOCyG+zXfLbcjZs2fVzxEREQDMnkdWagA4AIDZlk0qVKgAkjCZTGjWrBnWrVsHjUaDjh072qrZhU7//v0xbtw4rF27FjExMRBCYPr06ejWrRu8vLzynJ8B4EUAdwCsMBrhaTl+/vx5+Pj4gCTOnDmDfv36qaIybNgwdOzYEW+++Sb++OMPW3VNInl+edJUAsA+y78HYY57JAAcz8d1GgBnAZQBoIfZQP1CrnPKwey042b57v2keoti+eiHH35Ql0CWLl1KkoyKiiJgzjiWAXBKrqWShg0bqp/37NnDyMhI1q9f3+ZtL2xMJhP79u1LAOzatava5ylTpqib9nKXCID3AG4Vgtps+zocHBwYGBiYZ7f3nTt3GBwczICAAN68ebMIeyuRPL2gAHc0fw7AFcAQAGcsg/iifFxXB8DGbN/fAvBWrnOmARiQn4ZaS1GIgjUhDACeP3+eJNV8xDUt6+Qx2QY9Ozs76vV6CiHo7OzM69evUwjBCRMm2LzttiAjI4MtW7akRqNhx44dVZvC5MmTH7rbGQBftjy3eZZ8E4qiqCFEgoODqdFouHnzZvUeBw4coJ2dHZs0acLMzMwi7K1E8nRSIKJgmRUEZvteGkB4vio2LzvNz/a9J4BZuc5ZZRGGXQD2AGj5pHqLQhT69eunBoQzGo1MTU1VB7dXLYNbYLYBz/q2a80V8P3336szhmeVlJQUVq9enfb29mzcuLHZVVej4ZgxYx7pqvqh5dm9bjFUCyHo4+NDAPT396erqytjY2PVeyxevJgA+MYbbxRhTyWSp5P8isJjbQqWilZl+36B5NHHXZONvLuZzINBdrQwLyFFAegOYL4QwjVPRUIMEkIcEEIcuHnzZj5vX3DExcWBJFxcXKAoCs6cOaP+VgPA3wCuW0JqA4CDgwMAICMjA02bNsWmTZvg6uqKyMhIG7fcdjg5OWHt2rUoWbIkjh49isjISBiNRkyfPh1Dhgx56DWjAWwAMD0tDa0tu6Jv3LgBX19f3Lp1CyTRrl071SW1V69eGDZsGD755BMsWrTIdp2TSJ4j8mNo3iOEqPHk0/JwBUBgtu8BAK495JyfSWaSPA/gNMwikQOSX5CMJBn5MANmYRMbGwuSCAw0dye359F+AMhmZL59+7b62SoKzZo1U3MxPKuULFkSGzZsAADcunULoaGhyMzMxJIlS9CpU6c85xsBdIX5j770/n1UUsz/OV6/fh2Ojo7QarWIi4tD9+7dYTQaAQAffvghmjZtiiFDhmDv3r026plE8vyQH1FoDLMwnBVCHBVCHBNC5Ge2sB9AOSFEsBBCD6AbgNW5zlllqR9CCE8A5WHePV1sSE1Nxd9//w3A7FEEPMjAZg+gAswW+Kwsc4ALIQRu3rwJg8EAX19fCCFw5coVtGjRoghab3vKly+PX375BfHx8dBqtfDz80N6ejp+//13NGjQIM/5yQDaAEgHsNpkglXyb9++jbS0NAQEBGDDhg148803AZiTFy1btgz+/v7o1KkTrl3L/Z4hkUj+DfkRhVYwexA1AdAOQFvLv4+FZBbMrukbAZwEsJzkcSHEBCGEdZ/DRgAJQogTAH4HMJJkwj/vRuFx7twDjbIu/xw6dAgAUBHmB/hXtvMdHR0BmAPiNWnSBJs3bwYANG/e3BbNLRbUrl0bK1euxOnTp+Hn5wcXFxekpKTg9OnTCA8Pz3P+JZg3vpQE8DMAg+U4SVy6dAkhISH46KOP1CUjDw8P/Pzzz0hOTkZMTAzS09Nt1DOJ5DkgP4YHAPUB9LV89gIQnJ/rCqPY2tC8cuVK1TC6d+9ekmRwcDABsKfFUFoxm/HUz89P/bxo0SK2bt2a5cuXt2mbiwvWjHONGjWivb09NRoN/fz8WLp06YcanmMsz/MHjYbCcszFxYUAWKZMGep0Om7fvl2t3/q36dmzp9zxLJE8ARTwjuZRMLuUAoAOwDcFokhPAdljHoWHh8NkMuHSpUsAgEowRwu8bPdgL5/52ZupX78+tm7d+twsHeWmR48e+PTTT7Ft2zY0bNgQQgjEx8cjKyvroZvbfgQwEkBnoxEfWvJbJyUlwcPDA+fOnYOXlxc6dOigbhyMiYnB+PHj8fXXX2PChAk27JlE8uySn+WjTjDP7lMBgOQ1AM6F2ajiRFxcHIQQsLe3h8FgwKVLl1SjZxjMRtL7JpN6/q1bt6DRaFCmTBlcvnwZaWlpz60oAMBrr72GcePGYePGjWjdujVMJhOuX78OBwcHODk55Tl/OoB5AEZkZOC/BvNCUkJCAjw8PHD79m0IIdCqVStcv34dADBu3Dj06dMH7733HhYvXmzDnkkkzyb5EYUMy9TDvHFBCMfCbVLx4syZMyAJb29vADk9jyrBnFQnMzMTgNkImpmZCUVRUKtWLWzatAlarRZRUVG2b3gxYvz48Rg8eDBWr16NTp06wWg04sqVK/D09ISdXd6IKa/CbFv4MD0dPS2/JyQkwGAwgCTi4+PRpk0bpKSkQAiBefPmoWnTphgwYAB+++25y/8kkRQo+RGF5UKIeQBchRADAfwK4MvCbVbxITY2FgAQHBwM4IHnkSPMO/myG5lLlCgBwCwSERER2LRpE+rWrQtn5+dmYvVQhBCYPXs2unfvjpUrV6Jjx44wGo24fPkyfH19odVqc5xvhNlVbReA+ffvo5VlKSkxMRGZmZlwdXXFkSNH8NJLLyEzMxN6vR4rV65EhQoVEBMTg7/++itPGyQSSf54oiiQnA5gBYCVMLuMvkPys8JuWHHg3r17qsuj1Wtm/35z1oQXLOccz3a+3jJ4AUBISAgOHTr0XC8dZUej0WDx4sXo0KEDVq1ahXbt2qnCEBgYCEXJ+Z9iOsxrlqcBLMvIQE3LHo+7d+8iISEBQUFB2LBhA4YMGaJuLFy7di0cHR3Rpk0bdXlJIpH8M/IzUwCAYwB2ANhu+fxckN0dtW7dugCAI0eOADDbEwDgVLbBzJp0BzAvdwDPlyvqk9DpdFi2bBmio6OxZs0atGzZEkajEZcuXUJgYGCe5DyJAFoCSACwxmhEiOV4VlaW6qq6cOFCvPvuuwCAoKAgrF27FgkJCWjTpo0aglsikeSf/HgfDQCwD0AMzPGM9gghno2kAE8gu+dRo0aNAAAXLlwAYLYn3AMQl83IfPfuXdjb2yM0NBQ7d+6Em5ubGmpbYsbOzg4//vgjGjVqhE2bNqFp06YwGo24evUq/Pzypv6+BiAa5pgpGwH4wrwHRFEUxMXFoXz58pg4cSI++ugjAEC1atXwww8/4NixY2jXrh3u3btnw95JJM8AT/JZhXkG75HtuweA0/nxdy2MYst9CtOnT1cjeJJkQkKC6lO/HuDBbD72BoOBAGhvb8/u3bvT39+fXbp0sVlbnzaSk5NZu3ZtarVaNmvWTM1QZ83ElrtEAkwGeAKgV7ZotADUDHdffPGFWv93331HIQTbtGnDjIyMIuypRFI8QEHtU4A5PlFKtu8pAC7/WzF6GrDOFFxdzTH6Tp8+rf5m9TyyYt3JfO/ePQQGBuLq1avSnvAYnJ2dsX79elSuXBnbt29HkyZNkJmZidu3b6ueXtk5AHM4jCCYPR3cAdy/fx/29vY4efIkypcvj8GDB+P7778HAHTr1g1z587F2rVr0atXL9WNWCKRPJ78iMJVAHuFEO9ZNrLtARAnhPivEOK/hdu8osUaDTUgIADAA88jF5gj/WX3cTFlW0ZKS0sDIO0JT8LV1RW//vorwsLCsGPHDjRq1AgZGRlITEx86Oa2HTAbn8sD2Azz3+HevXtwdHREbGwsypQpg549e2LNmjUAgEGDBuGDDz7A999/j1dffTXHxkKJRPJw8iMKZ2EOXGf9P+pnANdh3sD2TPtaWnfOhoaGAgB2794N4IHn0Yls5yYmJkKv10MIgXPnzqFcuXIoVaqUDVv7dOLu7o5ff/0V4eHh2LVrFxo0aICMjAykpKSoM7TsbIF5N2UlmG0MzjAHLXR0dMTZs2cREBCALl26YMuWLQCAUaNGYfTo0Zg3bx7GjBljw55JJE8p+VljKk7FVjaF9PR0CiEIgFOnTiVJvvDCCwTAAZYYPWUVRU0OA4Bubm4MDQ2lt7c3e/XqZZN2PivcuXOHNWrUoFarZb169QiAjo6OdHJyeqiNoT3MaVB3AHS0HHNwcKAQggEBAXR0dOTWrVtJmlOGDhkyhAA4fvz4Iu6pRFI0oABjH0UKIX4SQhyyhM4+ms/Q2U8158+fV5cbrJ5HFy9eBGB2R00FcM6yZGRvbw/AvGmtYsWKiI+PR61atWze5qcZV1dXbN68GREREdizZw8aNWqE1NRUKIry0F3PqwH0gDnn6zoATjAv2xkMBly/fh1ubm5o1aoVtmzZAiEEZs2ahd69e+Pdd9/Fe++9Z9O+SSRPE/lZPloKYBGAF2EOmW0tzzTZ3VEjIyNx//59pKamAnhgZLaup+l0OgBml1TrrubatWvbsLXPBi4uLti4cSNq1KiBHTt2oEmTJkhOToa9vX2eXc+AeUdlDwB1kdPGoNfrcePGDXh5eaFNmzbYvHkzNBoNFixYgL59+2L8+PF45513pI1BInkI+RGFmyRXkzxP8qK1FHrLihirKBgMBmg0Gpw9e1b9LbfnUVpamrojNz09HQaDAZUrV7Zha58dXFxcsGnTJjRq1AhbtmxBkyZNkJiYCBcXlzyb2wBgOcxvK9Vgtjd4wCwMWq0WV69ehbe3N9q1a4eNGzdCo9Fg/vz5GDBgACZOnIhx48ZJYZBIcpEfUXhXCDFfCNFdCBFjLYXesiLGKgpWLxjrTmZ3mDdQZReFzMxMddC6cOECIiIi1NmD5J/j7OyMdevWoUOHDtiyZQuioqKQkJAAT09PAMgjDqsBdIA56dE2mJP13Lt3DzqdDleuXIGPjw86dOiAdevWQVEUzJs3DwMHDsTkyZPx9ttvS2GQSLKRH1HoC6AqzBEHrEtHbQuzUcUBq+eRNRDe+vXrAZhnCQBwItfA5OjoiNDQUBw5ckQuHRUABoMBK1asQM+ePbF161Y0aNAAN2/eRMmSJc3GsFyxkjYCaA2gFMzCEADzrE2r1eLy5cvw8fFBp06d8OOPP0JRFMydOxeDBw/GlClTMHLkSCkMEomVJ1miARzLj8XaVsVW3kfWnbVDhw4lSVaoUIEAOMTieVTK4nmk0WgIgK6urmzdujUBcPny5TZp4/OA0Wjka6+9RgCsV68eFUWhj49PjmefvdQBmAjwPHqKqWkAACAASURBVMDylmM6nY4AWLp0aSqKou58NhqNHDp0KAGwb9++zMzMLOLeSiSFBwpwR/MeIcQLTz7t2SEzMxM3btwAANSpUwcAcPXqVQDmmUISgIsWzyPrUkZiYqLqhSQ9jwoORVHwySef4J133sGuXbsQERGBpKQkeHp6wmg0QmOJnmrlD5iTiTvAHHq7Jsx/T41GgwsXLiAkJASDBg3C+++/DyEEPv30U7z77rtYtGgROnfuLPM9SyRPUg0AJwFkwBwD6SjMUVKP5kdxCqPYYqYQGxurvnleu3aNJNW9CFsA7sr1duro6EgAjI6OZsmSJWW+4ELi888/p6IorFChAt3c3Oji4kKtVkutVptnxlAW4BmAqQDbWo4pltld+fLlCYDDhg2j0WgkSX722WcUQrBRo0ZMTEws4p5KJAUPCnCm0BJAOQAt8MCekC+XVCFESyHEaSFEnBBi9GPO6yyEoBAiMj/1FjZWI7OiKPD19UVCQoK65hyGnEZmAPDw8ICiKDh79ixq1679UC8Zyb/nP//5D1atWoVLly7B3t4eDg4O0Ol0cHBwyDNjOAuzq+oJAD8B6IcHoUhiY2MREhKCjz/+GL169UJmZiaGDh2KpUuXYteuXYiKilJnihLJ80Z+kuxchDnUTxPL57T8XCeE0ACYDaAVzJEhuj9sGUoI4QzgdQB7/1nTCw+rKFgzpm3duhUA4GUpuUXBzs4O5cqVQ1xcnFw6KmTatWuHrVu3IisrC2lpafDz80NqaupDg+jdBBAFcwC9BQDGZvstLi4OwcHBWLp0Kdq0aYP/396dx1Vdpo0f/9zsIIiCCoIoCiruhJiik1ZmbpXpOI9WtoyZbfZMv5qnfdKpZn7NNM3U2PLY1K/Fdi1LLbXUct8VE1Nc2AUVEUEQ2c79++M+5+sBQcmAA3K9X6/vi3O+5+vhOnUO17m36y4oKOCWW25hyZIlJCcnM2zYMGvXPSFakrr8cZ8NPA48aT/lCXxYh+e+EjiktU7RWpcBn2JmDlb3PPB3zGZbTUL1QnjffPMNcG7m0b5qM19yc3OtayUpNLxBgwaxefNmQkNDyczMpH///uTk5NRYa6oY06x9H/NGewdw7I+Xnp5OeHg4q1evZujQoaSmpjJmzBhWr15NQUEBQ4YMsb4QCNFS1KX7aCKmOGUxgNY6m7oVwgunaontLPs5i1LqCiBCa720TtE2kp07dwLnCuFt3boVOLfbmmM6qmNa5KlTp/D09EQpRXx8k+gBu+x17dqVDRs2MHjwYHbt2sXgwYNJT0+vMTFUAHcBf8Z0I60E2mG6k44ePUrbtm3JyMhg8ODBbNq0iSFDhrBlyxZCQ0O5/vrreffddxvzpQnhUnVJCmX2QQoNoJRqVcfnrqlj3ZoMrpRyA/4FPHrRJ1JqplJqu1Jqe25ubh1//aXbv38/gLVrWmamyW19gJNAlr02v81ms/ZROHXqFH379rW6nETDCw4OZuXKlcyYMYMtW7YwYMAAjhw5QseOHXF3dz9vbGcOMAWIx2wl2AeorKzk1KlTlJeX4+bmxjXXXMMnn3xCt27d2LhxIyNGjGD69Ok89dRTVcqjC3HZuthINPBHYB6QAtyDmfX3UB3+XQKwwun+k8CTTvcDgRNAmv04i9l9Mf5Cz9vQs4/Onj1rzWBZt26d1vrczKO1oNdUm+USHh6ulVK6bdu2esaMGQ0am6iZzWbTc+fO1e7u7rpLly46MDBQBwYG6uDg4Fp3cTti38ltfLWZSV26dNGAnjNnjq6srNRlZWX63nvv1YCePHmyLi4udvXLFeKSUI+zj9pjao99AfQEnsUsGL2YbUB3pVRXpZQXMBVTkcCRjAq01u201pFa60jM5j03aa231+G5G0xKSop1Oz4+npycHGvmUfWaR2CK4UVHR5Ofny/jCS6ilGLWrFl89913nD59Gq01bdu2JT8/31qR7mw7MAgzx3oxZsDM0QpIT08nMjKSOXPmMGnSJM6cOcObb77Jyy+/zBdffMHQoUOrvEeEuNzUJSmM0lp/r7X+H631H7XW32NmFF2Q1roCmIWpQLAP+FxrvVcp9ZxS6qZfF3bDccw48fLywsfHh1WrVgGm3lEQVTfWATh58iQhISGADDK72rXXXsu2bduIiIggPT2dvn37kpqaSnR09HnXZgPDMQX1XsRMWw20P5aWlkZERARLly5l0KBB7N27l0ceeYRvv/2WjIwMBg4caJU9EeJyU2tSUErdr5TaA/R03kdBKZWKWcR2UVrrb7XWPbTWUVrrv9jPPau1XlzDtVe7upUA5/ZhdhTC+/7774FzM4/2V5sPX1hYCIC/vz+9e7eohd9NUrdu3di8eTNTpkzhp59+olevXqSlpdGxY0d8fHyqXFsC3AI8jNn/eTvQ3/7YkSNH8PPz48SJEwwZMoTPP/+cMWPGsH37diIjIxk/fjzPPfecjDOIy86FWgofY2bzLabqPgoDtdbTGiE2l0hMTASge/fuAGzZYpZPOJLCT05/BBybv+Tm5jJo0KDzFlAJ1/D39+fjjz/m9ddf5/DhwwQFBVFaWopSis6dO593/auY9Qy+mD7MOzHdSUVFRZw6dYq2bdsyZcoUHn30USIiItiwYQPTpk1j9uzZTJgwgVOnTjXq6xOiIdWaFOx9/mla61u00z4KWuuTjRlgY3NMR42NjQUgIyMDMEkhFziuz1XTDAoKAuDw4cPSddTEKKV44IEH2LBhA35+fhQWFhIWFkZGRgYxMTHnXb8RsyfDJuA94C3Axz7wlpWVRUREBP/85z8ZMWIEx48f5/3332fu3LksX76cuLg468uDEM1dXcYUWhRHEhg+fDgFBQWUlJQANQ8yBwQEEBQUREVFhSSFJio+Pp6dO3cyduxYDh8+TM+ePTlw4AAdO3bE39+/yrW5mFouf8VMs9sODLA/lpWVRdu2bdm9ezexsbEsWLCAWbNmsXbtWmw2G7/5zW948cUXpTtJNHuSFJw4J4GhQ4daXUlg6nRUTwqlpaXWxi+SFJqutm3b8tVXX/Gvf/2L1NRUAgMDqaio4OzZs3Tr1q3KtZXA08B1QBtM7ZWHAbQmPz+fkpIS/P39mTJlCvfccw/9+/cnMTGRiRMn8uSTT3L99deTk5PT2C9RiHojScGJo7yFm5sbHTp0YPt2M+4dhvkDUX1jnezsbKufumPHjo0crfgl3NzcePjhh9m2bRvh4eHk5ubSpUsXUlJS6NGjx3mb9qzCDDovx6ywXIbZ0U1rbS2Qe/vtt4mPjyctLY3PPvuMt99+m40bN9K/f3+rNIoQzY0kBSeOmUdt27ZFKcW6deuAmmseeXh4UF5ezokTJ6SV0Iz079+fbdu28Yc//IHDhw/TsWNHMjMz8fPzIyIiosq1ecDNwH3AVZgpdxPtj+Xk5ODt7U1OTg7x8fG88MIL3HHHHezYsYOwsDBuuOEGZs6cac1OE6K5kKTgxDHI7OhScLQUHElhj728BZybspqXl8eVV17ZeEGKX83Hx4dXXnmF5cuXo7WmrKyM4OBgMjMz6dGjx3nXzwMGYgp5fQl8gqmdVFpaSkFBAR07duTZZ58lISGByspKtmzZwmOPPcY777xDv379WLlyZaO+PiF+DUkKThyF7wYMGMDZs2er7LZ2HFOTwyEoKMiqreOokSSal9GjR5OUlMStt95Keno6oaGhpKen07p1a8LCwqpcux8YjCm/PQkzvjTZ/lhWVha+vr4kJycTFxfHK6+8wl/+8hc2bNiAr68vo0aN4v777+f06dON+vqEuBSSFJw4VjMPGzaMpKQk63xNM488PDysQeYrrriikSIU9S04OJgPPviApUuX4u7uTllZGf7+/mRnZ9OzZ88qa08qgL8AcUAGsMB+dABKSkooKioiKCiIJ598kmHDhuHr68uuXbt49NFHmTdvHv379+e7775zxcsUos4kKdhprcnLywPMvsy7du2yHqtp5lFubi7e3t5ERUXRpk2bxgtUNIjx48ezd+9eZsyYQXZ2Nu3btycjIwNPT8/zynHvBYYAT2BWc+4DZmLKAh87dgwvLy9+/vln4uLiePbZZ/nzn//MunXr8PLyYvTo0dxyyy0cPXq0sV+iEHUiScHu6NGjVNrHDKKioqykEI6pieNc88jNzY3s7GyKi4uJi4tr9FhFwwgMDOStt97i+++/p3Xr1pSUlBAcHEx6ejoRERHWCnYwU1f/hlnHsBsz7rAJiAXKysqsVsM//vEPevfuzcmTJ9m9ezdz5szhyy+/JCYmhjfffFPWNYgmR5KCnWPmkb+/Px4eHuzYsQM4N8jsPB3V0W2Un58v4wmXoeuuu46kpCTmzJnDiRMn8PHxobCwkLKyMmvjJYdk4FpgGhCJWfD2CmYXqhMnTqCUoqioiJtuuolbb72V3//+9+zZs4eBAwfywAMPnLceRghXk6Rg50gCnTp1orKykt27dwNOM4+cyluEhoZat6WlcHny8fFh9uzZJCUlMWLECAoKCmjXrh2HDh0iICCgynsA4CMgBtNieAiTLO4C0JqTJ0/i6+vL0qVL6dGjBx9++CFff/018+fPJyUlhbi4OGbOnMmxY8ca+VUKcT5JCnabN28GzMyjAwcOUFpaCpikcAyz45pDQEAAXl5mp19JCpe36Oholi1bxsKFC/H29qayshI/Pz+OHj1KWFhYlcqrp4AHMeMN6cC7mJbDcMxAdHl5Oa1ateL5558nJiYGpRTJyck8/PDDvPvuu3Tv3p2XXnrJeu8J4QqSFOwcLYPqg8w1zTwqKioiICCALl26EBwc3HhBCpdQSvHb3/6W5ORknn/+eU6fPo27uztnzpyxSmU4b/25DbPt4C2Y9QxrMDtUdcPsvwFQXFzMtGnTGDduHP/1X/9FUlISw4cP57HHHqNPnz4sWrTI2txJiMYkScHOsSahLjOPMjIyKC8vl1ZCC+Pn58czzzzDoUOHuOuuuzh16hS+vr5WraPIyMgq13+K2arwaUyhvX2Y8YYOmD29lVL89NNPJCQk8MQTT/C3v/2N5cuX4+XlxaRJkxg2bBhr1qxpzJcohCQFgPLycs6cOQNATEyMtYgtAmhN1aQQFBREfn4+hYWFMsjcQjnqHiUmJjJ06FBKSkpo1aoVWVlZeHp6Wjvxgdl4/K9Ad0xJ7gcxm52/ALTWmjNnzuDh4cGyZcvo168fn376KYsXL2bevHmkp6dz9dVXM3bs2CpfVIRoSJIUMNsvgtmC09vb20oKjn3UnJOC8yYt0lJo2QYMGMDKlStZuXIlffr0oaKiAh8fH44fP46fn1+VrsWjwL1AL8yuVU8DqZi1Dl4VFZSWluLh4cH8+fPp3bs3e/fuZe3atfz9739ny5YtxMXFMXXqVPbt2+eKlypaEEkKnKtx1KFDB3bu3MnZs2cBp+moTtd26NDBui1JQQCMHDmSTZs2sXjxYiIjI9Fa4+HhQV5eHn5+flUWNx4CbsWsb1gP/F8gDZMcfMrLrbUyc+fOpXfv3mRkZLBu3TqeeeYZli5dSp8+fZgyZQp79uxp7JcpWghJCsCGDRsA6N27N+vXr7fO98F8w3OeeVReXo6vry/h4eFVuglEy6aU4sYbbyQxMZFPPvmE8PBwwCx0LCgowNfXl4CAAOv6n4CbMAPSWzHJIR2YAwSUl6O1xmaz8cYbbxAbG0tubi5r1qzhiSeeYNmyZfTv35+JEydaRRyFqC8NmhSUUmOUUslKqUNKqSdqePwRpdTPSqmflFKrlFJdanqehuZoKQwbNowNGzbg6ekJ1Dzz6NixY7i5ucl4gqiRm5sbU6dOJSkpiYULFxIdHY3W2lrE5uXlRWBgoHX9ZuAGTBXW1cBsTMvhb0BIRQU2mw2bzcZ//vMfBg8eTGZmJkuXLmX27Nn8+OOPDBw4kLFjx7Jq1SqZrSTqh7bvQ1vfB+AOHMbMxPPCVAPoXe2aawA/++37gc8u9rwDBw7U9a1du3Ya0CtWrNDBwcHa3d1dA7oQ9KugsR9t27bVnp6eGtBz5syp9zjE5cdms+lvvvlGDx06VAPa29tbu7u7a6WUbtOmjfXechx9QH8EugJ0Gej5oGPtjymlrPfmqFGj9JdffqlfeOEFHRISogEdGxur58+fr8vKylz9skUTBGzXdfjb3ZAthSuBQ1rrFK11GWaG3gTnC7TWP2itz9jvbgY6NWA8tTp16hRgBprz8vKorKykM6ZUgXNLISYmhvLyckDKZYu6UUoxbtw41q9fz5o1axgzZgyVlZUopSguLgbMjDaHvcBtQBTwGuYDswvTihinNTb7mMOPP/7IpEmT+PTTT3nuued44403KC0t5fbbb6dbt2689NJL1poIIX6JhkwK4Zh9SRyy7Odqczdm18NGVVRUREVFBW5ubtZ2nHBukNk5KTj6iUEGmcUvo5Ri+PDhfPXVVxw4cID7778fDw8PAGs6tL+/v7UILh14BDMt+o9ANLAUOGi/39r+5WT//v3ce++9PP3004wbN463336bHj168NhjjxEeHs706dOt7lEh6qIhk4Kq4VyNnZ5KqWlAPPBSLY/PVEptV0ptz83NrccQz+221qZNGzZu3Iivry9wbjqq88wjd3d3lFKEhITInsziknXv3p3XXnuNrKwsXnzxRWtGW1lZmWm+u7lZFVkLgJcxfbBTMd+sXrL/fA+Iq6gATGv35ZdfZsaMGfj4+PDaa69x55138vnnnzNo0CAGDx7MBx98QElJSWO/XNHMNGRSyMJ80XHoBGRXv0gpdR1m2vZNWusai75ord/SWsdrreMd22DWF8eK0aioqPNmHuUA+U7Xnjx5Ei8vLwYOHFilrIEQlyIoKIjHH3+clJQUlixZwnXXXQeAzWazpqY6vqRUAJ8BVwN9gXcwO8BtwXQvzdIaRyfU8uXLmTVrFkuWLOG+++7j2WefpbCwkDvvvJOwsDAefPBBduzYIQPTomZ1GXi4lAPwwCze7Mq5geY+1a65AjMY3b2uz1vfA83jxo3TgJ4xY0aVAb8toL93ut++fXsdERGhlVL6mWeeqdcYhHBISUnRjz/+uG7fvr0GrIkNzoPMjsMf9H2gt4HWoM+C/gz0GNBu1QawR4wYoZ966ik9ZcoU7ePjowHdv39//eqrr+rc3FxXv2zRCKjjQHODJQUTA+OAA/Y//E/bzz2HaRUArMQUIU20H4sv9pz1nRQiIyM1oB955BHrA6RAnwb9itOHavTo0dbtRYsW1WsMQlRXVlamv/76az1hwgQrGXh4eFRJFM5HP9D/BH3cniCyQL8E+gqnmUuA9vPz01OnTtUPPfSQHjhwoPW8N9xwg/700091cXGxq1+6aCBNIik0xFHfScHX11cD+vbbb7c+dF3sH6x7nD509957r3U7PT29XmMQ4kKOHTumX375Zd23b1/rPejm5qaB81oPnqAngv4KdKn9fbwP9J9AR1dLEO3bt9e33XabnjZtmg4PDzetD39/fccdd+jly5fL1NbLTF2TgjLXNh/x8fG6vmZTaG0G9cBMMU1OTqaoqIhxwDfAMGCj/doHH3yQ119/naCgIGtHLSEaW1JSEp988gnz588nMzMTpZSjVV7lNkBb4LeYshojMAOIO4GFmFLeB5yeNzw8nKFDh3L27FnWrl1LQUEBQUFB3HzzzUyePJmRI0dae4iI5kkptUNrHX/R61pyUtiyZQtDhgyhVatWlJSUWJnyf4C/Yz5Up+zXTp06lQULFjBy5EhWrFhRL79fiEultWbz5s18/PHHfP755xw/fvy8pOCsEzAFkyQS7OeSMAliEabshkNoaCjx8fGcPXuWrVu3UlhYSGBgIBMmTODmm29m1KhR+Pv7N+CrEw2hrkmhRdc++uijjwBTCtlms1kfqL6YaVKOhNCxY0f27NmDzWaTRWuiSVBKkZCQwNy5c8nJyWH9+vU8/PDDVdbSOMvCTG0dikkQ/w3kAc9iZoCkYRbLjQZOHj3K0qVLWblyJZWVlQwfPpzY2Fi+/vprJk2aRLt27Rg/fjzz5s2z9iERl48W3VLo06cPP//8MwkJCWzatMk6n4JpZk+2358wYQJLlizBZrOxYMECJk+eXNPTCeFyWmt27tzJokWLWLRoET//bFba1NaK6ICpvXQjMApoBRQB3wHLgRVAhtP1ffv2JTAwkPT0dLKysgCzkHP06NGMGTOGhIQEq3aYaFqk+6gOvL29KSsrIzY2loMHD1JcXEwXzLemhzDfnODceAKYvRe6dHFJ3T4hfrGMjAy++eYbvvrqK1avXk2FfbFbTXwwxchuBMYDjp1D9mOSwwrM1qKOujSBgYFERUVx5swZDh48SGVlJQEBAVx33XWMHj2akSNHEhUVJeNvTYQkhYvIzMy0Nszx8/Pj7Nmz2Gw27gDeB/ph+lwBHnroIebOnUtoaCjZ2dnyJhfN0pkzZ1izZg3Lly9nyZIlpKamXvD6GEx30hjMQLUvUIZZMLca+AFTsMyx4jQ0NJTg4GCOHz+Oo/JA586dufbaa62jtu4t0fDqmhQ8GiOYpmj+/PkAhISEcOzYMev81cAJqtY8OnHiBG5ubowYMUISgmi2/Pz8GDt2LGPHjuXVV18lIyOD7777jhUrVrBy5UqrMKTDfvvxKqYV8RvgWvvxDKbMdwmwCVgHrDt6lM1Hj1Js//cdOnTA09OThQsX8t577wEQHR3NVVddxfDhw7nqqqvo1q2bfKaamBbbUkhISGDz5s0MHz6ctWvXWudTge3A7+z3O3fujJ+fH/v372fu3LnMmjXrV/9uIZoarTX79u3jhx9+YNmyZaxZs4aioqJar28NXAWMBIYDsZha+RWY8bj1mGSxGTPIDRAQEEBgYCD5+flWhdiwsDCGDRtGQkICCQkJXHHFFVbdJ1G/pPvoIvz9/SkuLubqq69m06ZNlJaWWuMJs4DX7ddNnjyZL774Aq01u3btIjY29lf/biGaOq01ycnJrFu3jlWrVvHjjz9WaVFXF4CZ6nqV/bgS090EcASTHDZjdpnbBZzGFJgMDAykoqKCwsJCwJSvj4uLIyEhgUGDBhEfH090dLS0JuqBdB9dQHFxMcXFxSilSE5OtvZIuNr++I9O13bs2BGtNb6+vvTr16+RIxXCNZRSxMTEEBMTwz333ANAdnY2W7ZsYf369axevZqkpCRr4Po0ZsbSd/Z/7wn0B4bYjwTMGgkAG5AMbK+sZPvJk+zETIs9jSkGuHfvXrZt22YVBWzTpg0DBw5k0KBBxMXFERsbS1RUlLXwVNSvFpkUFixYAED79u3Jycmxzl8D5HJ+uWyA+Ph467YQLVFYWBgTJ05k4sSJgNmvfM+ePWzfvp21a9eyadMmUlNT0VpTDuywH45Wd3vMtqPx9uNa4Han5z8EJFZUkHj6NInAHsx02IKCAjZu3Mjq1autabUBAQEMGDCA2NhYBgwYQL9+/ejTp48sqqsHLTIpfPbZZwAEBweTn59fpaXwI1U3fcjONtW+HWWNhRCGp6cncXFxxMXFMXPmTABKS0utRLFhwwa2bdtGSkoK5eXl5GLWPix3eo6OmPGIK+w/Yzm3PgigEEjSmqSSEvZgvrDtA3JOn2bTpk1s3LgRm81mXd+tWzf69etHv3796N27N7169aJnz55WCXJxcS1yTCEoKIj8/Hx8fX1p06YNOTk5RGIGmR8E3rBfFxUVhc1mIzU1ldWrV3PNNdf8uuCFaIFsNhuHDx8mMTGRbdu2sXnzZvbv309tG2YFYKaE93U6+gHtnK4pwCSHfZgZUsmcK8dcXm2hnlKKyMhI+vTpQ0xMDD169KBnz5706NGDkJCQFjNeIQPNtbDZbNYOalprPDw8qKio4C7gXczmOo7uo6lTp/LZZ5+hlKKwsJBWrVr9+hcghACgpKSEffv28fPPP7N161Z27NjBoUOHyM3NrXH1dQhmR8Re9sNx23kPxErMVqYHMd1Rh+3HIUylgrPVnrNVq1Z0796dnj17EhUVRXR0NNHR0URFRREaGnpZjVvIQHMtfvjhB8A0ff39/a3Nza8BjlN1PCE0NBStNVFRUZIQhKhnvr6+VvfTtGnTrPPl5eWkpKRw8OBBdu3axY4dO0hOTiYrK4sfior4odrztAa6Az2qHVdiilo6y8H0CKRiZhqmFheTlpjIzsREvqZq0vDw8KBTp05ERUXRo0cPunbtSpcuXejSpQudO3cmJCTkskoaDi0uKbz11luA2Q+3qKgILy8vysrKrPEEZ46+yhEjRjRmiEK0aJ6envTs2ZOePXtyww03VHns7NmzpKamcujQIRITE9m9ezeHDh0iIzubnSdOnNfCaAtEA1H2n5GYrSATMFVjq/8BPI5paWQCmRUVZKWlkZWWxp5Vq1iOKZTpWMHt5uZGSEgInTt3JjIykq5du9K5c2c6depEp06dCA8Pp127ds0ucbS4pODYkxnMNxKtNV0xdV5etJ/39PSkoqKCAwdMxfmxY8c2epxCiPP5+PjQq1cvevXqxY033ljlMa01eXl5pKWlceDAAfbs2UNycjJpaWmsysri87y8KoPS7piKsZ2BLvafjqMnZmFeYA0x5GGSQ7bNRk5ODtk5ORzbsoUUzP4rRzHbSRbYrw8MDCQkJISIiAgrgYSFhREaGmod7du3bzKL9lrcmIKbm5v1bcLPz48zZ84wHbMRem/MwFVoaCgdOnTgxIkTZGdnc/ToUUJCQuojfCGEC5WWlpKTk0NmZib79u1j3759pKamkpWVRVZWFnl5eVWKBgYA4ZjkEYEZvwhzOjraj5rqwpZiWh7HMUniOGbKey6mlM4Jp9t5QKFS+LZqRXBwMBEREYSFhREeHk5ERISVOLp3737JBTlloLkG+/fvp1evXued/wC4Hgh1OvfCCy/wpz/9iTZt2ljjDkKIluH06dMcO3aMrKwsDh48yKFDh0hLS+PIkSOkp6eTm5tLaanpSFKYbqrQakeHt+YYcgAAClhJREFUakcIZq1GbZNjK4F8TILIA07aj3ynn1uALZf4N1sGmmvw8ssvV7nv6+tLSUlJlfEEb29vKioqiI2NRWstZS2EaIECAgIICAggOjqaq6+++oLXlpaWkpeXx/Hjx0lLSyM1NZXMzEwSc3JISUkhMzOTvLw8ysrKAPDDTK9tbz/aAcE1HB0xsyHbAm3sv+uv9f5Kz9egSUEpNQZTZNEdeFtr/WK1x70xX9QHYpLjFK11WkPF8+2331a5X1JSQjdMs/BH+zlPT0+uv/56tm7dCsCYMWMaKhwhxGXA29ubsLAwwsLC6vwl0mazcfr0afLz88nOzraOI8eOsTE9nYMHD5KVlUV+fj4lJSW4YxKDTSmeatBX04BJQSnljlnhPgpTKHGbUmqx1tp51ufdQL7WOlopNRX4G2ZSQINwrE525liO9qP9Z1FREdOnT7daFTfddFNDhSOEaKHc3NwIDAwkMDCQyMhIV4dTRUPOlboSOKS1TtFalwGfAhOqXTMBs6cNmD3ER6oGWl5YW4XHqzGzBfZjKjR26NCB8ePHs3fvXmtqnBBCtBQN2X0Ujpnu65AFDK7tGq11hVKqANOddqK+g5kSGsoqp/tubm506NCBzseOsdQ+cFNeXs60adOorKwkPz+f7t27t5gl8EIIAQ2bFGr6a1p92Lwu16CUmgnMBKwtNC8lmCov1mbjxNGjHAf+1/GLteb3v/+9teo5ISHhkn6XEEI0Vw2ZFLIwY7gOnTBrPmq6Jksp5YFZK3Le/E+t9VvAW2CmpF5KMKtttvNWFjo22tFa4+bmRnx8PH379uXf//43ADfffPOl/CohhGi2GnJMYRvQXSnVVSnlBUwFFle7ZjFwp/32ZGC1bqCFE0opHnrooSrnioqKrIVsNpuN6dOnA1jbc44bN64hQhFCiCarwZKC1roCs7PlCsxC4c+11nuVUs8ppRxTet4BgpVSh4BHgCcaKh6A2bNn1/qYt7c3U6dOZf369SQnJxMVFdVklp0LIURjadB1Clrrb4Fvq5171un2WeB3DRmDs+DgYH7zm9+wfv16lFIopaxaKL/73e8ICAhg6tSpAHz44YeNFZYQQjQZzat8Xz14+umnAVNH3Waz4elpqpY41iYcOXKEwYMHM2TIEFeGKYQQLtGiylwAjBo1iqCgIE6ePElMTAz79+8nMjKSqKgoa/Xyf/7zHxdHKYQQrtHiWgru7u7cd999gJmCGhERwSOPPMLMmTMpKytj1KhR9OvXz8VRCiGEa7S4pABw9913A5CcnMy8efNo164dK1asAOCvf22MklNCCNE0tbjuI4Bu3boxYsQI1q9fz5w5c0hJScHDw4Nrr72W+PiLVpYVQojLVotsKQDMmDGDyspKtm7dysmTJ6moqOCZZ55xdVhCCOFSLTYpTJo0iYCAAFq3bo2/vz/Dhw/nqquucnVYQgjhUi02Kfj5+XHbbbdRWFhIYWGhNVVVCCFashabFACrrMWgQYMYNWqUi6MRQgjXa5EDzQ7x8fHMmTOH8ePHS4lsIYSghScFpdQF6yEJIURL06K7j4QQQlQlSUEIIYRFkoIQQgiLJAUhhBAWSQpCCCEskhSEEEJYJCkIIYSwSFIQQghhUVprV8fwiyilcoH0aqfbASdcEM4v1VziBIm1oTSXWJtLnCCx1lUXrXX7i13U7JJCTZRS27XWTX4jhOYSJ0isDaW5xNpc4gSJtb5J95EQQgiLJAUhhBCWyyUpvOXqAOqoucQJEmtDaS6xNpc4QWKtV5fFmIIQQoj6cbm0FIQQQtSDZp0UlFJjlFLJSqlDSqknXB1PbZRSEUqpH5RS+5RSe5VSf3B1TBeilHJXSu1SSi11dSwXopRqo5RaqJTab/9vm+DqmGqjlPo/9v/3SUqpT5RSPq6OyUEp9f+UUseVUklO54KUUt8rpQ7af7Z1ZYwOtcT6kv098JNSapFSqo0rY3SoKVanx/6olNJKqXauiO1Cmm1SUEq5A68DY4HewC1Kqd6ujapWFcCjWutewBDgwSYcK8AfgH2uDqIOXgWWa61jgAE00ZiVUuHAfwPxWuu+gDsw1bVRVfEeMKbauSeAVVrr7sAq+/2m4D3Oj/V7oK/Wuj9wAHiysYOqxXucHytKqQhgFJDR2AHVRbNNCsCVwCGtdYrWugz4FJjg4phqpLXO0VrvtN8+jfnjFe7aqGqmlOoEjAfednUsF6KUag0MB94B0FqXaa1PuTaqC/IAfJVSHoAfkO3ieCxa67XAyWqnJwDv22+/D9zcqEHVoqZYtdbfaa0r7Hc3A50aPbAa1PLfFeBfwGNAkxzQbc5JIRzIdLqfRRP9Q+tMKRUJXAFscW0ktXoF84a1uTqQi+gG5ALv2ru63lZKtXJ1UDXRWh8B/oH5ZpgDFGitv3NtVBcVorXOAfOlBujg4njqajqwzNVB1EYpdRNwRGu929Wx1KY5JwVVw7kmmXkdlFL+wBfAw1rrQlfHU51S6gbguNZ6h6tjqQMPIA54U2t9BVBM0+niqMLeHz8B6AqEAa2UUtNcG9XlRyn1NKar9iNXx1ITpZQf8DTwrKtjuZDmnBSygAin+51oQk3y6pRSnpiE8JHW+ktXx1OLYcBNSqk0THfctUqpD10bUq2ygCyttaPFtRCTJJqi64BUrXWu1roc+BIY6uKYLuaYUqojgP3ncRfHc0FKqTuBG4DbdNOdZx+F+WKw2/4Z6wTsVEqFujSqappzUtgGdFdKdVVKeWEG7ha7OKYaKaUUpu97n9b6n66OpzZa6ye11p201pGY/56rtdZN8hut1vookKmU6mk/NRL42YUhXUgGMEQp5Wd/L4ykiQ6KO1kM3Gm/fSfwtQtjuSCl1BjgceAmrfUZV8dTG631Hq11B611pP0zlgXE2d/LTUazTQr2gaVZwArMB+xzrfVe10ZVq2HA7Zhv3on2Y5yrg7oMPAR8pJT6CYgF/urieGpkb80sBHYCezCfuyazslUp9QmwCeiplMpSSt0NvAiMUkodxMyUedGVMTrUEutrQADwvf2z9b8uDdKullibPFnRLIQQwtJsWwpCCCHqnyQFIYQQFkkKQgghLJIUhBBCWCQpCCGEsEhSEOIi7NVYH7DfDlNKLXR1TEI0FJmSKsRF2OtVLbVXOBXisubh6gCEaAZeBKKUUonAQaCX1rqvUuouTPVQd6Av8DLghVmoWAqM01qfVEpFYcq8twfOAPdorfc3/ssQ4uKk+0iIi3sCOKy1jgX+p9pjfYFbMaXc/wKcsRfo2wTcYb/mLeAhrfVA4I/AG40StRCXQFoKQvw6P9j3yDitlCoAltjP7wH62yvjDgUWmLJHAHg3fphC1I0kBSF+nVKn2zan+zbM58sNOGVvZQjR5En3kRAXdxpTcO0Xs++bkaqU+h2YirlKqQH1GZwQ9UmSghAXobXOAzbYN2B/6RKe4jbgbqXUbmAvTXTbWCFApqQKIYRwIi0FIYQQFkkKQgghLJIUhBBCWCQpCCGEsEhSEEIIYZGkIIQQwiJJQQghhEWSghBCCMv/B3otj6nKSHq9AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(output.loc[output.index < 15,output.columns.str.startswith('ensemble')],'k')\n", + "plt.plot(output.loc[output.index < 15,'x'],'r')\n", + "plt.xlabel('time')\n", + "plt.ylabel('temperature')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Remember to remove the models from memory" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "x.finalize()\n", + "for ensembleMember in range(nEnsemble):\n", + " ensemble[ensembleMember].finalize()" + ] + } + ], + "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.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/run-model-from-bmi.ipynb b/examples/run-model-from-bmi.ipynb index ad1c650..c4ffdd0 100644 --- a/examples/run-model-from-bmi.ipynb +++ b/examples/run-model-from-bmi.ipynb @@ -23,7 +23,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -42,7 +42,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -58,9 +58,17 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The 2D Heat Equation\n" + ] + } + ], "source": [ "print(x.get_component_name())" ] @@ -74,16 +82,34 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "# Heat model configuration\n", + "shape:\n", + " - 6\n", + " - 8\n", + "spacing:\n", + " - 1.0\n", + " - 1.0\n", + "origin:\n", + " - 0.0\n", + " - 0.0\n", + "alpha: 1.0\n" + ] + } + ], "source": [ "cat heat.yaml" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -99,9 +125,21 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Start time: 0.0\n", + "End time: 1.7976931348623157e+308\n", + "Current time: 0.0\n", + "Time step: 0.25\n", + "Time units: s\n" + ] + } + ], "source": [ "print(\"Start time:\", x.get_start_time())\n", "print(\"End time:\", x.get_end_time())\n", @@ -326,7 +364,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -340,9 +378,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.5" + "version": "3.7.1" } }, "nbformat": 4, - "nbformat_minor": 1 + "nbformat_minor": 2 } diff --git a/heat/.ipynb_checkpoints/bmi_heat-checkpoint.py b/heat/.ipynb_checkpoints/bmi_heat-checkpoint.py new file mode 100644 index 0000000..63207ec --- /dev/null +++ b/heat/.ipynb_checkpoints/bmi_heat-checkpoint.py @@ -0,0 +1,399 @@ +#! /usr/bin/env python +"""Basic Model Interface implementation for the 2D heat model.""" + +import numpy as np +import json +from bmipy import Bmi + +from .heat import Heat + + +class BmiHeat(Bmi): + + """Solve the heat equation for a 2D plate.""" + + _name = "The 2D Heat Equation" + _input_var_names = ("plate_surface__temperature",) + _output_var_names = ("plate_surface__temperature",) + _state_var_names = ("plate_surface__temperature",) + + def __init__(self): + """Create a BmiHeat model that is ready for initialization.""" + self._model = None + self._values = {} + self._var_units = {} + self._var_loc = {} + self._grids = {} + self._grid_type = {} + + self._start_time = 0.0 + self._end_time = np.finfo("d").max + self._time_units = "s" + + def initialize(self, filename=None): + """Initialize the Heat model. + + Parameters + ---------- + filename : str, optional + Path to name of input file. + """ + if filename is None: + self._model = Heat() + elif isinstance(filename, str): + with open(filename, "r") as file_obj: + self._model = Heat.from_file_like(file_obj.read()) + else: + self._model = Heat.from_file_like(filename) + + self._values = {"plate_surface__temperature": self._model.temperature} + self._var_units = {"plate_surface__temperature": "K"} + self._var_loc = {"plate_surface__temperature": "node"} + self._grids = {0: ["plate_surface__temperature"]} + self._grid_type = {0: "uniform_rectilinear"} + + def update(self): + """Advance model by one time step.""" + self._model.advance_in_time() + + def update_frac(self, time_frac): + """Update model by a fraction of a time step. + + Parameters + ---------- + time_frac : float + Fraction fo a time step. + """ + time_step = self.get_time_step() + self._model.time_step = time_frac * time_step + self.update() + self._model.time_step = time_step + + def update_until(self, then): + """Update model until a particular time. + + Parameters + ---------- + then : float + Time to run model until. + """ + n_steps = (then - self.get_current_time()) / self.get_time_step() + + for _ in range(int(n_steps)): + self.update() + self.update_frac(n_steps - int(n_steps)) + + def finalize(self): + """Finalize model.""" + self._model = None + + def get_var_type(self, var_name): + """Data type of variable. + + Parameters + ---------- + var_name : str + Name of variable as CSDMS Standard Name. + + Returns + ------- + str + Data type. + """ + return str(self.get_value_ptr(var_name).dtype) + + def get_var_units(self, var_name): + """Get units of variable. + + Parameters + ---------- + var_name : str + Name of variable as CSDMS Standard Name. + + Returns + ------- + str + Variable units. + """ + return self._var_units[var_name] + + def get_var_nbytes(self, var_name): + """Get units of variable. + + Parameters + ---------- + var_name : str + Name of variable as CSDMS Standard Name. + + Returns + ------- + int + Size of data array in bytes. + """ + return self.get_value_ptr(var_name).nbytes + + def get_var_itemsize(self, name): + return np.dtype(self.get_var_type(name)).itemsize + + def get_var_location(self, name): + return self._var_loc[name] + + def get_var_grid(self, var_name): + """Grid id for a variable. + + Parameters + ---------- + var_name : str + Name of variable as CSDMS Standard Name. + + Returns + ------- + int + Grid id. + """ + for grid_id, var_name_list in self._grids.items(): + if var_name in var_name_list: + return grid_id + + def get_grid_rank(self, grid_id): + """Rank of grid. + + Parameters + ---------- + grid_id : int + Identifier of a grid. + + Returns + ------- + int + Rank of grid. + """ + return len(self._model.shape) + + def get_grid_size(self, grid_id): + """Size of grid. + + Parameters + ---------- + grid_id : int + Identifier of a grid. + + Returns + ------- + int + Size of grid. + """ + return int(np.prod(self._model.shape)) + + def get_value_ptr(self, var_name): + """Reference to values. + + Parameters + ---------- + var_name : str + Name of variable as CSDMS Standard Name. + + Returns + ------- + array_like + Value array. + """ + return self._values[var_name] + + def get_value(self, var_name, dest): + """Copy of values. + + Parameters + ---------- + var_name : str + Name of variable as CSDMS Standard Name. + dest : ndarray + A numpy array into which to place the values. + + Returns + ------- + array_like + Copy of values. + """ + dest[:] = self.get_value_ptr(var_name).flatten() + return dest + + def get_value_at_indices(self, var_name, dest, indices): + """Get values at particular indices. + + Parameters + ---------- + var_name : str + Name of variable as CSDMS Standard Name. + dest : ndarray + A numpy array into which to place the values. + indices : array_like + Array of indices. + + Returns + ------- + array_like + Values at indices. + """ + dest[:] = self.get_value_ptr(var_name).take(indices) + return dest + + def get_state(self): + outDict = {'time' : self.get_current_time()} + for var in self._state_var_names: + varValue = [] + varType = [] + varSize = [] + self.get_value(var, varValue) + outDict[var] = {'value': varValue, + 'type': self.get_var_type(var), + 'itemsize': self.get_var_itemsize(var), + 'nbytes': self.get_var_nbytes(var)} + + + return json.dumps(outDict, indent = 4) + + def get_state_ptr(self, state_ptr): + return "not implemented" + + + def set_value(self, var_name, src): + """Set model values. + + Parameters + ---------- + var_name : str + Name of variable as CSDMS Standard Name. + src : array_like + Array of new values. + """ + val = self.get_value_ptr(var_name) + val[:] = src.reshape(val.shape) + + def set_value_at_indices(self, name, inds, src): + """Set model values at particular indices. + + Parameters + ---------- + var_name : str + Name of variable as CSDMS Standard Name. + src : array_like + Array of new values. + indices : array_like + Array of indices. + """ + val = self.get_value_ptr(name) + val.flat[inds] = src + + def set_state(self, state): + inDict = json.loads(state) + for key in inDict: + if key == 'time': + self._model._time = inDict[key] + else: + # we don't have to check for type, since values are + # always numpy arrays in BMI-Python + self.set_value(key,np.array(inDict[key]['value'])) + + + def set_state_ptr(self, state_ptr): + return "not implemented" + + def get_component_name(self): + """Name of the component.""" + return self._name + + def get_input_item_count(self): + """Get names of input variables.""" + return len(self._input_var_names) + + def get_output_item_count(self): + """Get names of output variables.""" + return len(self._output_var_names) + + def get_input_var_names(self): + """Get names of input variables.""" + return self._input_var_names + + def get_output_var_names(self): + """Get names of output variables.""" + return self._output_var_names + + def get_grid_shape(self, grid_id, shape): + """Number of rows and columns of uniform rectilinear grid.""" + var_name = self._grids[grid_id][0] + shape[:] = self.get_value_ptr(var_name).shape + return shape + + def get_grid_spacing(self, grid_id, spacing): + """Spacing of rows and columns of uniform rectilinear grid.""" + spacing[:] = self._model.spacing + return spacing + + def get_grid_origin(self, grid_id, origin): + """Origin of uniform rectilinear grid.""" + origin[:] = self._model.origin + return origin + + def get_grid_type(self, grid_id): + """Type of grid.""" + return self._grid_type[grid_id] + + def get_start_time(self): + """Start time of model.""" + return self._start_time + + def get_end_time(self): + """End time of model.""" + return self._end_time + + def get_current_time(self): + return self._model.time + + def get_time_step(self): + return self._model.time_step + + def get_time_units(self): + return self._time_units + + def get_grid_edge_count(self, grid): + raise NotImplementedError("get_grid_edge_count") + + def get_grid_edge_nodes(self, grid, edge_nodes): + raise NotImplementedError("get_grid_edge_nodes") + + def get_grid_face_count(self, grid): + raise NotImplementedError("get_grid_face_count") + + def get_grid_face_nodes(self, grid, face_nodes): + raise NotImplementedError("get_grid_face_nodes") + + def get_grid_node_count(self, grid): + """Number of grid nodes. + + Parameters + ---------- + grid : int + Identifier of a grid. + + Returns + ------- + int + Size of grid. + """ + return self.get_grid_size(grid) + + def get_grid_nodes_per_face(self, grid, nodes_per_face): + raise NotImplementedError("get_grid_nodes_per_face") + + def get_grid_face_edges(self, grid, face_edges): + raise NotImplementedError("get_grid_face_edges") + + def get_grid_x(self, grid, x): + raise NotImplementedError("get_grid_x") + + def get_grid_y(self, grid, y): + raise NotImplementedError("get_grid_y") + + def get_grid_z(self, grid, z): + raise NotImplementedError("get_grid_z") diff --git a/heat/.ipynb_checkpoints/heat-checkpoint.py b/heat/.ipynb_checkpoints/heat-checkpoint.py new file mode 100644 index 0000000..e940511 --- /dev/null +++ b/heat/.ipynb_checkpoints/heat-checkpoint.py @@ -0,0 +1,184 @@ +"""The 2D heat model.""" + +import numpy as np +import yaml +from scipy import ndimage + + +def solve_2d(temp, spacing, out=None, alpha=1.0, time_step=1.0): + """Solve the 2D Heat Equation on a uniform mesh. + + Parameters + ---------- + temp : ndarray + Temperature. + spacing : array_like + Grid spacing in the row and column directions. + out : ndarray (optional) + Output array. + alpha : float (optional) + Thermal diffusivity. + time_step : float (optional) + Time step. + + Returns + ------- + result : ndarray + The temperatures after time *time_step*. + + Examples + -------- + >>> from heat import solve_2d + >>> z0 = np.zeros((3, 3)) + >>> z0[1:-1, 1:-1] = 1. + >>> solve_2d(z0, (1., 1.), alpha=.25) + array([[0. , 0. , 0. ], + [0. , 0.5, 0. ], + [0. , 0. , 0. ]]) + """ + dy2, dx2 = spacing[0] ** 2, spacing[1] ** 2 + stencil = ( + np.array([[0.0, dy2, 0.0], [dx2, -2.0 * (dx2 + dy2), dx2], [0.0, dy2, 0.0]]) + * alpha + * time_step + / (2.0 * (dx2 * dy2)) + ) + + if out is None: + out = np.empty_like(temp) + + ndimage.convolve(temp, stencil, output=out) + out[(0, -1), :] = 0.0 + out[:, (0, -1)] = 0.0 + return np.add(temp, out, out=out) + + +class Heat(object): + + """Solve the Heat equation on a grid. + + Examples + -------- + >>> heat = Heat() + >>> heat.time + 0.0 + >>> heat.time_step + 0.25 + >>> heat.advance_in_time() + >>> heat.time + 0.25 + + >>> heat = Heat(shape=(5, 5)) + >>> heat.temperature = np.zeros_like(heat.temperature) + >>> heat.temperature[2, 2] = 1. + >>> heat.advance_in_time() + + >>> heat = Heat(alpha=.5) + >>> heat.time_step + 0.5 + >>> heat = Heat(alpha=.5, spacing=(2., 3.)) + >>> heat.time_step + 2.0 + """ + + def __init__( + self, shape=(10, 20), spacing=(1.0, 1.0), origin=(0.0, 0.0), alpha=1.0 + ): + """Create a new heat model. + + Parameters + --------- + shape : array_like, optional + The shape of the solution grid as (*rows*, *columns*). + spacing : array_like, optional + Spacing of grid rows and columns. + origin : array_like, optional + Coordinates of lower left corner of grid. + alpha : float + Alpha parameter in the heat equation. + """ + self._shape = shape + self._spacing = spacing + self._origin = origin + self._time = 0.0 + self._alpha = alpha + self._time_step = min(spacing) ** 2 / (4.0 * self._alpha) + + self._temperature = np.random.random(self._shape) + self._next_temperature = np.empty_like(self._temperature) + + @property + def time(self): + """Current model time.""" + return self._time + + @property + def temperature(self): + """Temperature of the plate.""" + return self._temperature + + @temperature.setter + def temperature(self, new_temp): + """Set the temperature of the plate. + + Parameters + ---------- + new_temp : array_like + The new temperatures. + """ + self._temperature[:] = new_temp + + @property + def time_step(self): + """Model time step.""" + return self._time_step + + @time_step.setter + def time_step(self, time_step): + """Set model time step.""" + self._time_step = time_step + + @property + def shape(self): + """Shape of the model grid.""" + return self._shape + + @property + def spacing(self): + """Spacing between nodes of the model grid.""" + return self._spacing + + @property + def origin(self): + """Origin coordinates of the model grid.""" + return self._origin + + @classmethod + def from_file_like(cls, file_like): + """Create a Heat object from a file-like object. + + Parameters + ---------- + file_like : file_like + Input parameter file. + + Returns + ------- + Heat + A new instance of a Heat object. + """ + config = yaml.safe_load(file_like) + return cls(**config) + + def advance_in_time(self): + """Calculate new temperatures for the next time step.""" + solve_2d( + self._temperature, + self._spacing, + out=self._next_temperature, + alpha=self._alpha, + time_step=self._time_step, + ) + np.copyto(self._temperature, self._next_temperature) + + self._time += self._time_step diff --git a/heat/bmi_heat.py b/heat/bmi_heat.py index 310855a..63207ec 100644 --- a/heat/bmi_heat.py +++ b/heat/bmi_heat.py @@ -2,6 +2,7 @@ """Basic Model Interface implementation for the 2D heat model.""" import numpy as np +import json from bmipy import Bmi from .heat import Heat @@ -14,6 +15,7 @@ class BmiHeat(Bmi): _name = "The 2D Heat Equation" _input_var_names = ("plate_surface__temperature",) _output_var_names = ("plate_surface__temperature",) + _state_var_names = ("plate_surface__temperature",) def __init__(self): """Create a BmiHeat model that is ready for initialization.""" @@ -236,6 +238,25 @@ def get_value_at_indices(self, var_name, dest, indices): dest[:] = self.get_value_ptr(var_name).take(indices) return dest + def get_state(self): + outDict = {'time' : self.get_current_time()} + for var in self._state_var_names: + varValue = [] + varType = [] + varSize = [] + self.get_value(var, varValue) + outDict[var] = {'value': varValue, + 'type': self.get_var_type(var), + 'itemsize': self.get_var_itemsize(var), + 'nbytes': self.get_var_nbytes(var)} + + + return json.dumps(outDict, indent = 4) + + def get_state_ptr(self, state_ptr): + return "not implemented" + + def set_value(self, var_name, src): """Set model values. @@ -264,6 +285,20 @@ def set_value_at_indices(self, name, inds, src): val = self.get_value_ptr(name) val.flat[inds] = src + def set_state(self, state): + inDict = json.loads(state) + for key in inDict: + if key == 'time': + self._model._time = inDict[key] + else: + # we don't have to check for type, since values are + # always numpy arrays in BMI-Python + self.set_value(key,np.array(inDict[key]['value'])) + + + def set_state_ptr(self, state_ptr): + return "not implemented" + def get_component_name(self): """Name of the component.""" return self._name