diff --git a/docs/source/api.rst b/docs/source/api.rst index c37f9e6d..c59c6b08 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -21,6 +21,7 @@ Modules, methods, classes and attributes are explained here. fatigue humidity letid + montecarlo spectral standards temperature diff --git a/docs/source/whatsnew/releases/v0.3.0.rst b/docs/source/whatsnew/releases/v0.3.0.rst new file mode 100644 index 00000000..d2ee6892 --- /dev/null +++ b/docs/source/whatsnew/releases/v0.3.0.rst @@ -0,0 +1,13 @@ +v0.3.0 (2024-02-12) +======================= + +Enhancements +------------ +* Initial integration of the Monte-Carlo analysis package + +Contributors +~~~~~~~~~~~~ +* Tobin Ford (:ghuser:`tobin-ford`) +* Silvana Ovaitt (:ghuser:`shirubana`) +* Martin Springer (:ghuser:`martin-springer`) +* Mike Kempe (:ghuser:`MDKempe`) diff --git a/pvdeg/__init__.py b/pvdeg/__init__.py index 199087b1..301c9abb 100644 --- a/pvdeg/__init__.py +++ b/pvdeg/__init__.py @@ -11,6 +11,7 @@ from . import geospatial from . import humidity from . import letid +from . import montecarlo from .scenario import Scenario from . import spectral from . import standards diff --git a/pvdeg/degradation.py b/pvdeg/degradation.py index e776e90f..9ef54acf 100644 --- a/pvdeg/degradation.py +++ b/pvdeg/degradation.py @@ -3,7 +3,7 @@ import numpy as np import pandas as pd -from numba import jit +from numba import jit, njit from rex import NSRDBX from rex import Outputs from pathlib import Path @@ -722,3 +722,55 @@ def degradation( degradation = C * data["dD"].sum(axis=0) return degradation + +# change it to take pd.DataFrame? instead of np.ndarray +@njit +def vecArrhenius( + poa_global : np.ndarray, + module_temp : np.ndarray, + ea : float, + x : float, + lnr0 : float + ) -> float: + + """ + Calculates degradation using :math:`R_D = R_0 * I^X * e^{\\frac{-Ea}{kT}}` + + Parameters + ---------- + poa_global : numpy.ndarray + Plane of array irradiance [W/m^2] + + module_temp : numpy.ndarray + Cell temperature [C]. + + ea : float + Activation energy [kJ/mol] + + x : float + Irradiance relation [unitless] + + lnR0 : float + prefactor [ln(%/h)] + + Returns + ---------- + degredation : float + Degradation Rate [%/h] + + """ + + mask = poa_global >= 25 + poa_global = poa_global[mask] + module_temp = module_temp[mask] + + ea_scaled = ea / 8.31446261815324E-03 + R0 = np.exp(lnr0) + poa_global_scaled = poa_global / 1000 + + degredation = 0 + # refactor to list comprehension approach + for entry in range(len(poa_global_scaled)): + degredation += R0 * np.exp(-ea_scaled / (273.15 + module_temp[entry])) * np.power(poa_global_scaled[entry], x) + + return (degredation / len(poa_global)) diff --git a/pvdeg/montecarlo.py b/pvdeg/montecarlo.py new file mode 100644 index 00000000..74fd3732 --- /dev/null +++ b/pvdeg/montecarlo.py @@ -0,0 +1,325 @@ +""" +Collection of functions for monte carlo simulations. +""" +import numpy as np +import pandas as pd +from numba import njit +from scipy.linalg import cholesky +from scipy import stats +from typing import Callable +import inspect + +class Corr: + """ + corrlation class : + stores modeling constants and corresponding correlation coefficient to access at runtime + """ + + # modeling constants : str + mc_1 = '' + mc_2 = '' + # corresonding corelation coefficient : float + correlation = 0 + + def __init__(self, mc_1_string, mc_2_string, corr): + """parameterized constructor""" + self.mc_1 = mc_1_string + self.mc_2 = mc_2_string + self.correlation = corr + + def getModelingConstants(self)->list[str, str]: + """ + Helper method. Returns modeling constants in string form. + + Parameters + ---------- + self : Corr + Reference to self + + Returns + ---------- + modeling_constants : list[str, str] + Both modeling constants in string from from their corresponding correlation coefficient object + """ + + modeling_constants = [self.mc_1, self.mc_2] + return modeling_constants + +def _symettric_correlation_matrix(corr: list[Corr])->pd.DataFrame: + """ + Helper function. Generate a symmetric correlation coefficient matrix. + + Parameters + ---------- + corr : list[Corr] + All correlations between appropriate modeling constants + + Returns + ---------- + identity_df : pd.DataFrame + Matrix style DataFrame containing relationships between all input modeling constants + Index and Column names represent modeling constants for comprehensibility + """ + + if not corr: + return None + + # unpack individual modeling constants from correlations + modeling_constants = [mc for i in corr for mc in i.getModelingConstants()] + + uniques = np.unique(modeling_constants) + + # setting up identity matrix, labels for columns and rows + identity_matrix = np.eye(len(uniques)) + identity_df = pd.DataFrame(identity_matrix, columns = uniques, index=uniques) + + # walks matrix to fill in correlation coefficients + # make this a modular standalone function if bigger function preformance is not improved with @njit + for i in range(len(uniques)): + for j in range(i): # only iterate over lower triangle + x, y = identity_df.index[i], identity_df.columns[j] + + # find the correlation coefficient + found = False + for relation in corr: + if set([x, y]) == set(relation.getModelingConstants()): + # fill in correlation coefficient + identity_df.iat[i, j] = relation.correlation + found = True + break + + # if no matches in all correlation coefficients, they will be uncorrelated (= 0) + if not found: + identity_df.iat[i, j] = 0 + + # mirror the matrix + # this may be computationally expensive for large matricies + # could be better to fill the original matrix in all in one go rather than doing lower triangular and mirroring it across I + identity_df = identity_df + identity_df.T - np.diag(identity_df.to_numpy().diagonal()) + + # identity_df should be renamed more appropriately + return identity_df + +def _createStats( + stats : dict[str, dict[str, float]], + corr : list[Corr] + ) -> pd.DataFrame: + + """ + helper function. Unpacks mean and standard deviation for modeling constants into a DataFrame + + Parameters + ---------- + stats : dict[str, dict[str, float]] + contains mean and standard deviation for each modeling constant + example of one mc: {'Ea' : {'mean' : 62.08, 'stdev' : 7.3858 }} + + Returns + ---------- + stats_df : pd.DataFrame + contains unpacked means and standard deviations from dictionary + """ + + # empty correlation list case + if not corr: + stats_df = pd.DataFrame(stats) + return stats_df + + + # incomplete dataset + for mc in stats: + if 'mean' not in stats[mc] or 'stdev' not in stats[mc]: + raise ValueError(f"Missing 'mean' or 'stdev' for modeling constant") + + # unpack data + modeling_constants = list(stats.keys()) + mc_mean = [stats[mc]['mean'] for mc in modeling_constants] + mc_stdev = [stats[mc]['stdev'] for mc in modeling_constants] + + stats_df = pd.DataFrame({'mean' : mc_mean, 'stdev' : mc_stdev}, index=modeling_constants).T + + + # flatten and reorder + modeling_constants = [mc for i in corr for mc in i.getModelingConstants()] + uniques = np.unique(modeling_constants) + + # what happens if columns do not match? + if len(uniques) != len(corr): + raise ValueError(f"correlation data is insufficient") + + # should match columns from correlation matrix + stats_df = stats_df[uniques] + + return stats_df + +def _correlateData( + samples_to_correlate : pd.DataFrame, + stats_for_correlation : pd.DataFrame + ) -> pd.DataFrame: + + """ + helper function. Uses meaningless correlated samples and makes meaningful by + multiplying random samples by their parent modeling constant's standard deviation + and adding the mean + + Parameters + ---------- + samples_to_correlate : pd.DataFrame + contains n samples generated with N(0, 1) for each modeling constant + column names must be consistent with all modeling constant inputs + stats_for_correlation : pd.DataFrame + contains mean and stdev each modeling constant, + column names must be consistent with all modeling constant inputs + + Returns + ---------- + correlated_samples : pd.DataFrame + correlated samples in a tall dataframe. column names match modeling constant inputs, + integer indexes. See generateCorrelatedSamples() references section for process info + """ + + # accounts for out of order column names, AS LONG AS ALL MATCH + # UNKNOWN CASE: what will happen if there is an extra NON matching column in stats + columns = list(samples_to_correlate.columns.values) + ordered_stats = stats_for_correlation[columns] + + means = ordered_stats.loc['mean'] + stdevs = ordered_stats.loc['stdev'] + + correlated_samples = samples_to_correlate.multiply(stdevs).add(means) + + return correlated_samples + +def generateCorrelatedSamples( + corr : list[Corr], + stats : dict[str, dict[str, float]], + n : int, + seed = None + ) -> pd.DataFrame: + + # columns are now named, may run into issues if more mean and stdev entries than correlation coefficients + # havent tested yet but this could cause major issues (see lines 163 and 164 for info) + + """ + Generates a tall correlated samples numpy array based on correlation coefficients and mean and stdev + for modeling constants. Values are correlated from cholesky decomposition of correlation coefficients, + and n random samples for each modeling constant generated from a standard distribution with mean = 0 + and standard deviation = 1. + + Parameters + ---------- + corr : List[Corr] + list containing correlations between variable + stats : dict[str, dict[str, float]] + dictionary storing variable mean and standard deviation. Syntax : ` : {'mean' : , 'stdev' : }` + n : int + number of samples to create + seed : Any, optional + reseed the numpy BitGenerator, numpy legacy function (use cautiously) + + Returns + ---------- + correlated_samples : pd.Dataframe + tall dataframe of dimensions (n by # of modeling constants). + Columns named as modeling constants from Corr object inputs + + References + ---------- + Burgess, Nicholas, Correlated Monte Carlo Simulation using Cholesky Decomposition (March 25, 2022). + Available at SSRN: https://ssrn.com/abstract=4066115 + """ + + if seed: + np.random.seed(seed=seed) + + # base case + if corr: + coeff_matrix = _symettric_correlation_matrix(corr) # moved inside + + decomp = cholesky(coeff_matrix.to_numpy(), lower = True) + + # list of correlations + # using to check if all r = 0 + values = [] + for i in corr: + values.append(i.correlation) + + # check if all zero + all_zeros = all(value == 0 for value in values) + + samples = np.random.normal(loc=0, scale=1, size=(len(stats), n)) + + stats_df = _createStats(stats, corr) + + # no correlation data given, only stats + # OR, all correlations are 0 + if (not corr) or (all_zeros): + nocorr_df = pd.DataFrame(samples.T, columns=stats_df.columns.tolist()) + + meaningful_nocorr_df = _correlateData(nocorr_df, stats_df) + + return meaningful_nocorr_df + + if corr: + precorrelated_samples = np.matmul(decomp, samples) + + precorrelated_df = pd.DataFrame(precorrelated_samples.T, columns=coeff_matrix.columns.to_list()) + + correlated_df = _correlateData(precorrelated_df, stats_df) + + return correlated_df + +# monte carlo function +# model after - https://github.com/NREL/PVDegradationTools/blob/main/pvdeg_tutorials/tutorials/LETID%20-%20Outdoor%20Geospatial%20Demo.ipynb + +def simulate( + func : Callable, + correlated_samples : pd.DataFrame, + **function_kwargs + ) -> pd.Series: + + """ + Applies a target function to data to preform a monte carlo simulation. If you get a key error and the target function has default parameters, + try adding them to your ``func_kwargs`` dictionary instead of using the default value from the target function. + + Parameters + ---------- + func : function + Function to apply for monte carlo simulation + correlated_samples : pd.DataFrame + Dataframe of correlated samples with named columns for each appropriate modeling constant, can be generated using generateCorrelatedSamples() + function_kwargs : dict + Keyword arguments to pass to func, only include arguments not named in your correlated_samples columns + + Returns + ------- + res : pandas.Series + Series with monte carlo results from target function + """ + + ### NOTES ### + # func modeling constant parameters must be lowercase in function definition + # dynamically construct argument list for func + # call func with .apply(lambda) + + args = {k.lower(): v for k, v in function_kwargs.items()} # make lowercase + + func_signature = inspect.signature(func) + + func_args = set(func_signature.parameters.keys()) + + def prepare_args(row): + return {arg: row[arg] if arg in row else function_kwargs.get(arg) for arg in func_args} + + args = prepare_args(correlated_samples.iloc[0]) + + def apply_func(row): + row_args = {**args, **{k.lower(): v for k, v in row.items()}} + + return func(**row_args) + + # this line is often flagged when target function is not given required arguments + # problems also arise when target function parameter names are not lowercase + result = correlated_samples.apply(apply_func, axis=1) + + return result \ No newline at end of file diff --git a/pvdeg/standards.py b/pvdeg/standards.py index 3c8f41d6..a41f098e 100644 --- a/pvdeg/standards.py +++ b/pvdeg/standards.py @@ -242,7 +242,7 @@ def standoff( conf_inf : str, optional Model for the lowest temperature module on the exponential decay curve. Default: 'open_rack_glass_polymer' - x0 : float, optional + x_0 : float, optional Thermal decay constant (cm), [Kempe, PVSC Proceedings 2023] wind_factor : float, optional Wind speed correction exponent to account for different wind speed measurement heights diff --git a/pvdeg/temperature.py b/pvdeg/temperature.py index 5b381a13..690090a4 100644 --- a/pvdeg/temperature.py +++ b/pvdeg/temperature.py @@ -177,7 +177,7 @@ def cell( # this one does a linear conversion from the other models, faiman, pvsyst, noct_sam, sapm_module and generic_linear. # An appropriate facter will need to be figured out. else: - wind_speed_factor = 1 # this is just hear for completeness. + wind_speed_factor = 1 # this is just here for completeness. parameters = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS[temp_model][conf] if poa is None: diff --git a/pvdeg_tutorials/tutorials/Monte Carlo - Arrhenius.ipynb b/pvdeg_tutorials/tutorials/Monte Carlo - Arrhenius.ipynb new file mode 100644 index 00000000..0f99797f --- /dev/null +++ b/pvdeg_tutorials/tutorials/Monte Carlo - Arrhenius.ipynb @@ -0,0 +1,452 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Monte Carlo - Arrhenius Degradation \n", + "\n", + "Author: Tobin Ford | tobin.ford@nrel.gov\n", + "\n", + "2023\n", + "***\n", + "short intro here\n", + "process and relevance (why shold the reader care, justification)\n", + "\n", + "A monte carlo simulation can be used to predict results of an event with a certain amount of uncertainty. This will be introduced to our use case via mean and standard deviation for each modeling constant. Correlated multivariate monte carlo simulations expand on this by linking the behavior of multiple input variables together with correlation data, in our case we will use correlation coefficients but \n", + "\n", + "**Objectives**\n", + "1. Define necessary monte carlo simulation parameters : correlation coefficients, mean and standard standard deviation, number of trials, function to apply, requried function input\n", + "2. Define process for creating and utilizing modeling constant correlation data\n", + "3. Preform simple monte carlo simulation using arrhenius equation to calculate degredation and plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# if running on google colab, uncomment the next line and execute this cell to install the dependencies and prevent \"ModuleNotFoundError\" in later cells:\n", + "# !pip install pvdeg==0.2.0" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import pvlib\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pvdeg\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Working on a Windows 10\n", + "Python version 3.11.4 | packaged by Anaconda, Inc. | (main, Jul 5 2023, 13:38:37) [MSC v.1916 64 bit (AMD64)]\n", + "Pandas version 2.1.0\n", + "Pvlib version 0.9.5\n", + "pvdeg version 0.2.0+40.g968e483.dirty\n" + ] + } + ], + "source": [ + "# This information helps with debugging and getting support :)\n", + "import sys, platform\n", + "print(\"Working on a \", platform.system(), platform.release())\n", + "print(\"Python version \", sys.version)\n", + "print(\"Pandas version \", pd.__version__)\n", + "print(\"Pvlib version \", pvlib.__version__)\n", + "print(\"pvdeg version \", pvdeg.__version__)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Correlated Monte Carlo Simulation (parameters)\n", + "\n", + "For this simulation we will be using an arrhenius equation to calculate degredation rate given by $R_D = R_0 * I ^ X * e ^ {\\frac{-Ea}{kT}}$, where R0 is prefactor degredation, I is irradiance, X is the irridiance relation, Ea is activation energy and T is degrees K\n", + "\n", + "We will use R0, X and Ea to preform a 3 variable monte carlo simulation to calculate degredation.\n", + "\n", + "### Required inputs\n", + "To run a monte carlo simulation with pvdeg.montecarlo the following inputs will be required\n", + "- function (currently only works with pvdeg.montecarlo.vecArrhenius() but will eventually work with most pvdeg calculation functions)\n", + "- function arguments (ex: metadata, weather, cell temperature, solar position, etc.)\n", + "- mean and standard deviation (Required for all correlation constants)\n", + "- correlation constants (if not entered, default = 0)\n", + "- number of trials to run" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Defining Correlation Coefficients\n", + "pvdeg.montecarlo stores correlation coefficients in a ``Corr`` object. To represent a given correlation coefficient follow the given syntax below, replacing the values in the brackets with your correlation coefficients\n", + "\n", + " {my_correlation_object} = Corr('{variable1}', '{variable2}', {correlation coefficient})\n", + "\n", + "note: ordering of `variable1` and `variable2` does not matter\n", + "\n", + "After defining the all known correlations add them to a list which we will feed into our simulation later" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "corr_Ea_X = pvdeg.montecarlo.Corr('Ea', 'X', 0.0269)\n", + "corr_Ea_LnR0 = pvdeg.montecarlo.Corr('Ea', 'LnR0', -0.9995)\n", + "corr_X_LnR0 = pvdeg.montecarlo.Corr('X', 'LnR0', -0.0400)\n", + "\n", + "corr_coeff = [corr_Ea_X, corr_Ea_LnR0, corr_X_LnR0]" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "pvdeg.montecarlo.Corr" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "type(corr_Ea_X)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Defining Mean and Standard Deviation\n", + "We will store the mean and correlation for each variable, expressed when we defined the correlation cefficients. If a variable is left out at this stage, the monte carlo simulation will throw errors." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "stats_dict = {\n", + " 'Ea' : {'mean' : 62.08, 'stdev' : 7.3858 },\n", + " 'LnR0' : {'mean' : 13.7223084 , 'stdev' : 2.47334772},\n", + " 'X' : {'mean' : 0.0341 , 'stdev' : 0.0992757 }\n", + "}\n", + "\n", + "# and number of monte carlo trials to run\n", + "n = 20000" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generating Monte Carlo Input Data\n", + "Next we will use the information collected above to generate correlated data from our modeling constant correlations, means and standard deviations." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Ea LnR0 X\n", + "0 43.310964 20.082424 -0.076470\n", + "1 63.858966 13.259692 0.024310\n", + "2 59.968897 14.428516 0.042583\n", + "3 73.746862 9.652391 0.295242\n", + "4 67.954485 11.666279 0.029110\n", + "... ... ... ...\n", + "19995 57.022304 15.367897 0.105352\n", + "19996 69.850571 11.244541 0.012642\n", + "19997 55.831637 16.003457 0.004917\n", + "19998 51.458692 17.205392 0.125170\n", + "19999 52.547601 16.893998 0.240680\n", + "\n", + "[20000 rows x 3 columns]\n" + ] + } + ], + "source": [ + "mc_inputs = pvdeg.montecarlo.generateCorrelatedSamples(corr=corr_coeff, stats=stats_dict, n=n)\n", + "print(mc_inputs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Sanity Check\n", + "We can observe the mean and standard deviation of our newly correlated samples before using them for calculations to ensure that we have not incorrectly altered the data. The mean and standard deviation should be the similar (within a range) to your original input (the error comes from the standard distribution of generated random numbers)\n", + "\n", + "This also applies to the correlation coefficients originally inputted, they should be witin the same range as those orginally supplied." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ea : mean 62.00571434068475, stdev 7.2945584039777325\n", + "LnR0 : mean 13.747061024151163, stdev 2.4430079919211964\n", + "X : mean 0.0341427590574392, stdev 0.09887113563275718\n", + "\n", + "Ea_X 0.026648561044721943\n", + "Ea_lnR0 -0.9994904431591315\n", + "X_lnR0 -0.03999643337451283\n" + ] + } + ], + "source": [ + "# mean and standard deviation match inputs \n", + "for col in mc_inputs.columns:\n", + " print(f\"{col} : mean {mc_inputs[col].mean()}, stdev {mc_inputs[col].std()}\")\n", + "\n", + "print()\n", + "\n", + "# come up with a better way of checking \n", + "print('Ea_X', np.corrcoef(mc_inputs['Ea'], mc_inputs['X'])[0][1])\n", + "print('Ea_lnR0', np.corrcoef(mc_inputs['Ea'], mc_inputs['LnR0'])[0][1])\n", + "print('X_lnR0', np.corrcoef(mc_inputs['X'], mc_inputs['LnR0'])[0][1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Other Function Requirements \n", + "Based on the function chosen to run in the monte carlo simulation, various other data will be required. In this case we will need cell temperature and total plane of array irradiance.\n", + "\n", + "
\n", + "Please use your own API key: The block below makes an NSRDB API to get weather and meta data and then calculate cell temperature and global poa irradiance. This tutorial will work with the DEMO Key provided, but it will take you less than 3 minutes to obtain your own at https://developer.nrel.gov/signup/ so register now.) \n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Column \"relative_humidity\" not found in DataFrame. Calculating...\n", + "{'Source': 'NSRDB', 'Location ID': '1060699', 'City': '-', 'State': '-', 'Country': '-', 'Time Zone': -5, 'Dew Point Units': 'c', 'DHI Units': 'w/m2', 'DNI Units': 'w/m2', 'GHI Units': 'w/m2', 'Temperature Units': 'c', 'Pressure Units': 'mbar', 'Wind Direction Units': 'Degrees', 'Wind Speed Units': 'm/s', 'Surface Albedo Units': 'N/A', 'Version': '3.2.0', 'latitude': 25.77, 'longitude': -80.18, 'altitude': 0, 'timezone': -5}\n" + ] + } + ], + "source": [ + "weather_db = 'PSM3'\n", + "weather_id = (25.783388, -80.189029)\n", + "weather_arg = {'api_key': 'DEMO_KEY',\n", + " 'email': 'user@mail.com',\n", + " 'names': 'tmy',\n", + " 'attributes': [],\n", + " 'map_variables': True}\n", + "\n", + "weather_df, meta = pvdeg.weather.get(weather_db, weather_id, **weather_arg)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Calculate the sun position, poa irradiance, and module temperature. " + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "sol_pos = pvdeg.spectral.solar_position(weather_df, meta)\n", + "poa_irradiance = pvdeg.spectral.poa_irradiance(weather_df, meta)\n", + "temp_mod = pvdeg.temperature.module(weather_df=weather_df, meta=meta, poa=poa_irradiance, conf='open_rack_glass_polymer')\n", + "\n", + "# the function being used in the monte carlo simulation takes numpy arrays so we need to convert from pd.DataFrame to np.ndarray with .tonumpy()\n", + "poa_global = poa_irradiance['poa_global'].to_numpy()\n", + "cell_temperature = temp_mod.to_numpy()" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "# must already be numpy arrays\n", + "function_kwargs = {'poa_global': poa_global, 'module_temp': cell_temperature}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Runs monte carlo simulation for the example `pvdeg.montecarlo.vecArrhenius` function, using the correlated data dataframe created above and the required function arguments. \n", + "\n", + "We can see the necessary inputs by using the help command:" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Help on function simulate in module pvdeg.montecarlo:\n", + "\n", + "simulate(func: Callable, correlated_samples: pandas.core.frame.DataFrame, **function_kwargs)\n", + " Applies a funtion to preform a monte carlo simulation\n", + " \n", + " Parameters\n", + " ----------\n", + " func : function\n", + " Function to apply for monte carlo simulation\n", + " correlated_samples : pd.DataFrame \n", + " Dataframe of correlated samples with named columns for each appropriate modeling constant\n", + " trials : int\n", + " Number of monte carlo iterations to run\n", + " func_kwargs : dict\n", + " Keyword arguments to pass to func.\n", + " \n", + " Returns\n", + " -------\n", + " res : pandas.DataFrame\n", + " DataFrame with monte carlo results\n", + "\n" + ] + } + ], + "source": [ + "help(pvdeg.montecarlo.simulate)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Running the Monte Carlo Simulation\n", + "We will pass the target function, `pvdeg.degredation.vecArrhenius()`, its required arguments via the correlated_samples and func_kwargs. Our fixed arguments will be passed in the form of a dictionary while the randomized monte carlo input data will be contained in a DataFrame. \n", + "\n", + "All required target function arguments should be contained between the column names of the randomized input data and fixed argument dictionary, \n", + "\n", + "(You can use any data you want here as long as the DataFrame's column names match the required target function's parameter names NOT included in the kwargs)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "results = pvdeg.montecarlo.simulate(\n", + " func=pvdeg.degradation.vecArrhenius,\n", + " correlated_samples=mc_inputs, \n", + " **function_kwargs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Viewing Our Data\n", + "Let's plot the results using a histogram" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArcAAAINCAYAAAAkzFdkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACXc0lEQVR4nOzdeVxU5f4H8M8AwzIgoLKaLOaKO+KGqbkguGRp3muae5ZlaKlZZrmbS13TNtN7+6mYaJZlyzVTEbdU3HBXciFttMsSbogIzHJ+f5zmyMAwMMPALHzer9e8OHPW7zycOfOdZ57zPDJBEAQQERERETkAJ2sHQERERERkKUxuiYiIiMhhMLklIiIiIofB5JaIiIiIHAaTWyIiIiJyGExuiYiIiMhhMLklIiIiIofB5JaIiIiIHIaLtQOwBVqtFv/73/9Qq1YtyGQya4dDRERERCUIgoD79++jXr16cHIqu36WyS2A//3vfwgJCbF2GERERERUjhs3bqB+/fplLmdyC6BWrVoAxMLy9vau9P5UKhV27dqF2NhYyOXySu/P0bB8jLOF8mn2WTNk3M9AcK1g/DbpN6vEYIwtlJEtq0z53NVo0Tz1BgDgYlQIfJ0dr/Uazx/jVCoVPm78MVR3VKgVXAuTfptk7ZBsDs8h46qqfHJzcxESEiLlbWVhcgtITRG8vb0tltwqFAp4e3vzpDeA5WOcLZSPk7sToBL/WuI9YWm2UEa2rDLlU1ikht9uJQDAo3sEvF0d72OC549xKpUK7k7ucIYz3J3cbfIaYG08h4yr6vIprwmp430lJyIiIqIai8ktERERETkMJrdERERE5DAcrzEVERGRHRAEAWq1GhqNxtqh6FGpVPAM8YTKSwWPQA8UFBRYOySbo1Kp4OLigoKCApv7/9kCc8vH2dkZLi4ule6WlcktERFRNSsqKkJGRgby8/OtHUopgiCg0/udIGgFODk74dq1a9YOyeYIgoCgoCDcuHGD/eMbUJnyUSgUCA4Ohqurq9nHZ3JLRERUjbRaLa5duwZnZ2fUq1cPrq6uNpUgabVa5BTmABpA5iKDfwN/a4dkc7RaLfLy8uDl5WV0MIGaypzyEQQBRUVF+Ouvv3Dt2jU0btzY7LJlcktERBJXZydkj20vTZPlFRUVQavVIiQkBAqFwtrhlKLVaiGXySFAgJPMCe7u7tYOyeZotVoUFRXB3d2dya0B5paPh4cH5HI5/vjjD2l7czC5JSIiiYuzEx42CxSnrRyLo2NSRFSaJd4XfGcRERERkcOwmS/mS5cuxcyZM/H666/jo48+AgAUFBTgjTfewObNm1FYWIi4uDh8/vnnCAwMlLZTKpWYOHEi9u7dCy8vL4wZMwZLliyBi4vNvDQiIruh0mjheepPcTryMYBNE6qNUgnk5FTf8fz8gNDQ6jueNe3btw89e/bEnTt34OvrW+n93Lp1q0pr3nv06IG2bdtK+VB1qK5jVsdxbCIDPH78OP7973+jdevWevOnTp2Kn3/+GVu2bIGPjw8mTZqEZ599FocOHQIAaDQaDBgwAEFBQTh8+DAyMjIwevRoyOVyLF682BovhYjIrqk0Wvh9e1acbh3M5LaaKJVARARQnZ0nKBRAWlrFE9yxY8di/fr1AAAXFxfUqVMHrVu3xvDhwzF27FiHa2ZhKAnr0qULMjIy4OPjg/v371fZsbdu3VrpYWt1/6+XX34Zq1ev1lsWHx+Pzz//HGPGjEFCQoLFjlkR1XEcqye3eXl5GDFiBL744gu899570vx79+5hzZo12LRpE3r16gUAWLduHSIiInDkyBF07twZu3btwsWLF7F7924EBgaibdu2WLhwIWbMmIF58+ZVqhsJIiKi6pKTIya2iYliklvV0tKAkSPF45pSe9u3b1+sW7cOGo0GWVlZ2LFjB15//XV8++23+Omnn6r0V9OioiKrf667uroiKCgIWq22So9Tp04di+wnJCQEmzdvxooVK+Dh4QFA/FV806ZNCC3xj7fUMctTHcexenIbHx+PAQMGICYmRi+5TU1NhUqlQkxMjDSvWbNmCA0NRUpKCjp37oyUlBS0atVKr5lCXFwcJk6ciAsXLiAyMtLgMQsLC1FYWCg9z83NBSB2OqxSqSr9mnT7sMS+HBHLxzhbKx9biaM4WysjW1OZ8lGr1MWmVVDJBIvFZSusff6oVCoIggCtVislSeIfJzRtqkXbtlUfg+54Ygz6ywRBKLGuVprv6uqKgIAAAEBwcDDatm2Ljh07ok+fPli7di1efPFFAMDdu3fx5ptv4qeffkJhYSHat2+PDz/8EG3atJH2u2jRInz66ad4+PAhhg4dCj8/P+zcuRMnT54EAIwbNw53795Fhw4d8Pnnn8PNzQ3p6enYsGEDPv30U1y6dAmenp7o2bMnVqxYIcUFANu3b8e0adNw48YNdO7cGaNGjZJei1arxa1btzB58mT8+uuvuHPnDho2bIi3334bw4cPl469f/9+7N+/Hx9//DEAID09HdevX0fv3r2Rk5MDZ2dnCIKALVu2YN68ebh69SqCg4MxadIkTJs2TYrl8ccfx0svvYSrV6/i22+/Re3atfHOO+9gwoQJZf5/evXqhTZt2mDFihVm70MQBERGRuL333/Ht99+ixEjRgAAvv32W4SGhiI8PFw6Dw0ds7xy3rdvH3r37o3t27fjnXfewW+//Ybo6Ghs2rQJJ06cwBtvvIGMjAwMGDAAX3zxhdQzSMnjlKTVaiEIAlQqFZydnfWWVfQ9a9XkdvPmzTh58iSOHz9eallmZiZcXV1LtY0JDAxEZmamtE7xxFa3XLesLEuWLMH8+fNLzd+1a5dFu2VJSkqy2L4cEcvHOGuWj25EooKCAmzfvt1qcZSH55Bx5pRPLpwBiH2u7t6dDG847uhL1jp/XFxcEBQUhLy8PBQVFQEAHjxwBlALDx48QG5u1Zd5RY8nCIJeBZBarZae67Rv3x4tW7bEli1bMHToUADAkCFD4O7ujm+++Qbe3t5ISEhATEwMTpw4gdq1a+Obb77B4sWLsWzZMnTq1Albt27FZ599hrCwML3j7dmzBx4eHvjuu+8AiJVR9+/fx4wZM9C4cWP89ddfePfddzFq1Chs2bIFAHDz5k384x//wIsvvogxY8bg1KlTmDlzJgDg/v37cHJywl9//YUWLVogPj4etWrVwq5duzBmzBgEBQUhKioKCxYsQFpaGpo3by5t6+PjIw26kZeXBx8fH/z6668YNmwY3n77bQwePBjHjh3D9OnToVAo8PzzzwMQk7UPP/wQ77zzDiZPnowff/wR8fHxiIqKQuPGjQ2Wu1qtRlFRkVQW5uxD9/8aPnw41qxZg4EDBwIA/u///g/Dhg3DwYMHoVKppGOUPGZ55awri7lz52LJkiVQKBQYN24c/vGPf8DNzQ1ffPEF8vLyMGrUKCxbtgxTpkwxeJySioqK8PDhQxw4cABqtVpvWUUHPbFacnvjxg28/vrrSEpKqvY+9GbOnKn3rSo3NxchISGIjY2Ft7d3pfevUqmQlJSEPn36VEv7FXvD8jHOFsrHPd0dUAHu7u7o37+/VWIwxhbKyJZVpnxyitSYnbIHABAT0xt+rlb/gc/irH3+FBQU4MaNG/Dy8pI+/zw98fdfT1jgY6hcxo4nCAIe4iEAQCaTSZ+LcrkcLi4uBj8nmzdvjnPnzsHb2xsHDx7EyZMnkZmZCTc3NwBAZGQkfvnlF+zcuRMTJkzA2rVr8cILL2DixIkAgHbt2uHAgQPIy8vTO56npycSEhL0miO8+uqresf28fFBp06d4OTkBC8vL2zcuBENGzbEJ598AgCIiopCeno6PvjgA9SqVQve3t7w9vbGu+++K+2jdevW2L9/P7Zv346ePXvC29sbCoUCPj4+esmjrgLMy8sLAPCf//wHvXr1wsKFC6XXce3aNaxcuRKvvPIKALFrq/79+0t5R5s2bbB69WocP34cUVFRBv8/Li4ucHV1lcrCnH3o/l/jx4/HggULcOfOHQDA0aNH8c033+DIkSOQy+XSMUoes7xy1pXFokWL0Lt3bwDAiy++iHfeeQeXL1+Gv78/atWqhX/84x9ISUnBnDlzDB6npIKCAnh4eKB79+6l8sOyEuJS5VehtapAamoqsrOz0a5dO2meRqPBgQMH8Nlnn2Hnzp0oKirC3bt39Wpvs7KyEBQUBAAICgrCsWPH9PablZUlLSuLm5ub9IYrTi6XW/RCZ+n9ORqWj3G2Uj62EENZbKWMbJU55eMiPBopy0Uuh1zueMmtjrXOH41GA5lMBicnJ+kmLN29WOK8qo/B2PFKtifVxSiTyaS4DdEtO3fuHPLy8uDvrz+y2cOHD3Ht2jU4OTnh0qVLePXVV/X21bFjR+zZs0fveK1atSqV4KSmpmLevHk4c+YM7ty5I8V78+ZNNG/eHL/99puUhOl06dKl2Ot1gkajweLFi/HNN9/gzz//RFFREQoLC+Hp6am3XcnXWzw2APjtt9/wzDPP6K3TtWtXfPzxxxAEQfpZvU2bNnrrBAUFIScnx+hNeCWPbeo+dP+vwMBADBgwAF9++SUEQcCAAQMQEBBg8P9Z/Hl55axbr23bttJ0UFAQFAoFGjZsiNzcXMhkMgQFBeH48eNGy7U4JycnyGQyg+/Pir5frXbV6t27N86dO6c3b9y4cWjWrBlmzJiBkJAQyOVyJCcnY8iQIQCAS5cuQalUIjo6GgAQHR2NRYsWITs7W2oDkpSUBG9vbzRv3rx6XxAREVENlZaWhgYNGgAQf7IPDg7Gvn37Sq1najdcnroq5r89ePAAcXFxiIuLw8aNG+Hv7w+lUom4uDipiUdF/Otf/8LHH3+Mjz76CK1atYKnpyemTJli0j5MUTIpk8lkJt+UVpl9vPDCC5g0aRIAYOXKleWub0o5F49Ll5SaG6elWC25rVWrFlq2bKk3z9PTE3Xr1pXmjx8/HtOmTUOdOnXg7e2NyZMnIzo6Gp07dwYAxMbGonnz5hg1ahQ++OADZGZmYtasWYiPjzdYM0tERMa5Ojvhr+fbSdNE5dmzZw/OnTuHqVOnAhB/ms/MzISLiwvCw8MNbtO0aVMcP34co0ePluYZuv+mpN9++w23bt3C0qVLERISAgA4ceKE3joRERH46aef9OYdOXJE7/mhQ4fwzDPPYOTIkQDE2urLly/rVYy5urpCozHe/rlZs2ZS96TF992kSZNSN0NZU9++fVFUVASZTIa4uLhy169IOdsym75yrVixAk899RSGDBmC7t27IygoCFu3bpWWOzs7Y9u2bXB2dkZ0dDRGjhyJ0aNHY8GCBVaMmojIfrk4OyG/dTDyWwfDpZzkVqkETp4UH0plNQVIVlVYWIjMzEz8+eefOHnyJBYvXoxnnnkGTz31lJSoxsTEIDo6GoMGDcKuXbtw/fp1HD58GO+++66UIE2ePBlr1qzB+vXrceXKFbz33ns4e/as9HN/WUJDQ+Hq6opPP/0Uv//+O3766SepvavOK6+8gitXruDNN9/EpUuXsGnTJqkvV53GjRsjKSkJhw8fRlpaGl5++WWpWaNOeHg4jh49iuvXryMnJ8dg7eO0adOQnJyMhQsX4vLly1i/fj0+++wzTJ8+3dSirVLOzs5IS0vDxYsXK5R0V6ScbZlNNaYq+ROGu7s7Vq5cabQKPSwszKbv5iYickQlBx0wdUAAMiwtzbaPs2PHDgQHB8PFxQW1a9dGmzZt8Mknn2DMmDF67VG3b9+Od999F+PGjcNff/2FoKAgdO/eXerRaMSIEfj9998xffp0FBQUYOjQoRg7dmyp+2hK8vf3R0JCAt555x188sknaNeuHZYtW4ann35aWic0NBTfffcdpk6dik8//RQdO3bE4sWL8cILL0jrzJo1C7///jvi4uKgUCgwYcIEDBo0CPfu3ZPWmT59OsaMGYPmzZtL7YVLateuHb755hvMmTMHCxcuRHBwMBYsWICxY8eaV8BVyJQb5itSzrZMJpTs0K4Gys3NhY+PD+7du2ex3hK2b9+O/v3782YXA1g+xtlC+dRfXh9/3v8Tj9V6DDen3bRKDMbYQhnZssqUzz2NFvUuiDVY/2sRCJ8yam9PngSiosRBBwBxQIDUVKDYPcI2y9rnT0FBAa5du4YGDRpIN0vZ0ghlWq0WWWezIKgFOMmdENSm7Bu0LalPnz4ICgrChg0bquV4laHVapGbmwtvb2+HG5nNEipTPobeHzoVzddsquaWiIisq0ijhf8msRP9ogVx5Q6/Wx2jadUEoaFiopmTU33H9POzXk17fn4+Vq9ejbi4ODg7O+Orr77C7t272Xc1WQSTWyIiIhsQGlpzmnXomi4sWrQIBQUFaNq0Kb777ju9UUmJzMXkloiIiKqVh4cHdu/ebe0wyEGxoQgREREROQwmt0RERETkMJjcEhEREZHDYHJLRERERA6DN5QREZFE7uyEnH+0lqaJiOwNk1siIpLInZ3woL04ljyHxyAie8Sv5URERGRTEhIS4Ovra+0wbIJMJsMPP/wAALh+/TpkMhlOnz5t1ZhsHZNbIiKSqDVaePyWBY/fsqDWaK0dDtmQJUuWoEOHDqhVqxYCAgIwaNAgXLp0yeg2CQkJkMlkeo+SQ6qGh4fjo48+qsLI7cO8efPQtm3bUvMzMjLQr1+/6g/IjjG5JSIiSZFGi4CEEwhIOIEiJrdUzP79+xEfH48jR44gKSkJKpUKsbGxePDggdHtvL29kZGRIT3++OOPaorYMoqKiqx6/KCgILi5uVk1BnvD5JaIyIEplcDJk+JDqbR2NGTPduzYgbFjx6JFixZo06YNEhISoFQqkZqaanQ7mUyGoKAg6REYGCgt69GjB/744w9MnTpVqtktbufOnYiIiICXlxf69u2LjIyMMo+zb98+yGQy/Pzzz2jdujXc3d3RuXNnnD9/Xm+9gwcPolu3bvDw8EBISAhee+01vQQ9PDwcCxcuxOjRo+Ht7Y0JEyYAAA4dOoQePXpAoVCgbt26GDJkCO7cuQMA0Gq1WLJkCRo0aAAPDw+0adMG3377banYkpOT0b59eygUCnTp0kWq+U5ISMD8+fNx5swZqRwSEhKk8tM1SzDk/Pnz6NevH7y8vBAYGIhRo0YhJyfHyH/E8TG5JSJyUEolEBEBREWJj4gIJrhkOffu3QMA1KlTx+h6eXl5CAsLQ0hICJ555hlcuHBBWrZ161bUr18fCxYskGp2dfLz87Fs2TJs2LABBw4cgFKpxPTp08uN680338SHH36I48ePw9/fHwMHDoRKpQIApKeno2/fvhgyZAjOnj2Lr7/+GgcPHsSkSZP09rFs2TK0adMGp06dwuzZs3H69Gn07t0bzZs3R0pKCg4cOIC4uDhoNBoAYpONL7/8EqtXr8aFCxcwdepUjBw5Evv379fb77vvvosPP/wQJ06cgIuLC1544QUAwHPPPYc33ngDLVq0kMrhueeeK/e13r17F7169UJkZCROnDiBHTt2ICsrC0OHDi13W0fG3hKIiBzUrVtAfj6QmCg+HzkSyMkBQkOtGxeV1v4/7ZGZl1ntxw3yCsKJCSdM3k6r1WLKlCl44okn0LJlyzLXa9q0KdauXYvWrVvj3r17WLZsGbp06YILFy6gfv36qFOnDpydnVGrVi0EBQXpbatSqbB69Wo0bNgQADBp0iQsWLCg3Njmzp2LPn36AADWr1+P+vXr4/vvv8fQoUOxZMkSjBgxAlOmTAEANG7cGJ988gmefPJJrFq1SmoP3KtXL7zxxhvSPp9//nm0b98en3/+ufT6Q0JC4O3tjcLCQixevBi7d+9GdHQ0AODxxx/HwYMH8e9//xtPPvmktJ9FixZJz99++20MGDAABQUF8PDwgJeXF1xcXEqVgzGfffYZIiMjsXjxYmne2rVrERISgsuXL6NJkyYV3pcjYXJLROTgIiKsHQGVJzMvE3/e/9PaYVRYfHw8zp8/j4MHDxpdLzo6Wkr4AKBLly6IiIjAv//9byxcuNDotgqFQkpsASA4OBjZ2dnlxlb8eHXq1EHTpk2RlpYGADhz5gzOnj2LjRs3SusIggCtVotr164h4u83S/v27fX2efr0afzzn/80eLyrV68iPz9fSqh1ioqKEBkZqTevdevWeq8HALKzsxFq5jfOM2fOYO/evfDy8iq1LD09ncktERERWUeQV8Vr66x93EmTJmHbtm04cOAA6tevb9K2crkckZGRuHr1aoXWLU4mk0EQBJOOV1JeXh5efvllvPbaa6WWFU8wPT099ZZ5eHgY3ScA/Pzzz3jsscf0lpW8Eaz4a9K1L9Zqzb9xMy8vDwMHDsT7779fapkuea6JmNwSERFZmTlNA6qbIAiYPHkyvv/+e+zbtw8NGjQweR8ajQbnzp1D//79pXmurq5S21VLOHLkiJSo3rlzB5cvX5ZqZNu1a4eLFy+iUaNGJu2zdevWSE5Oxvz580sta968Odzc3KBUKvWaIJjKnHJo164dvvvuO4SHh8PFhSmdDm8oIyIiidzZCbeeaYFbz7Tg8LukJz4+HomJidi0aRNq1aqFzMxMZGZm4uHDh9I6o0ePxsyZM6XnCxYswK5du/D777/j5MmTGDlyJP744w+8+OKL0jrh4eE4cOAA/vzzT4vc5b9gwQIkJyfj/PnzGDt2LPz8/DBo0CAAwIwZM3D48GFMmjQJp0+fxpUrV/Djjz+WuqGspJkzZ+L48eN49dVXcfbsWfz2229Ys2YNcnJyUKtWLUyfPh1Tp07F+vXrkZ6ejpMnT+LTTz/F+vXrKxx3eHg4rl27htOnTyMnJweFhYXlbhMfH4/bt29j+PDhOH78ONLT07Fz506MGzfOol8Y7A2vXEREJJE7OyEvOhx50eFMbknPqlWrcO/ePfTo0QPBwcHS4+uvv5bWUSqVej0e3LlzBy+99BIiIiLQv39/5Obm4vDhw2jevLm0zoIFC3D9+nU0bNgQ/v7+lY5z6dKleP311xEVFYXMzEz897//haurKwCxBnb//v24fPkyunXrhsjISMyZMwf16tUzus8mTZpg165dOHPmDDp27IgnnngCv/zyi1RbunDhQsyePRtLlixBREQE+vbti59//tmk2u0hQ4agb9++6NmzJ/z9/fHVV1+Vu029evVw6NAhaDQaxMbGolWrVpgyZQp8fX3h5FRz37+swyYiIqJyVaS96759+/Ser1ixAitWrDC6TefOnXHmzBm9eWPHjsXYsWP15g0aNKhCMXTt2rVU37bFdejQAbt27Spz+fXr1w3Of/LJJ3Ho0CEAYjvZ3NxceHt7AxDbz77++ut4/fXXDW7bo0ePUrG3bdtWb56bm5te37g6xdcJDw8vtZ/GjRtj69atZb6emojJLRERSTRaAW7XbovTDeoATrJytiAisi01t86aiIhKKVRrEPTFEQR9cQSF6prbZo+I7BdrbomIiMjuGfrpn2om1twSERERkcNgcktEREREDoPJLRERERE5DLa5JSKqgZRKQNdfvp8fYObQ9kRENofJLRFRDaNUAhERQH6++FyhANLSxAT3z5uP1vvzJuD/uHViJCIyF5NbIqIaJidHTGwTE8XnI0c+qsXtEu0E+TfNpOnfjrNWl4jsC9vcEhHVUBER4kMnJwd4mOeE3CcbIvfJhniY5yQlvUS2Zt68eWjbtq21w7C669evQyaT4fTp0wDEUeJkMhnu3r1r1bisicktERERlev+/fuYMmUKwsLC4OHhgS5duuD48eN668hkMoOPf/3rX2Xud968eaXWb9asWan9/vDDD1XxsuzK2LFjMWjQIL15ISEhyMjIQMuWLa0TlA1iswQiInpEJsD1xj0AQJHMBwCH3yXRiy++iPPnz2PDhg2oV68eEhMTERMTg4sXL+Kxxx4DAGRkZOht88svv2D8+PEYMmSI0X23aNECu3fvlp67uNhXeqJSqSCXy61ybGdnZwQFBVnl2LaKNbdERCSROWsQvPIQglcegsyZw++S6OHDh/juu+/wwQcfoHv37mjUqBHmzZuHRo0aYdWqVdJ6QUFBeo8ff/wRPXv2xOOPG78z0cXFRW87Pz8/aVl4eDgAYPDgwZDJZNJznQ0bNiA8PBw+Pj4YNmwY7t+/X+ZxEhIS4Ovrix9++AGNGzeGu7s74uLicOPGDb31fvzxR7Rr1w7u7u54/PHHMX/+fKjVamm5s7Mz1qxZg2eeeQaenp5YtGgRAOC///0vOnToAHd3d/j5+WHw4MHSNoWFhZg+fToee+wxeHp6olOnTti3b1+p2Hbu3ImIiAh4eXmhb9++0heGefPmYf369fjxxx+lGu59+/aVapZgyMGDB9GtWzd4eHggJCQEr732Gh48eFDm+vaOyS0REREZpVarodFo4O7urjffw8MDBw8eNLhNVlYWfv75Z4wfP77c/V+5cgX16tXD448/jhEjRkCpVErLdE0f1q1bh4yMDL2mEOnp6fjhhx+wbds2bNu2Dfv378fSpUuNHis/Px+LFi3Cl19+iUOHDuHu3bsYNmyYtPzXX3/F6NGj8frrr+PixYv497//jYSEBCmB1Xn//fcxaNAgnDt3Di+88AJ+/vlnDB48GP3798epU6eQnJyMjh07SutPmjQJKSkp2Lx5M86ePYt//vOf6Nu3L65cuaIX27Jly7BhwwYcOHAASqUS06dPBwBMnz4dQ4cOlRLejIwMdOnSpdyyTU9PR9++fTFkyBCcPXsWX3/9NQ4ePIhJkyaVu629sq96fyIiqnbsE7fq/af9f5CXmVftx/UK8sKEExPKXa9WrVqIjo7GwoULERERgcDAQHz11VdISUlBo0aNDG6zfv161KpVC88++6zRfXfq1AkJCQlo2rQpMjIyMH/+fHTr1g3nz59HrVq14O/vDwDw9fUt9fO7VqtFQkICatWqBQAYNWoUkpOTSyWixalUKnz22Wfo1KmTFGdERASOHTuGjh07Yv78+Xj77bcxZswYAMDjjz+OhQsX4q233sLcuXOl/fzjH//AuHHj4OQk1hMOGzYMw4YNw/z586V12rRpAwBQKpVYt24dlEol6tWrB0BMVnfs2IF169Zh8eLFUmyrV69Gw4YNAYgJ8YIFCwAAXl5e8PDwQGFhoUnNEJYsWYIRI0ZgypQpAIDGjRvjk08+wZNPPolVq1aV+sLiCJjcEhFRmZRKoFu30n3ikmXlZebh/p9l/5xuCzZs2IAXXngBjz32GJydndGuXTsMHz4cqampBtdfu3YtRowYUW7y1K9fP2m6devW6NSpE8LCwvDNN9+UW+sbHh4uJbYAEBwcjOzsbKPbuLi4oEOHDtLzZs2awdfXF2lpaejYsSPOnDmDQ4cO6SXIGo0GBQUFyM/Ph0KhAIBSPTWcPn0aL730ksFjnjt3DhqNBk2aNNGbX1hYiLp160rPFQqFlNhW9PWU58yZMzh79iw2btwozRMEAVqtFteuXUNE8S5THASTWyIiKlNqatl94lYEa30rxivIy+aP27BhQ+zfvx8PHjxAbm4ugoOD8dxzzxlsT/vrr7/i0qVL+Prrr02OydfXF02aNMHVq1fLXbfkTVwymQxardbkYxaXl5eH+fPnG6xxLp6oe3p66i3z8PAwuk9nZ2ekpqbC2dlZb5mX16P/gaHXIwiCSfEbOvbLL7+M1157rdSyUAd9QzK5JSIig9wVwHvvibW13bqZltQCxkdCI30VaRpgKzw9PeHp6Yk7d+5g586d+OCDD0qts2bNGkRFRUk/y5siLy8P6enpGDVqlDRPLpdDo7HMDY5qtRonTpyQ2sNeunQJd+/elWow27Vrh0uXLpXZ3KIsrVu3RnJyMsaNG1dqWWRkJDQaDbKzs9GtWzezY3d1dTW5HNq1a4eLFy+a/HrsGW8oIyIigw4fEmtuzU1Ii4+ElpgoTnNQCPu1c+dO7NixA9euXUNSUhJ69uyJZs2alUrmcnNzsWXLFrz44osG99O7d2989tln0vPp06dj//79uH79Og4fPozBgwfD2dkZw4cPl9YJDw9HcnIyMjMzcefOnUq9DrlcjsmTJ+Po0aNITU3F2LFj0blzZynZnTNnDr788kvMnz8fFy5cQFpaGjZv3oxZs2YZ3e/cuXPx1VdfYe7cuUhLS8O5c+fw/vvvAwCaNGmCESNGYPTo0di6dSuuXbuGY8eOYcmSJfj5558rHHt4eDjOnj2LS5cuIScnByqVqtxtZsyYgcOHD2PSpEk4ffo0rly5gh9//NGhbyhjcktERBJB64S7vRvjbu/GCAt1Qrt2la9pLTkSGtmne/fuIT4+Hs2aNcPo0aPRtWtX7Ny5s9RP6Zs3b4YgCHrJaXHp6enIKfYt5+bNmxg+fDiaNm2KoUOHom7dujhy5Ih0IxkAfPjhh0hKSkJISAgiIyMr9ToUCgVmzJiB559/Hk888QS8vLz0mk/ExcVh27Zt2LVrFzp06IDOnTtjxYoVCAsLM7rfHj16YMuWLfjpp5/Qtm1b9OrVC8eOHZOWr1u3DqNHj8Ybb7yBpk2bYtCgQTh+/LhJTQNeeuklNG3aFO3bt4e/vz8OHTpU7jatW7fG/v37cfnyZXTr1g2RkZGYM2eOdGObI2KzBCIiG2WV9qpaJ9zrI9704loNhyP7MXToUAwdOrTc9SZMmIAJE8puZnH9+nW955s3by53nwMHDsTAgQP15s2bNw/z5s3TmzdlyhSpVwBjnn32WaO9OMTFxSEuLq7M5RqNBrm5uSbtVy6XY/78+Xq9KRQ3duxYjB07Vm/eoEGD9Nrc+vv7Y9euXaW2Lb5Ojx49SrXT7dChg8HtHBWTWyIiG1SR9qq8WYuIqDQmt0RENqh4e1XgUS8FugS26m7WEiDPEvtb1fp7AU4cfpeI7ItV29yuWrUKrVu3hre3N7y9vREdHY1ffvlFWt6jRw9piDnd45VXXtHbh1KpxIABA6BQKBAQEIA333xTb4g8IiJ7VlZ71aq6WUvmokG9FQdQb8UBFKg5/C45lrFjx+Lu3bvWDoOqmFVrbuvXr4+lS5eicePGEAQB69evxzPPPINTp06hRYsWAMTG07rROQBInScDYpuXAQMGICgoCIcPH0ZGRgZGjx4NuVwujfZBROTIeKMWEZE+qya3JRuHL1q0CKtWrcKRI0ek5FahUJQ5zNyuXbtw8eJF7N69G4GBgWjbti0WLlyIGTNmYN68eXB15e0QRERERDWJzbS51Wg02LJlCx48eIDo6Ghp/saNG5GYmIigoCAMHDgQs2fPlmpvU1JS0KpVKwQGBkrrx8XFYeLEibhw4UKZ3YUUFhaisLBQeq6741GlUlWoz7jy6PZhiX05IpaPcbZWPrYSR3G2VkZVQWxdJYdarXuN4rTuJRtb/qh81AbXMT79qI2tWqWCSvboruuyjmkovoq8Bmux9vmjUqmk4U8rO5pWVSh5p70txmhtujLS/R9JX2XKR6vVQhAEqFSqUqO5VfQ9a/Xk9ty5c4iOjkZBQQG8vLzw/fffo3nz5gCA559/HmFhYahXrx7Onj2LGTNm4NKlS9i6dSsAIDMzUy+xBSA9z8zMLPOYS5YsMdgVx65du/SaPVRWUlKSxfbliFg+xlmzfAoKCqS/27dvt1oc5XHkcyg93QdADxw8qOvHUpzOyLhXoeUAcPToUYPrGJ/uKm2/e3cyvPGo3W1Zxyzr+BWJ0Zqsdf64uLggKCgIeXl5KCoqskoMFSUIgsEur0h0//59a4dg08wpn6KiIjx8+BAHDhwodQ9Vvu4O2nJYPblt2rQpTp8+jXv37uHbb7/FmDFjsH//fjRv3lyvn7xWrVohODgYvXv3Rnp6Oho2bGj2MWfOnIlp06ZJz3NzcxESEoLY2Fh4e3tX6vUA4jeLpKQk9OnTp1Tn1sTyKY8tlI97ujugEsdR79+/v1ViMMYWyqiqnTol/u3a9QlpXteuT0D3g5Sx5bry6dSpk8F1ypvWiYnpDT/XRx8TZR3TUHwVeQ3WYu3zp6CgADdu3ICXlxfc3d2r/fjlEQQBD/EQACCTySzyuehoBEHA/fv3UatWLchk7FGkpMqUT0FBATw8PNC9e/dS74+KftGyenLr6uoqjXccFRWF48eP4+OPP8a///3vUuvqLtRXr15Fw4YNERQUpDf6BwBkZWUBQJntdAHAzc0Nbm5upebL5XKLXugsvT9Hw/IxzlbKxxZiKIutlFFVcHHR/ZUXmyeH7uWWtxwA5HIXg+sYn35UU+Iil0v7MHbMso5fkRityVrnj0ajgUwmg5OTE5ycbG+g0JI/I9tijNamKyPd/5H0VaZ8nJycIJPJDL4/K/p+tbn/iFar1WsPW9zp06cBAMHBwQCA6OhonDt3DtnZ2dI6SUlJ8Pb2lpo2EBFRxQlaJ9zr/jjudX8cLvzQJiI7ZNWa25kzZ6Jfv34IDQ3F/fv3sWnTJuzbtw87d+5Eeno6Nm3ahP79+6Nu3bo4e/Yspk6diu7du6N169YAgNjYWDRv3hyjRo3CBx98gMzMTMyaNQvx8fEGa2aJiGq6tLRyVtA64W5/sX8x9jdDRPbIqsltdnY2Ro8ejYyMDPj4+KB169bYuXMn+vTpgxs3bmD37t346KOP8ODBA4SEhGDIkCGYNWuWtL2zszO2bduGiRMnIjo6Gp6enhgzZoxev7hERCQOz6tQiCOdAeK0n59lBn4gIrIlVv3Nac2aNbh+/ToKCwuRnZ2N3bt3o0+fPgCAkJAQ7N+/H7du3UJBQQGuXLmCDz74oFTD9rCwMGzfvh35+fn466+/sGzZMri4WL0pMRGRTQkNFWttU1PFR9lD9Qpwvp0P59v50GoFQytQDdajRw9MnjwZU6ZMQe3atREYGIgvvvgCDx48wLhx41CrVi00atRIb7TR8+fPo1+/fvDy8kJgYCBGjRqFnGLfqnbs2IGuXbvC19cXdevWxVNPPYX09HRp+fXr1yGTybB161b07NkTCoUCbdq0QUpKSrW+drIfbFBFRFRDhIYC7dqJD8OJrTj8bv0P9qL+B3s5/G41EQA8sNLDnK8v69evh5+fH44dO4bJkydj4sSJ+Oc//4kuXbrg5MmTiI2NxahRo5Cfn4+7d++iV69eiIyMxIkTJ7Bjxw5kZWVh6NCh0v4ePHiAadOm4cSJE0hOToaTkxMGDx5c6sa2d999F9OnT8fp06fRpEkTDB8+vFRXUUSADfSWQEREVJPlA/Cy0rHzAHiauE2bNm2kJoIzZ87E0qVL4efnh5deegkAMGfOHKxatQpnz57F7t27ERkZicWLF0vbr127FiEhIbh8+TKaNGmCIUOG6O1/7dq18Pf3x8WLF9GyZUtp/vTp0zFgwAAAwPz589GiRQtcvXoVzZo1M/2Fk0NjzS0RERFVmO6mbkC896Vu3bpo1aqVNE83mFJ2djbOnDmDvXv3wsvLS3roklFd04MrV65g+PDhePzxx+Ht7Y3w8HAAgFKpLPO4ul6TiveWRKTDmlsiIiIrUkCsQbXWsU1Vsq9RXZ+kxZ8DYteeeXl5GDhwIN5///1S+9ElqAMHDkRYWBi++OIL1KtXD1qtFi1btiw1eltZxyAqicktERGRFclgetMAe9GuXTt89913CA8PN3iz961bt3Dp0iV88cUX6NatGwDg4MGD1R0mORg2SyAiIpujVAInT4qPEr9Okx2Jj4/H7du3MXz4cBw/fhzp6enYuXMnxo0bB41Gg9q1a6Nu3br4z3/+g6tXr2LPnj2YNm2atcMmO8fkloiIbIpSCUREAFFR4iMiggmuvapXrx4OHToEjUaD2NhYtGrVClOmTIGvr680/PDmzZuRmpqKli1bYurUqfjXv/5l7bDJzrFZAhERSQStDPc7hwEAXJxkVokhJwfIzwcSE8XnI0eK88rqvoyqz759+0rNu379eql5gvCok7HGjRtj69atZe4zJiYGFy9eLHP78PBwvecA4OvrW2oekQ6TWyIiekTrjNuDxO6XrD38bkSElQMgIrvEZglERERE5DBYc0tERMUIcMoTu2ASPF0BmXWaJhARmYvJLRERSWQuGoS8txsA8HBBHLxc+TFBRPaFVy0iIjuRlib+9fOzbhxERLaMyS0RkY3z8wMUCrHXAECc3rjRujEREdkq3lBGRGTjQkPFWtvUVLF7rPx8sWssIiIqjTW3RER2IDS0/H5edc0WfHyqPh5zFW9aoXs9SqWYrBefR0RkLia3RER2rnSzBRd8/LEHgoOtG1dxhppW6BLdiAixNrr4PCIic7FZAhGRnSvdbEGG3FxrD8Ggr6ymFbrRyGbNYnMLeqRHjx6YMmWK9Dw8PBwfffSR1eIh+8KaWyIiB1CRZgsVIWhlyGtXH4Dlh981FmNYmEUPRQ7m+PHj8PT0tHYYZCeY3BIR0SNaZ9wa2gaA9YffJdLx9/e3dghkR9gsgYiIkJbG9q5Uvh49emDy5MmYMmUKateujcDAQHzxxRd48OABxo0bh1q1aqFRo0b45ZdfpG3Onz+Pfv36wcvLC4GBgRg1ahRyirU/efDgAUaPHg0vLy8EBwfjww8/LHXcks0Sli9fjlatWsHT0xMhISF49dVXkZeXJy1PSEiAr68vdu7ciYiICHh5eaFv377IyMiomoIhm8LkloioBit+o9fIkYCHQoCsSA1ZkRqCIFg7vBolv0hd5qNApbH4uuZav349/Pz8cOzYMUyePBkTJ07EP//5T3Tp0gUnT55EbGwsRo0ahfz8fNy9exe9evVCZGQkTpw4gR07diArKwtDhw6V9vfmm29i//79+PHHH7Fr1y7s27cPJ0+eNBqDk5MTPvnkE1y4cAHr16/Hnj178NZbb+m/7vx8LFu2DBs2bMCBAwegVCoxffp0s1832Q82SyAiqsF0N3rpKtKcfDUYNGcnAA6/W92a/13uhvRs6o914zpKz6MW7sbDEkmsTqcGdfD1y9HS867v78XtB0Wl1ru+dIBZcbZp0wazZs0CAMycORNLly6Fn58fXnrpJQDAnDlzsGrVKpw9exa7d+9GZGQkFi9eLG2/du1ahISE4PLly6hXrx7WrFmDxMRE9O7dG4CYPNevX99oDCVvNnvvvffwyiuv4PPPP5fmq1QqrF69Gg0bNgQATJo0CQsWLDDrNZN94VWLiMgB3bxZC2lpFbshrPiNXn+VzoFsiq5PXID94lpL69atpWlnZ2fUrVsXrVq1kuYFBgYCALKzs3HmzBns3bsXXl5epfaTnp6Ohw8foqioCJ06dZLm16lTB02bNjUaw+7du7FkyRL89ttvyM3NhVqtRkFBAfLz86FQKAAACoVCSmwBIDg4GNnZ2ea9aLIrTG6JiByI2MxAwIoVUQDEJgd+flYOykKUykd94gKP+sV1lAT34oK4Mpc5yfS/qKTOjqnwugdn9KxcYCXI5XK95zKZTG+e7O/ja7Va5OXlYeDAgXj//fdL7Sc4OBhXr141+fjXr1/HU089hYkTJ2LRokWoU6cODh48iPHjx6OoqEhKbg3FyaY2NQOTWyIiBxIaCpw9q8aPPx5C165PIChI7jDJn65P3MRE8fnIkeI8R3l9ChOagFTVupbWrl07fPfddwgPD4eLS+k4GjZsCLlcjqNHjyL073/knTt3cPnyZTz55JMG95mamgqtVosPP/wQTk7irUPffPNN1b0Isju8oYyIyMGEhgING95DZKTjJH7FRUSID7J98fHxuH37NoYPH47jx48jPT0dO3fuxLhx46DRaODl5YXx48fjzTffxJ49e3D+/HmMHTtWSloNadSoEVQqFT799FP8/vvv2LBhA1avXl2Nr4psHZNbIiIiqhL16tXDoUOHoNFoEBsbi1atWmHKlCnw9fWVEth//etf6NatGwYOHIiYmBh07doVUVFRZe6zTZs2WL58Od5//320bNkSGzduxJIlS6rrJZEdYLMEIiIiqpB9+/aVmnf9+vVS84q3bW3cuDG2bt1a5j69vLywYcMGbNiwQZr35ptvGj3G1KlTMXXqVL15o0aNkqbHjh2LsWPH6i0fNGgQ29zWEExuiYhI4iyT4UGrIGmaiMjeMLklIiKJm9wZOSPEn4TdTNxWN8JZZbro4ihpRFRZTG6JiKhSio9yBjzqoquy+/Dze9SnLRFRRTG5JSKiSik+ylla2qMuuszdB/Co9pfJLRGZisktEZEN0Y3AZa2f5/OL1Aj7exjY/AVx8KxgH6nFRzkzlyX2QUTE5JaIyEYYGoHLFkcXY7tYy+Cd+0SlWeJ9weSWiMhGFB+BKyKicjdmVYWy2sWSaXTDwubn58PDw8PK0RDZlvy/v92XHD7ZFExuiYhsTEQE0K6dtaMorax2sWQaZ2dn+Pr6Ijs7GwCgUCggs6Fu17RaLVSCCgAgE2QoKCiwckS2R6vVoqioCAUFBUZHU6upzCkfQRCQn5+P7Oxs+Pr6wtnZ2ezjM7klIqIKY7tYywgKEvsS1iW4tkQQBOT+lQtBK8DJ2Ql5bnnWDsnmCIKAhw8fwsPDw6a+mNiKypSPr6+v9P4wF5NbIiKiaiaTyRAcHIyAgACoVCprh6NHpVJh7Yi1UN1RwTPQE+P2j7N2SDZHpVLhwIED6N69e6V+PndU5paPXC6vVI2tDpNbIiIiK3F2drbIh7klOTs748GNB1DdUsFJ7QR3d3drh2RznJ2doVar4e7uzuTWAGuXD5NbIiKSOMtkyG/qL00TEdkbJrdERCRxkzvjr3EdxWkrx0JEZA4mt0REZBXsL5eIqgKTWyIisjhjiSv7yyWiqmTVztlWrVqF1q1bw9vbG97e3oiOjsYvv/wiLS8oKEB8fDzq1q0LLy8vDBkyBFlZWXr7UCqVGDBgABQKBQICAvDmm29CrVZX90shInII+UVqhMzegZDZO5BfZPq1tHjiOnKk4cRV119uaqr4SEtj92JEZDlWrbmtX78+li5disaNG0MQBKxfvx7PPPMMTp06hRYtWmDq1Kn4+eefsWXLFvj4+GDSpEl49tlncejQIQCARqPBgAEDEBQUhMOHDyMjIwOjR4+GXC7H4sWLrfnSiIjslpNKY/a2FR3owVL95epqiDmgBBHpWDW5HThwoN7zRYsWYdWqVThy5Ajq16+PNWvWYNOmTejVqxcAYN26dYiIiMCRI0fQuXNn7Nq1CxcvXsTu3bsRGBiItm3bYuHChZgxYwbmzZsHV1dXa7wsIqIarToGejDUtIE1wEQEWLlZQnEajQabN2/GgwcPEB0djdTUVKhUKsTExEjrNGvWDKGhoUhJSQEApKSkoFWrVggMDJTWiYuLQ25uLi5cuFDtr4GIiKpH8aYNiYlAfv6j2mIiqtmsfkPZuXPnEB0djYKCAnh5eeH7779H8+bNcfr0abi6usLX11dv/cDAQGRmZgIAMjMz9RJb3XLdsrIUFhaisLBQep6bmwtAHFHDEiPF6PZha6PO2AqWj3G2Vj62EkdxtlZGliLeLiCHWq1CWS9Nt45GowbgYnDdypSPWqUuNq2CSiaYvA9LKl4mokflExwsPipSbsU56vljKSXLheVUGs8h46qqfCq6P6snt02bNsXp06dx7949fPvttxgzZgz2799fpcdcsmQJ5s+fX2r+rl27oFAoLHacpKQki+3LEbF8jLNm+RQUFEh/t2/fbrU4yuNo51B6ug+AHjh48BAyMu4ZXefcuXMAIo2ua0755MIZgDh4w+7dyfCG+e1vLaF4mYhKl09Fys0QRzt/qoKtXwOsjeeQcZYun/z8/AqtZ/Xk1tXVFY0aNQIAREVF4fjx4/j444/x3HPPoaioCHfv3tWrvc3KykJQUBAAICgoCMeOHdPbn643Bd06hsycORPTpk2Tnufm5iIkJASxsbHw9vau9GtSqVRISkpCnz59OCyfASwf42yhfNzT3QEV4O7ujv79+1slBmNsoYyqwqlT4t+uXZ9AZKTxdVq1alXmupUpn5wiNWan7AEAxMT0hp+rdT8mipeJTsnXXJFyK85Rzx9LUalUuACxaZ+tXgOsjeeQcVVVPrpf2stj9eS2JK1Wi8LCQkRFRUEulyM5ORlDhgwBAFy6dAlKpRLR0dEAgOjoaCxatAjZ2dkICAgAIH5L8Pb2RvPmzcs8hpubG9zcSo+9I5fLLfpPsPT+HA3LxzhbKR9biKEstlJGluLiovsrR1kvS7eOs7NLueuaUz6ucEJBgzritKsr5HJnk7a3tOJl8mie/muuSLkZ4mjnT1VhGZWN55BxVZFXVYRVk9uZM2eiX79+CA0Nxf3797Fp0ybs27cPO3fuhI+PD8aPH49p06ahTp068Pb2xuTJkxEdHY3OnTsDAGJjY9G8eXOMGjUKH3zwATIzMzFr1izEx8cbTF6JiMg4d7kzsl4WKxDcrRwLEZE5rJrcZmdnY/To0cjIyICPjw9at26NnTt3ok+fPgCAFStWwMnJCUOGDEFhYSHi4uLw+eefS9s7Oztj27ZtmDhxIqKjo+Hp6YkxY8ZgwYIF1npJRERERGRFVk1u16xZY3S5u7s7Vq5ciZUrV5a5TlhYGBu7ExEREREAG2xzS0RE1pNfpEb99/eK0zN6wtPKN5QREZmKVy0iItLj/KDI2iEQEZnNZkYoIyIiIiKqLNbcEhFVE6Xy0RCxfn7iELJERGRZTG6JiKqBUglERAC6AXYUCiAtjQkuEZGlsVkCEVE1yMkRE9vERPGRn/+oFpeIiCyHNbdERNUoIsLaERAROTYmt0REVqZri5uWVvFt/vijamJxkslQWN9HmiYisjdMbomIrMhQW1w/v7LX9/MT13nvvfLXNYe73BmZk7qK05bdNRFRtWByS0RkRcXb4kZElN+LQmioWMObk8MeF4iIDGFyS0RkIbrmBeYknRERQLt2FVs3NJRJLRFRWZjcEhFZQPHmBfbczdfDIg0eW75fnJ72JDxdna0cERGRadgVGBGRBeiaF8yaZd/dfAkQ4HL3IVzuPoQAwdrhEBGZjMktEZEFhYVZOwIiopqNzRKIiMihcJhjopqNyS0RETkMDnNMRGyWQEREDoPDHBMRa26JiMjhcJhjopqLyS0RURUyZ2hda5JBhqIAL2maiMjeMLklIqoiZQ2ta8s/k3u4OiNj2pPitJVjISIyB5NbIqIqUtbQurac3BIR2Tsmt0REVcyUoXWJiKhyzE5uVSoVMjMzkZ+fD39/f9SpU8eScRERkRU8LNIg+LOD4vSkrhx+l4jsjkldgd2/fx+rVq3Ck08+CW9vb4SHhyMiIgL+/v4ICwvDSy+9hOPHj1dVrEREVMUECHDNzoNrdh6H3yUiu1Th5Hb58uUIDw/HunXrEBMTgx9++AGnT5/G5cuXkZKSgrlz50KtViM2NhZ9+/bFlStXqjJuIiIiIqJSKtws4fjx4zhw4ABatGhhcHnHjh3xwgsvYPXq1Vi3bh1+/fVXNG7c2GKBEhHZE3vp+stesDyJqKIqnNx+9dVXFVrPzc0Nr7zyitkBERHZMz8/scuvkSPF57ruv8g8LE8iMlWleksoLCwEICa0REQkdvWVlvaouy9d919kHpYnEZnKpBvKACApKQn9+/dH7dq1oVAooFAoULt2bfTv3x+7d++uihiJiOxKaKjY9Ve7dkzELIHlSUSmMKnmdv369XjxxRfxj3/8AytWrEBgYCAAICsrC7t27UL//v2xZs0ajBo1qkqCJSKiqiWDDGpfD2nanrBdLhEBJia3ixYtwkcffYT4+PhSy8aOHYuuXbtiwYIFTG6JiOyUh6sz/ny7lzht5Vgqqqx2ubqmDLqkl00aiGoGk5JbpVKJmJiYMpf37t0bb7zxRqWDIiIiqihj7XJLJr1paUBwsHXiJKLqYVKb2xYtWmDNmjVlLl+7di2aN29e6aCIiIhMYahdri7pTU0FEhOB/PxHCTAROS6Tam4//PBDPPXUU9ixYwdiYmL02twmJyfj999/x88//1wlgRIRUdUrUGkQ9O8UcfrlaHjK7Xv43dBQNkUgqmlMSm579OiB8+fPY9WqVThy5AgyMzMBAEFBQejXrx9eeeUVhIeHV0WcRERUDbSCALeb96RpIiJ7Y3I/t+Hh4Xj//ferIhYiIiIiokoxaxAHtVqNCxcuSDW3wcHBiIiIgFwut2hwRERERESmMCm51Wq1mDNnDlauXIl79+7pLfPx8cGkSZMwf/58ODmZPDYEEREREVGlmZTcvv3220hISMDSpUsRFxdXahCH2bNno6ioiM0WiIiIiMgqTEpuv/zyS2zYsAFxcXF688PDwzFhwgSEhYVh9OjRTG6JiMhmKZVAeroPTp0CgoLYmwKRozEpub1//z7q1atX5vLg4GA8ePCg0kEREZH1aDxdrR1ClVEqgREjXJCf3wPAo4EdmOASOQ6TGsf26NED06dPR46BXrBzcnIwY8YM9OjRw1KxERFRNVO4uuDm7D64ObsPFK5m3XNs03JygPx8GaZOTUVCgpoDOxA5IJOuXKtXr0b//v0RHByMVq1a6bW5PXfuHJo3b45t27ZVSaBERLZIqRSTo7Q0a0dCpqhf/z4iItiPL5EjMim5DQkJwZkzZ7Bz5069QRw6duyIxYsXIzY2lj0lEFGNoVQCERHisK6A+BO3n591YyIiqulM/s3JyckJ/fr1Q79+/aoiHiIiuyH+xA0kJopJrp+f/bfdLFBpELj2mDj9Qke7H36XiGoes6pZjx07ho8//hgzZ87EzJkz8fHHH+P48eMm72fJkiXo0KEDatWqhYCAAAwaNAiXLl3SW6dHjx6QyWR6j1deeUVvHaVSiQEDBkChUCAgIABvvvkm1Gq1OS+NiMhkERFAu3b2n9gC4pC77tduw/3abQ6/S0R2yaSa2+zsbAwZMgSHDh1CaGioXpvbqVOn4oknnsB3332HgICACu1v//79iI+PR4cOHaBWq/HOO+8gNjYWFy9ehKenp7TeSy+9hAULFkjPFQqFNK3RaDBgwAAEBQXh8OHDyMjIwOjRoyGXy7F48WJTXh4RERER2TmTkttXX30VGo0GaWlpaNq0qd6yS5cu4YUXXkB8fDy2bNlSof3t2LFD73lCQgICAgKQmpqK7t27S/MVCgWCgoIM7mPXrl24ePEidu/ejcDAQLRt2xYLFy7EjBkzMG/ePLi6Om6XNkRERESkz6TkdufOnThw4ECpxBYAmjZtik8++aRSXYHphvStU6eO3vyNGzciMTERQUFBGDhwIGbPni3V3qakpOj13AAAcXFxmDhxIi5cuIDIyMhSxyksLERhYaH0PDc3FwCgUqmgUqnMjl9Htw9L7MsRsXyMs7XysZU4irOVMhJbP8mhVqtgSijFtxOZvg9jKlM+apW62LQKKpljNE3QlblGo4buo0+lUsPSZe8ISp431n6f2SJbuQbZqqoqn4ruz6Tk1s3NTUoEDbl//z7c3NxM2aVEq9ViypQpeOKJJ9CyZUtp/vPPP4+wsDDUq1cPZ8+exYwZM3Dp0iVs3boVAJCZmamX2AKQnut6cyhpyZIlmD9/fqn5u3bt0mvyUFlJSUkW25cjYvkYZ83yKSgokP5u377danGUx9rnUHq6D4AeOHjwEDIy7pm1ncj0fVSEOeWTC2cAMgDA7t3J8IbGojFZi67Mz507B0Cs9Dh69Ciqquwdha1fA6zN2tcgW2fp8snXdU1TDpOS2+eeew5jxozBihUr0Lt3b3h7ewMQaz6Tk5Mxbdo0DB8+3PRoAcTHx+P8+fM4ePCg3vwJEyZI061atUJwcDB69+6N9PR0NGzY0KxjzZw5E9OmTZOe5+bmIiQkBLGxsdJrqgyVSoWkpCT06dMHcrm80vtzNCwf42yhfNzT3QEV4O7ujv79+1slBmNsoYwA4NQp8W/Xrk/AwI9E5W5Xt25XaZ6p+zCmMuWTU6TG7JQ9AICYmN7wc5CBHHRl3qpVK2lep06dAFi27B2BSqXCBVwAYLvXAGuzlWuQraqq8jFWwVqcSVet5cuXQ6vVYtiwYVCr1VJ71qKiIri4uGD8+PFYtmyZycFOmjQJ27Ztw4EDB1C/fn2j6+ouRlevXkXDhg0RFBSEY8eO6a2TlZUFAGW203VzczNYwyyXyy36T7D0/hwNy8c4WykfW4ihLNYuIxcX3V85TAkjKEjsE3fsWHEHCgUQFGTaPirCnPJxEWTQ/t39l4tcDrncMZJb3f/K2fnR69G9NlP/fzWNLV8DrM3a1yBbVxV5VUWY3Cxh1apVeP/993HixAm9JDIqKsrkWk9BEDB58mR8//332LdvHxo0aFDuNqdPnwYABAcHAwCio6OxaNEiZGdnS700JCUlwdvbG82bNzcpHiKi6hAaKo5ophv21Zb6x1W4uuDGwr7itJVjISIyh1lfyb29vdGrV69KHzw+Ph6bNm3Cjz/+iFq1akltZH18fODh4YH09HRs2rQJ/fv3R926dXH27FlMnToV3bt3R+vWrQEAsbGxaN68OUaNGoUPPvgAmZmZmDVrFuLj481u/0tEVNVCQ20noSUiciQmJ7c5OTlYu3YtUlJSpGQ0KCgIXbp0wdixY+Hv71/hfa1atQoASvWwsG7dOowdOxaurq7YvXs3PvroIzx48AAhISEYMmQIZs2aJa3r7OyMbdu2YeLEiYiOjoanpyfGjBmj1y8uEREREdUMJiW3x48fR1xcHBQKBWJiYtCkSRMAYhvXTz75BEuXLsXOnTvRvn37Cu1PKGf0m5CQEOzfv7/c/YSFhfFuTiIiCyhUaeCfmCpOj4zi8LtEZHdMSm4nT56Mf/7zn1i9ejVkMpneMkEQ8Morr2Dy5MlISUmxaJBERFQ9NIIAxaW/pGkiIntjUnJ75swZJCQklEpsAUAmk2Hq1KkGB00gIrJ3SqVt3gBGRET6TEpudd1uNWvWzODyY8eOlRpQgYjI3imVQEQEoOs/XKEQezsgIiLbY1JyO336dEyYMAGpqano3bu3lMhmZWUhOTkZX3zxhVn93BIR2bKcHDGxTUwUn48c+agWl4iIbItJyW18fDz8/PywYsUKfP7559BoxGEZnZ2dERUVhYSEBAwdOrRKAiUisraICGtHQERE5TG5K7DnnnsOzz33HFQqFXL+rrrw8/PjCB1EREREZHVmj6sol8ulUcKIiIiIiGyBycntmTNn8N///hd16tTB0KFD4efnJy3Lzc3FlClTsHbtWosGSURE1UPh6oI/lg4Qp60cS1X44w9rR0BEVc3JlJV37dqFjh07YvPmzXj//ffRrFkz7N27V1r+8OFDrF+/3uJBEhERVYafn9jLxXvvAQqFAG/vImuHRERVxKTkdt68eZg+fTrOnz+P69ev46233sLTTz+NHTt2VFV8RERElRYaKnbflpoKnD2rhr//Q2uHRERVxKRmCRcuXMCGDRsAiIM2vPXWW6hfvz7+8Y9/YPPmzejQoUOVBElERNWjUKWB3zenxemhbR1q+N3QUPGhUgHnz1s7GiKqKiYlt25ubrh7967evOeffx5OTk547rnn8OGHH1oyNiIiqmYaQYDnuUxx+p8cfpeI7I9JyW3btm2xd+9eREVF6c0fNmwYBEHAmDFjLBocEREREZEpTEpuJ06ciAMHDhhcNnz4cAiCgC+++MIigRERERERmcqk5Hbw4MEYPHhwmcuff/55PP/885UOioiIiIjIHCb1lkBEREREZMuY3BIRERGRwzB7+F0ioposLc3aERARkSEVTm5zc3Ph7e1dlbEQEdk83UhXI0eKzxUKcZ6j8JA7Q7kgTpquSZRKICdHnPbzE/vEJSL7U+Hktnbt2sjIyEBAQAB69eqFrVu3wtfXtwpDIyKyPbqRrhw1CZLJZBBcxY8GmZVjqU5KJRARAeTni88VCvH/7Ej/W6KaosLJrZeXF27duoWAgADs27cPKpWqKuMiIqo2ptbY6Ua6IseRkyMmtomJ4vORI8V5/D8T2Z8KJ7cxMTHo2bMnIiIiAIjdgrm6uhpcd8+ePZaJjoioirHGTl+RWoO6W8WxaYuebQlPl5rVNOHvjzgismMVTm4TExOxfv16pKenY//+/WjRogUUCkVVxkZEVOVYY6dPrRXgdfKmOD2ohZWjISIyXYWTWw8PD7zyyisAgBMnTuD9999nm1sichjFa+x0PSE4WntaIqKawKyuwPbu3StNC4IAQLwJgYjInhnqCYFdfhER2RezB3H48ssv0apVK3h4eMDDwwOtW7fGhg0bLBkbEVG10vWEkJoqNlPIz390oxk5rrQ0fokhciRm1dwuX74cs2fPxqRJk/DEE08AAA4ePIhXXnkFOTk5mDp1qkWDJCKqLiV7QmDS47jq1jXcZzG/0BDZN7OS208//RSrVq3C6NGjpXlPP/00WrRogXnz5jG5JSK75+iDNVDZfRYzuSWyb2YltxkZGejSpUup+V26dEFGRkalgyIisjZHH6yBROyzmMjxmJXcNmrUCN988w3eeecdvflff/01GjdubJHAiIisrSYmPh5yZ9yYFSNNExHZG7OS2/nz5+O5557DgQMHpDa3hw4dQnJyMr755huLBkhERNVHJpNB6+UmTls5FiIic5jVW8KQIUNw9OhR+Pn54YcffsAPP/wAPz8/HDt2DIMHD7Z0jEREREREFWJWzS0AREVFIVE3pA8RETmEIrUGdbaJXUQUPRVR44bfJSL7Z3Y/t0RE5HjUWgG1jvyBWkf+gForWDscIiKTMbklIiIiIofB5JaIiIiIHIbZbW6JiIgcmW50OvZxTGRfzKq5feGFF3D//v1S8x88eIAXXnih0kERERFZS/HR6aKigIgIQKm0dlREVFFmJbfr16/Hw4cPS81/+PAhvvzyy0oHRUREZC260elSU4HERCA/n0PyEtkTk5ol5ObmQhAECIKA+/fvw93dXVqm0Wiwfft2BAQEWDxIIiKi6lQTR6cjchQmJbe+vr6QyWSQyWRo0qRJqeUymQzz58+3WHBERFS93F2ccfOtntI0EZG9MSm53bt3LwRBQK9evfDdd9+hTp060jJXV1eEhYWhXr16Fg+SiMjSlErxp2bdTUMkcnKSQVNHIU5bORYiInOYlNw++eSTAIBr164hNDQUMhlHHici+6NUijcJ5eeLzxUK8SYiIiKyf2Z1BfbHH3/gjz/+KHN59+7dzQ6IiKiq5eSIiW1iopjksqunR4rUWvjuuiROxzaFpwvrb4nIvpiV3Pbo0aPUvOK1uBqNxuyAiIiqS0QE0K6dtaOwLWqtFj4HfhenYxqDjROIyN6YddW6c+eO3iM7Oxs7duxAhw4dsGvXrgrvZ8mSJejQoQNq1aqFgIAADBo0CJcuXdJbp6CgAPHx8ahbty68vLwwZMgQZGVl6a2jVCoxYMAAKBQKBAQE4M0334RarTbnpRERERGRHTOr5tbHx6fUvD59+sDV1RXTpk1Dampqhfazf/9+xMfHo0OHDlCr1XjnnXcQGxuLixcvwtPTEwAwdepU/Pzzz9iyZQt8fHwwadIkPPvsszh06BAAsZZ4wIABCAoKwuHDh5GRkYHRo0dDLpdj8eLF5rw8IiKiculuSgTYtIXIllh0+N3AwMBSNa/G7NixQ+95QkICAgICkJqaiu7du+PevXtYs2YNNm3ahF69egEA1q1bh4iICBw5cgSdO3fGrl27cPHiRezevRuBgYFo27YtFi5ciBkzZmDevHlwdXW15EskIiIyeFNiWhoTXCJbYFZye/bsWb3ngiAgIyMDS5cuRdu2bc0O5t69ewAgdTGWmpoKlUqFmJgYaZ1mzZohNDQUKSkp6Ny5M1JSUtCqVSsEBgZK68TFxWHixIm4cOECIiMjSx2nsLAQhYWF0vPc3FwAgEqlgkqlMjt+Hd0+LLEvR8TyMc7WysdW4iiusmUktlqSQ61WwQZfXqVVpnzUKnWxaRVUMsFicdkKU8qnrHMlMxPIz5cjIUEsr7FjXZCZqUJwcBUEXM1KlostXgOszdau07amqsqnovszK7lt27YtZDIZBEH/ote5c2esXbvWnF1Cq9ViypQpeOKJJ9CyZUsAQGZmJlxdXeHr66u3bmBgIDIzM6V1iie2uuW6ZYYsWbLE4GATu3btgkKhMCt+Q5KSkiy2L0fE8jHOmuVTUFAg/d2+fbvV4iiPuWWUnu4DoAcOHjyEjIx7lg3KhphTPrlwBiDeILx7dzK84bg3CFekfMo6V3Tzb906+PccxzyfbP0aYG38HDPO0uWTr/uppBxmJbfXrl3Te+7k5AR/f3+94XhNFR8fj/Pnz+PgwYPlr1xJM2fOxLRp06Tnubm5CAkJQWxsLLy9vSu9f5VKhaSkJPTp0wdyubzS+3M0LB/jbKF83NPdARXg7u6O/v37WyUGYypbRqdOiX+7dn0CBn7csXuVKZ+cIjVmp+wBAMTE9Iafq0Vbr9kEU8qnrHOl+HwdRzmfVCoVLuACANu9BlibLVynbVlVlY/ul/bymHXVCgsLM2ezMk2aNAnbtm3DgQMHUL9+fWl+UFAQioqKcPfuXb3a26ysLAQFBUnrHDt2TG9/ut4UdOuU5ObmBjc3t1Lz5XK5Rf8Jlt6fo2H5GGcr5WMLMZTF3DJycdH9lcOGX16lmVM+Xs4u+N9Usa9yLw93yJ0cd7CeipRPWedK8fmP5jnm+WTL1wBrs5XrtK2qiryqIszuwHD//v0YOHAgGjVqhEaNGuHpp5/Gr7/+atI+BEHApEmT8P3332PPnj1o0KCB3vKoqCjI5XIkJydL8y5dugSlUono6GgAQHR0NM6dO4fs7GxpnaSkJHh7e6N58+bmvjwiohrJyUkGVWAtqAJrwcmBE1siclxmJbeJiYmIiYmBQqHAa6+9htdeew0eHh7o3bs3Nm3aVOH9xMfHIzExEZs2bUKtWrWQmZmJzMxMPHz4EIDY5dj48eMxbdo07N27F6mpqRg3bhyio6PRuXNnAEBsbCyaN2+OUaNG4cyZM9i5cydmzZqF+Ph4g7WzREREROS4zGqWsGjRInzwwQeYOnWqNO+1117D8uXLsXDhQjz//PMV2s+qVasAlB7xbN26dRg7diwAYMWKFXBycsKQIUNQWFiIuLg4fP7559K6zs7O2LZtGyZOnIjo6Gh4enpizJgxWLBggTkvjYioRitSa+Gz96o43bMRh98lIrtjVnL7+++/Y+DAgaXmP/3003jnnXcqvJ+SvS0Y4u7ujpUrV2LlypVlrhMWFsa7OYmILECt1cI3+Yo4/eTj4PC7RGRvzLpqhYSE6LWD1dm9ezdCQkIqHRQRERERkTnMqrl944038Nprr+H06dPo0qULAODQoUNISEjAxx9/bNEAiYiIiIgqyqzkduLEiQgKCsKHH36Ib775BgAQERGBr7/+Gs8884xFAyQiIiIiqiize+cePHgwBg8ebMlYiIiIbFJamvjXzw8IDbVuLERkXIWTW0EQIJOxz0MiIqo5/PwAhQIYOVJ8rlA8SnSJyDZV+IayFi1aYPPmzSgqKjK63pUrVzBx4kQsXbq00sERERFZU2iomMympgKJiUB+PpCTY+2oiMiYCtfcfvrpp5gxYwZeffVV9OnTB+3bt0e9evXg7u6OO3fu4OLFizh48CAuXLiASZMmYeLEiVUZNxERVQE3F2dkxD8hTZOY4LIpApH9qHBy27t3b5w4cQIHDx7E119/jY0bN+KPP/7Aw4cP4efnh8jISIwePRojRoxA7dq1qzJmIiKqIs5OMhSF+IrT1g2FiMgsJt9Q1rVrV3Tt2rUqYiEiIiIiqhSze0sgIiLHU6TWwvvQNXH6iQYcfpeI7A6TWyKqMZRK8WYg3u1eNrVWi9q//CZOR4eBw++WxvOHyLYxuSWiGkGpBCIixLvdAbFLJz8/68ZE9sVQt2B+fuw9gcjWMLklohohJ0dMbBMTxSSXnfGTqXTdgumSWd05xOSWyLYwuSWiGiUiAmjXztpRkL1it2BEts+sxlQnT57EuXPnpOc//vgjBg0ahHfeeafcQR6IiIiIiKqKWcntyy+/jMuXLwMAfv/9dwwbNgwKhQJbtmzBW2+9ZdEAiYiIiIgqyqzk9vLly2jbti0AYMuWLejevTs2bdqEhIQEfPfdd5aMj4iIiIiowsxqcysIArRaLQBg9+7deOqppwAAISEhyGHLeiIiu+Xm4ozMlzpL00RE9sas5LZ9+/Z47733EBMTg/3792PVqlUAgGvXriEwMNCiARIRUfVxdpKhsGFdcdrKsdgrXX/K7JGDyDrMSm5XrFiBkSNH4ocffsC7776LRo0aAQC+/fZbdOnSxaIBEhER2Yvi/SkrFGLXYUxwiaqXWcltmzZt9HpL0PnXv/4FFxf2LkZEZK9UGi28jinF6Y6hgDNHKDOFrj/lWbOA994TnzO5JapeZl21Hn/8cdy6davU/IKCAjRp0qTSQRERkXWoNFrU/fEC6v54ASqN1trh2K2wMGtHQFRzmZXcXr9+HRqNptT8wsJC3Lx5s9JBERERERGZw6Q2BD/99JM0vXPnTvj4+EjPNRoNkpOT0aBBA8tFR0RERERkApOS20GDBgEAZDIZxowZo7dMLpcjPDwcH374ocWCIyIiIiIyhUnJra5v2wYNGuD48ePw8/OrkqCIiCpK1+0SwK6XiIjIzN4Srl27Zuk4iIhMplQCrVuLd6cD7HqJiIjMTG4BIDk5GcnJycjOzpZqdHXWrl1b6cCIiMpz65aY2CYmis9HjmTXS0RENZ1Zye38+fOxYMECtG/fHsHBwZDJZJaOi4iowiIirB2B43B1dkL22PbSNBGRvTEruV29ejUSEhIwatQoS8dDRGQRJdviUsW4ODvhYTNxGHUOyUNE9sisa1dRURGH2SUim1V8CFRAbIu7caN1YyIiouph1m9OL774IjZt2mTpWIiILEI3BGpiovjIz39Ui0vGqTRaeJ64Ac8TNzhCGRHZJbNqbgsKCvCf//wHu3fvRuvWrSGXy/WWL1++3CLBERFVRvG2uH/8Yb047IlKo4Xft2fF6dbBANvdVlhamrUjICLAzOT27NmzaNu2LQDg/Pnzest4cxkR2RI/P7FZwnvviX/Z/pYsTXeOjRwpPud5RmRdZiW3e/futXQcREQWUbL2LDRUnJeTw0EeqGoUP8cA8TxjMxgi6+HNsETkEIzVnoWGMqmlqlXyHGNyS2Q9ZiW3PXv2NNr8YM+ePWYHRERkDkO1Z0xoiYhqHrOSW117Wx2VSoXTp0/j/PnzGDNmjCXiIiIyGWtoiYjIrOR2xYoVBufPmzcPeXl5lQqIiIiIiMhcFu3jZeTIkVi7dq0ld0lERNXI1dkJfz3fDn89347D7xKRXbLoDWUpKSlwd3e35C6JiKgauTg7Ib91sDht5ViIiMxh1rXr2Wef1XsuCAIyMjJw4sQJzJ492yKBERERERGZyqzk1sfHR++5k5MTmjZtigULFiA2NtYigRERUfVTa7RQXMgSp1sEcoQyIrI7ZiW369ats3QcRERkA4o0WvhvOilOL4hjcltJukFF2DUdUfWp1FUrNTUViYmJSExMxKlTp0ze/sCBAxg4cCDq1asHmUyGH374QW/52LFjIZPJ9B59+/bVW+f27dsYMWIEvL294evri/Hjx7PHBiIisqrig4pERQEREYBSae2oiGoGs2pus7OzMWzYMOzbtw++vr4AgLt376Jnz57YvHkz/P39K7SfBw8eoE2bNnjhhRdKtePV6du3r15NsZubm97yESNGICMjA0lJSVCpVBg3bhwmTJiATZs2mfPSiIiIKq34oCJpaWKSm5PD2lui6mBWcjt58mTcv38fFy5cQEREBADg4sWLGDNmDF577TV89dVXFdpPv3790K9fP6PruLm5ISgoyOCytLQ07NixA8ePH0f79u0BAJ9++in69++PZcuWoV69eia8KiIiIsvhoCJE1mFWs4QdO3bg888/lxJbAGjevDlWrlyJX375xWLBAcC+ffsQEBCApk2bYuLEibh165a0LCUlBb6+vlJiCwAxMTFwcnLC0aNHLRoHEREREdk+s2putVot5HJ5qflyuRxarbbSQen07dsXzz77LBo0aID09HS888476NevH1JSUuDs7IzMzEwEBATobePi4oI6deogMzOzzP0WFhaisLBQep6bmwtAHEZYpVJVOm7dPiyxL0fE8jHO1srHVuIo7lEZqQHIoVarYINhWk1lziG1Sl1sWgWVTLBYXLaiut9jajVgT+dpyXKxxWuAtdnaddrWVFX5VHR/ZiW3vXr1wuuvv46vvvpK+un/zz//xNSpU9G7d29zdmnQsGHDpOlWrVqhdevWaNiwIfbt21ep4yxZsgTz588vNX/Xrl1QKBRm77ekpKQki+3LEbF8jLNm+RQUFEh/t2/fbrU4yiP+QtMDBw8eQkbGPWuHY3PMOYdy4QxABgDYvTsZ3tBYOCrbUV3vsfR0H9jreWrr1wBr4+eYcZYun/z8/AqtZ1Zy+9lnn+Hpp59GeHg4QkJCAAA3btxAy5YtkZiYaM4uK+Txxx+Hn58frl69it69eyMoKAjZ2dl666jVaty+fbvMdroAMHPmTEybNk16npubi5CQEMTGxsLb27vScapUKiQlJaFPnz4Ga7hrOpaPcbZQPu7p7oAKcHd3R//+/a0SgzG6MurUqRMAoGvXJxAZaeWgbEhlzqG7Gi0mPiZeV/u2DICvA3YFVt3vMV1nQvZynqpUKlzABQC2ew2wNlu4Ttuyqiof3S/t5TEruQ0JCcHJkyexe/du/PbbbwCAiIgIxMTEmLO7Crt58yZu3bqF4GBxaMjo6GjcvXsXqampiIqKAgDs2bMHWq1W+tAzxM3NrVSvC4DYrMKS/wRL78/RsHyMs5XysYUYyiKXi5cwFxc5bDhMqzHnHFLIgQftxUoLBQBHLtbqeo+5/P1Je+WKHC4u9tfnrS1fA6zNVq7Ttqoq8qqKMHvocJlMhj59+qBPnz7m7gJ5eXm4evWq9PzatWs4ffo06tSpgzp16mD+/PkYMmQIgoKCkJ6ejrfeeguNGjVCXFwcADGh7tu3L1566SWsXr0aKpUKkyZNwrBhw9hTAhER2YTifd4C4nRamn0luET2xKTfm/bs2YPmzZsbrBa+d+8eWrRogV9//bXC+ztx4gQiIyMR+ffvNNOmTUNkZCTmzJkDZ2dnnD17Fk8//TSaNGmC8ePHIyoqCr/++qterevGjRvRrFkz9O7dG/3790fXrl3xn//8x5SXRUREf1NrtPD4LQsev2VBrbHcDcI1ma7P29RUIDERyM8X+7wloqphUs3tRx99hJdeeslgu1QfHx+8/PLLWL58Obp161ah/fXo0QOCUPaduDt37ix3H3Xq1OGADUREFlKk0SIg4YQ4zeF3LYZ93hJVH5OuWmfOnCk1/G1xsbGxSE1NrXRQRERERETmMCm5zcrKMtqY18XFBX/99VelgyIiIiIiModJye1jjz2G8+fPl7n87NmzUk8GRERERETVzaTktn///pg9e7bUwXtxDx8+xNy5c/HUU09ZLDgiIiIiIlOYdEPZrFmzsHXrVjRp0gSTJk1C06ZNAQC//fYbVq5cCY1Gg3fffbdKAiUiIiIiKo9JyW1gYCAOHz6MiRMnYubMmVJPBzKZDHFxcVi5ciUCAwOrJFAiIiIiovKYPIhDWFgYtm/fjjt37uDq1asQBAGNGzdG7dq1qyI+IiKqRnJnJ9x6poU0TURkb8weoax27dro0KGDJWMhIiIrkzs7IS86XJy2bihERGbh13IiIiIichhm19wSEZHj0WgFuF27LU43qAM4yawcERGRaZjcEpHdUSqB9HQf3L3LxMvSCtUaBH1xRJxeEAe48mOCiOwLr1pEZFeUSqB1axfk5/cAACgUgJ+fdWMiIiLbwTa3RGRXcnKA/HwZpk5NxdGjKqSlAaGh1o6KiIhsBWtuicgu1a9/H5GRgJy39BMRUTFMbomIiKxEqRR/jQDE5jX8FYKo8pjcEhERWYFSCUREAPn54nOFAmxmQ2QBTG6JyGaxVoscmdh+HEhMFJ+PHCnO43lOVDlMbonIJpVVq0VVy8XJCXf6NZOmqepFRFg7AiLHwuSWiGxSWbVaVLVcXZyQ+2RDcdrKsRARmYPJLRHZNNZqERGRKZjcEhGRRKMV4PrnPXH6MR8Ov0tEdocNqoiISFKo1iB45SEErzyEQrXG2uEQEZmMyS0REREROQw2SyAiu8HeEoiIqDxMbonI6srrz9bPT+wKbORI8blCIcDbu6h6gySqRuzjmch8TG6JyKoq0p9taKg4T/dh7+OjxvnzD6s3UKJqwpHLiCqHyS0RWVVF+7MNDX304a5SAefPV1+MRJZmrIkNRy4jqhwmt0RkE4r3Z8u2teSoSjexEecZ+kLHPp6JzMPklohsRlkf/FR9XJyccLd3Y2maLKtkExtde1qOvkdkOUxuichmlPXBT9XH1cUJ9/o0EaetHIujKt7Ehogsj8ktEdkUfvATEVFlMLklIiKJVitA/leeOO3vxeF3icjusEEVERFJCtQa1FtxAPVWHEABh98lIjvE5JaIiIiIHAaTWyIiIiJyGExuiYiIiMhhMLklIiIiIofB5JaIiIiIHAa7AiMiIrIRZQ09rZvPgU2IysfkloiIJC5OTrjX/XFpmqqHsaGnS85PS2OCS2QMk1siIpK4ujjhbv8IcdrKsdQkxoae1s1PSxOT3JwcJrdExjC5JSIisgFlDT3NIamJTMPkloiIJFqtAOe7D8VpXw8Ov0tEdocNqoiISFKg1qD+B3tR/4O9HH6XiOwSk1siIiIichhWTW4PHDiAgQMHol69epDJZPjhhx/0lguCgDlz5iA4OBgeHh6IiYnBlStX9Na5ffs2RowYAW9vb/j6+mL8+PHIy8urxldBRERUfdLSgJMnAaXS2pEQ2SarJrcPHjxAmzZtsHLlSoPLP/jgA3zyySdYvXo1jh49Ck9PT8TFxaGgoEBaZ8SIEbhw4QKSkpKwbds2HDhwABMmTKiul0BEZVAqxQ9gfggTWUbx7sKiooCICL63iAyx6g1l/fr1Q79+/QwuEwQBH330EWbNmoVnnnkGAPDll18iMDAQP/zwA4YNG4a0tDTs2LEDx48fR/v27QEAn376Kfr3749ly5ahXr161fZaiOgRpVL84M3PF5+zb06iyiveXRi7BSMqm832lnDt2jVkZmYiJiZGmufj44NOnTohJSUFw4YNQ0pKCnx9faXEFgBiYmLg5OSEo0ePYvDgwQb3XVhYiMLCQul5bm4uAEClUkGlUlU6dt0+LLEvR8TyMc7WysecODIzgfx8ORIS1ACAsWNdkJmpQnBw6XXVagCQQ61WoaKHsrUysjWVKR+1Sl1sWgWVTLBYXLbCns+f4GDxYc77pqJKlos9llNVs+dzqDpUVflUdH82m9xmZmYCAAIDA/XmBwYGSssyMzMREBCgt9zFxQV16tSR1jFkyZIlmD9/fqn5u3btgkKhqGzokqSkJIvtyxGxfIyzZvnomv4UFBRg+/btJm+fnu4DoAdu3Tr495weOHjwEDIy7pW5blnLjeE5ZJw55ZMLZwBi91+7dyfDG47bY4I9nz+6983mzWdx9uwt+Ps/rJLjmHsNqCns+RyqDpYun3zdz4HlsNnktirNnDkT06ZNk57n5uYiJCQEsbGx8Pb2rvT+VSoVkpKS0KdPH8jl8krvz9GwfIyzhfJxT3cHVIC7uzv69+9v8vanTol/u3Z9QprXtesTiIw0vq6h5YbYQhnZssqUzx21FlME8cbduNjGqO3ieJ3qOML5o1QCs2cLWLEiCgqFgLNn1RZrnqBSqXABFwCYfw1wdI5wDlWlqiof3S/t5bHZ5DYoKAgAkJWVheBiv2VmZWWhbdu20jrZ2dl626nVaty+fVva3hA3Nze4ubmVmi+Xyy36T7D0/hwNy8c4Wykfc2JwcdH9lRebJ4ehXRVf19RD2UoZ2SpzysdTDtwe1FKcBuDIpWvP50/DhmK7219/BUaOlOHePdPfPxVlr2VUHez5HKoOVZFXVYTNfiVv0KABgoKCkJycLM3Lzc3F0aNHER0dDQCIjo7G3bt3kZqaKq2zZ88eaLVadOrUqdpjJiIiqi6hoeKNm0Skz6o1t3l5ebh69ar0/Nq1azh9+jTq1KmD0NBQTJkyBe+99x4aN26MBg0aYPbs2ahXrx4GDRoEAIiIiEDfvn3x0ksvYfXq1VCpVJg0aRKGDRvGnhKIiMwgCAKcHhSJ056ugIzD7xKRfbFqcnvixAn07NlTeq5rBztmzBgkJCTgrbfewoMHDzBhwgTcvXsXXbt2xY4dO+Du7i5ts3HjRkyaNAm9e/eGk5MThgwZgk8++aTaXwsRkSN4qNIg5L3d4vSCOHi52mzrNSIig6x61erRowcEoexuZmQyGRYsWIAFCxaUuU6dOnWwadOmqgiPiIiIiOwMv5ITUbVIS9N/7ufHzueJLEX3/uL7iojJLRFVseJDhhanG7WMiMxX8v3F0QCJmNwSURUrPmSoTvGhQ4nIfBUdklepfPR+Y+0uOTomt0RU5UJD+WFKVFXKe38plWKXYbrBnVi7S47OZvu5JSIiosrLyRET28RE8ZGfz19NyLGx5paIiCQuTjLktasvTZPj4IAPVFMwuSUiIomrizNuDW0jTls5FiIiczC5JaJK4Y0qRERkS5jcEpHZyrpRheyXIAiQqTTitNyZw+8Skd1hcktEZit+owrA7r0cwUOVBqFzdorTHH6XiOwQr1pEVGm8UYWIiGwFuwIjIiIiIofB5JaIiIiIHAabJRARETkQ3U2d7L2Eaiomt0RERA7Az0/ssWTkSPE5ey+hmorNEoiIiBxAaKiYzKamcphdqtlYc0tERBJnmQwPWgVJ02RfQkP1myKw5pZqIia3REQkcZM7I2dElDht5VjIfIaaKPj5sSaXagYmt0RUIRxml8h+6JoolHzPMrmlmoDJLRGVi8PsEtmfkk0UiGoKJrdEVC5Thtk1Jellgmx78ovUCPt7+N38BXHw5PC7DondhZEj41WLiCrM2DC7ZbXxq+y6RGQ5ZXUXxgSXHAmTWyKyiLLa+FV2XSKynOLvvbS0R7/C8P1HjoTJLRFZjClt/NgekMg6+N4jR8dBHIiIiIjIYTC5JSIiIiKHweSWiIiIiBwG29wSEZHEWSZDflN/aZqIyN4wuSUiIomb3Bl/jesoTls5FiIic7BZAhERERE5DCa3REREROQw2CyBiIgk+UVqhCzcLU7PjuHwu0Rkd3jVIiIiPU4qjbVDICIyG5slEBEREZHDYHJLRERERA6DzRKIiIioTEolkJMjTvv5AaGh1o2HqDxMbomIiMggjQaIiADy88XnCgWQlsYEl2wbmyUQkR6lElCpxGmVCjh5UvwwI6KaR6sVE9vERPGRn/+oFpfIVrHmlogkSuXftTSvAPAGsrOBqChxmUIh/iTJDzbH5iSToaBBHWmaaiaNRvzf677oRkRYMRgiEzG5JSJJTo5YM1O7NnBHAwQEAL+kist0be2Y3Do2d7kzsl6OFqetHAtZh1IJ3L3rBm8U4c4dfrEl+8PklohKkcsBaMS/7doZXodNFYgc061bj6YDAoC04/xiS/aFyS0RmcTPT6zJGTlSfK6r1SEixyOX8+Yxsj9MbokIgPhTZEVqY0NDxfXYNZBjyi9So/77e8XpGT05/C4R2R1etYjo0Y1k+WJNrFM5/aiEhjKhdWTOD4qsHQIRkdnYFRgRSTeSJSaKtbLOztaOiIiIyDw2ndzOmzcPMplM79GsWTNpeUFBAeLj41G3bl14eXlhyJAhyMrKsmLERPYtIoI1skREZN9sOrkFgBYtWiAjI0N6HDx4UFo2depU/Pe//8WWLVuwf/9+/O9//8Ozzz5rxWiJiIiIyJpsvs2ti4sLgoKCSs2/d+8e1qxZg02bNqFXr14AgHXr1iEiIgJHjhxB586dqztUIiIiIrIym09ur1y5gnr16sHd3R3R0dFYsmQJQkNDkZqaCpVKhZiYGGndZs2aITQ0FCkpKUaT28LCQhQWFkrPc3NzAQAqlQoq3XAslaDbhyX25YhYPsZZo3zUagCQQ61WoeRhbfH/xHPIuMqUj1qlLjatgkomWCwuW8HzR6R73587p4ZaLaBuXbFZkqrYOQA8Kidj14mahueQcVVVPhXdn00nt506dUJCQgKaNm2KjIwMzJ8/H926dcP58+eRmZkJV1dX+Pr66m0TGBiIzMxMo/tdsmQJ5s+fX2r+rl27oFAoLBZ/UlKSxfbliFg+xlVn+aSn+wDogYMHDyEj4x4KCgoAiO3at2/fXm1xmIrnkHHmlM99OKOwvi8AIHl3MmpBY+GobEdNP3/++ssDbm69MHasmAq4uanx2Wd7kJvrKq1T/BpQ8jpBPIfKY+nyyc/Pr9B6Np3c9uvXT5pu3bo1OnXqhLCwMHzzzTfw8PAwe78zZ87EtGnTpOe5ubkICQlBbGwsvL29KxUzIH6zSEpKQp8+fSCXyyu9P0fD8jGuOstHqRRHI7p7VxxHvmvXJxAZCbinuwMqwN3dHf3796/SGMzBc8i4ypTPAwCj/t7mKZUKnlUQn7Xx/HmkZ08Bt26pkJYmw9ixLmjduidUKjW24iQA/WvAqVPiNrrrRE3Gc8i4qiof3S/t5bHp5LYkX19fNGnSBFevXkWfPn1QVFSEu3fv6tXeZmVlGWyjW5ybmxvc3NxKzZfL5Rb9J1h6f46G5WNcVZePUgm0bi12AQaI/dsGBclR8pC2/D/iOWScOeVTfG25XA5HLl2eP0DDhuLD5e9swMWldHnoyqj4OjW82CQ8h4yriryqImy+t4Ti8vLykJ6ejuDgYERFRUEulyM5OVlafunSJSiVSkRHR1sxSiL7ULxv29RUsX9bdgNGRET2zqZrbqdPn46BAwciLCwM//vf/zB37lw4Oztj+PDh8PHxwfjx4zFt2jTUqVMH3t7emDx5MqKjo9lTApEJIiKAdu2sHQXZiodFGjy2fL84Pe1JeLpyRA8isi82ndzevHkTw4cPx61bt+Dv74+uXbviyJEj8Pf3BwCsWLECTk5OGDJkCAoLCxEXF4fPP//cylETEdkvAQJc7j6UpomI7I1NJ7ebN282utzd3R0rV67EypUrqykiIiIix5WWBqjVMmuHQVQpNp3cEhERUdXz8xNvKh05EgBcMK2c9ZVKsd2+blu21ydbwuSWiIiohgsNFWttc3IAtVqFnXGF0N41vK5SCXTrpt/TCm9IJVvC5JaIiIgQGqoboQzY7SxAW8Z6xXtaAcTa3pwcJrdkO5jcEhERkckiIqwdAZFhTG6JiEgigwxFAV7SNBGRvWFyS1TD6G4ESUuzdiRkizxcnZEx7Ulx2sqxEBGZg8ktUQ2iVIo/JRa/EcTPz7oxERERWRKTWyIHZairnuI3gkREsAsfIjLdH39YOwIi45jcEjkgQzW0xZshcMhdKsvDIg2CPzsoTk/qyuF3SaLrC/e99x796qP7Ak1kS5jcEjmgsrrqISqPAAGu2XnSNJFO8b5wi/8aRGRrmNwSOTB21UNElqTrC5fIljG5JSIiokrRNXtiO36yBUxuiYiIyCy6drgjR4rPORQv2QImt0RERGSW4u1w09LEJPfXX0s3iWKNLlUnJrdERERkNl073JK1uMUpFMDWrYC/PxNdqnpMbomISCKDDGpfD2maqKKK1+IW99dfwLPPAn37is/ZdIGqGpNbIiKSeLg648+3e4nTVo6F7E9ZvSmUbLqQk8PklqoOk1siIiKqUuxCjKoTk1siIiKqVuw6jKoSk1siIpIUqDQI+neKOP1yNDzlHH6XLIddh1F1cLJ2AEREZDu0ggC3m/fgdvMetAKH3yXL0t10lpoqDg+en88hfMnyWHNLRERE1Ybtb6mqseaWiIiIiBwGa26JagjdDRxERESOjMktkY1TKsU2aRW5q1i3bvFE1tANHH5+VRcvERGRNTG5JbJhSqU4Rnt+ftl3FesSWt0oQPn54nxdElty1CB2vUNERI6MyS2RDcvJEZPVWbOA994rPapP8eQXEBPaHTtKj9/OGzjIFBpPV2uHQERkNia3RFakq3UFjNeohoWJf0t2fK5LfhMTxSSXtbJUWQpXF9yc3UectnIsRETmYHJLVEXKS1wN1bpu3QrUrg389ZeH3rpldXyuExEBtGtXRS+EiKgKcbQysjQmt0QWZKz9a8n2ssVrXf38xPX79gUAOdzceqFnz0cd6BdvN5uWJia57PiciOyZoS/tW7eWblZFZComt0QWUlb715ycR8mooYu1rtZVl7yeO6fG2LEuOHhQDZdi71C2m6XqUKDSIHDtMXH6hY4cfpeqTPEv7boKAfELPoflpcphcktkIWW1fz15smLb65JXHx8Bbm5igguU3XUX+62lqqAVBLhfuy1NE1Wl4l/aS/469euvQLduj5ZX9B4FIia3RBZW2favoaHAZ5/tQevWPeHiIi91EWe/tUTkiHSJbvFrXPH7C0r+MsaaXSoLk1siG+Tv/xCRkYBcXnoZ+60lIkemu8b9+qv+/QW6X8YA4029dFjTW3MxuSWqZoZGETMV298SkSMLDRVraoHSPcOUZCiJNXQPBGt6aw4mt0TVyNAFl00KiIhKK6sJli6RTUsru2ea4vdAABWr6SXHweSWqBpx0AUiooox1gSrZNJbsmcaHUM1veT4mNwSVZI5zQw46ALZMi27/yIbYagJVllJr65nGvYkQ0xuiSrBlGYGvOCSPVC4uuDGQrGzUQ6/S7bKUNJbkWYMuvX4i5ljY3JLVAkVaWbArruIiKqeKc0YeHOZY2NyS2QBxpoZsOsuIqLqUV4zhuLDl/M67LiY3BJVA3bdRfaiUKWBf2KqOD0yisPvkkMo6xqsu2eieKVD8a7FimPFhP1gcks1WlmdfFek82+lku1oyfFoBAGKS39J00SOqvg9EwoFsHWrOL9412LFsTmD/WBySw7BnJFoyurkGyg9f+tWwN/fcAfhbENLRGR/dPdMfPQR8M47QF/xPkqpazF//0frlmzOoFQC6ek+OHUKcCmWSbF21zYwuSW7ZyhJ3bpVnGfsImOok+9ffxWndfP9/MRv8cUvelu36m/brRsvZkRE9qpbN9Pui1AqgdatXZCf36PUsrJqdw01f6gK1XUcW+cwye3KlSvxr3/9C5mZmWjTpg0+/fRTdOzY0dph1RjVNYa3obZQaWmGk9HyLjLFa2kN9WigS1p1Fz3dSDjFE10mtkRE9s+U+yLEyg0Zpk5NxbBhreHiIgfwqHb311/1e88p+Uuf7nPJ0p+bZR2nKth6Eu0Qye3XX3+NadOmYfXq1ejUqRM++ugjxMXF4dKlSwgICLB2eA6vusbwLnmc4komo7/+aviO2LL6pTXWo0Hxix57PSAiIgCoX/8+IiMBuZjbGqwkKT4U8KxZwHvvPfoMsfTnpqHjVMVnVHUm0eZyiOR2+fLleOmllzBu3DgAwOrVq/Hzzz9j7dq1ePvtt60cnWGmfGOr6M1Nlry705T4jI3hbWg/ZbVVKk/xGtqSQyqWTEZ1y0ve8FVyH2UlsWVhrwdERPbP3JuBjW1nqMsxXVM3AAgLezRt6ucmUPHP+eLHqQrVlURXht0nt0VFRUhNTcXMmTOleU5OToiJiUFKSorBbQoLC1FYWCg9v3fvHgDg9u3bUKlUlY5JpVIhPz8ft27dglz3la6YGzeALl1c8PChDADg4SFg/XoN6tYtfWfyrVsyjBnjbHTdkusUZ2zfZanIMYu7fFkGwAXBwbqyk+PYMTXS01FqP8uXazBtmjMePmwH4OHfj4rz8BDQooUaISGG4n407eQEeHi4YORIw2VSfB/Ft7MF5Z0/1UFboAUKAK1ci1u2VkCwjTKyZZUpnztFamgLxeqkO7duwcnV7j8mSuH5Y5xKpUKBtgAqqCDXym3yGmCukp8NHh4CnJzU5X4OGNpOLr9T6hzy9BQfhtZ3ddUAcMGxY+q/167Y5+b69RoAKPdzXnwNLsjLU0vHuXvX8j2e6D7z69YVj3P3rqpU+VXVe+z+/fsAAKG8nlwEO/fnn38KAITDhw/rzX/zzTeFjh07Gtxm7ty5AgA++OCDDz744IMPPuzscePGDaO5oeN9Ja+AmTNnYtq0adJzrVaL27dvo27dupDJSn8rMlVubi5CQkJw48YNeHt7V3p/joblYxzLp3wsI+NYPsaxfIxj+ZSPZWRcVZWPIAi4f/8+6tWrZ3Q9u09u/fz84OzsjKysLL35WVlZCAoKMriNm5sb3Nzc9Ob5+vpaPDZvb2+e9EawfIxj+ZSPZWQcy8c4lo9xLJ/ysYyMq4ry8fHxKXcdJ4se0QpcXV0RFRWF5ORkaZ5Wq0VycjKio6OtGBkRERERVTe7r7kFgGnTpmHMmDFo3749OnbsiI8++ggPHjyQek8gIiIioprBIZLb5557Dn/99RfmzJmDzMxMtG3bFjt27EBgYKBV4nFzc8PcuXNLNX0gEcvHOJZP+VhGxrF8jGP5GMfyKR/LyDhrl49MEMrrT4GIiIiIyD7YfZtbIiIiIiIdJrdERERE5DCY3BIRERGRw2ByS0REREQOg8mtGRYtWoQuXbpAoVCUOfiDTCYr9di8ebPR/d6+fRsjRoyAt7c3fH19MX78eOTl5VXBK6ha5ZXPmTNnMHz4cISEhMDDwwMRERH4+OOPy91veHh4qTJdunRpFbyCqleRc0ipVGLAgAFQKBQICAjAm2++CbVabXBdHUc5h4rbt2+fwfeTTCbD8ePHy9yuR48epdZ/5ZVXqjHy6mXO+6OgoADx8fGoW7cuvLy8MGTIkFID4jiC69evY/z48WjQoAE8PDzQsGFDzJ07F0VFRUa3c+RzaOXKlQgPD4e7uzs6deqEY8eOGV1/y5YtaNasGdzd3dGqVSts3769miKtfkuWLEGHDh1Qq1YtBAQEYNCgQbh06ZLRbRISEkqdK+7u7tUUcfWaN29eqdfarFkzo9tU9/nD5NYMRUVF+Oc//4mJEycaXW/dunXIyMiQHoMGDTK6/ogRI3DhwgUkJSVh27ZtOHDgACZMmGDByKtHeeWTmpqKgIAAJCYm4sKFC3j33Xcxc+ZMfPbZZ+Xue8GCBXplOnnyZEuHXy3KKyONRoMBAwagqKgIhw8fxvr165GQkIA5c+YY3a+jnEPFdenSRe9/npGRgRdffBENGjRA+/btjW770ksv6W33wQcfVFPU1mHq+2Pq1Kn473//iy1btmD//v343//+h2effbaaoq0+v/32G7RaLf7973/jwoULWLFiBVavXo133nmn3G0d8Rz6+uuvMW3aNMydOxcnT55EmzZtEBcXh+zsbIPrHz58GMOHD8f48eNx6tQpDBo0CIMGDcL58+erOfLqsX//fsTHx+PIkSNISkqCSqVCbGwsHjx4YHQ7b29vvXPljz/+qKaIq1+LFi30XuvBgwfLXNcq549AZlu3bp3g4+NjcBkA4fvvv6/wvi5evCgAEI4fPy7N++WXXwSZTCb8+eeflYzUOoyVT0mvvvqq0LNnT6PrhIWFCStWrKh8YDakrDLavn274OTkJGRmZkrzVq1aJXh7ewuFhYUG9+WI55AhRUVFgr+/v7BgwQKj6z355JPC66+/Xj1B2QBT3x93794V5HK5sGXLFmleWlqaAEBISUmpgghtywcffCA0aNDA6DqOeg517NhRiI+Pl55rNBqhXr16wpIlSwyuP3ToUGHAgAF68zp16iS8/PLLVRqnrcjOzhYACPv37y9zHVM+7+zd3LlzhTZt2lR4fWucP6y5rULx8fHw8/NDx44dsXbtWghGuhROSUmBr6+vXk1UTEwMnJyccPTo0eoI16ru3buHOnXqlLve0qVLUbduXURGRuJf//pXuT/T26uUlBS0atVKbyCSuLg45Obm4sKFC2VuUxPOoZ9++gm3bt2q0AiEGzduhJ+fH1q2bImZM2ciPz+/GiK0HlPeH6mpqVCpVIiJiZHmNWvWDKGhoUhJSamOcK2qotccRzuHioqKkJqaqvd/d3JyQkxMTJn/95SUFL31AfF6VBPOE0A8VwCUe77k5eUhLCwMISEheOaZZ8q8VjuCK1euoF69enj88ccxYsQIKJXKMte1xvnjECOU2aIFCxagV69eUCgU2LVrF1599VXk5eXhtddeM7h+ZmYmAgIC9Oa5uLigTp06yMzMrI6Qrebw4cP4+uuv8fPPPxtd77XXXkO7du1Qp04dHD58GDNnzkRGRgaWL19eTZFWn8zMzFIj7Omel3U+1JRzaM2aNYiLi0P9+vWNrvf8888jLCwM9erVw9mzZzFjxgxcunQJW7duraZIq5ep74/MzEy4urqWavMdGBjoUOeLIVevXsWnn36KZcuWGV3PEc+hnJwcaDQag9eX3377zeA2ZV2PHP08AQCtVospU6bgiSeeQMuWLctcr2nTpli7di1at26Ne/fuYdmyZejSpQsuXLhQ7rXK3nTq1AkJCQlo2rQpMjIyMH/+fHTr1g3nz59HrVq1Sq1vlfOnyuqE7cyMGTMEAEYfaWlpetuY8jPE7Nmzhfr165e5fNGiRUKTJk1Kzff39xc+//xzk15LVaiq8jl37pzg5+cnLFy40OSY1qxZI7i4uAgFBQUmb1sVLFlGL730khAbG6s378GDBwIAYfv27QaPb+vnUEnmlNeNGzcEJycn4dtvvzX5eMnJyQIA4erVq5Z6CVXOnDLSKe/9sXHjRsHV1bXU/A4dOghvvfWWRV9HVTGnfG7evCk0bNhQGD9+vMnHs8dzqKQ///xTACAcPnxYb/6bb74pdOzY0eA2crlc2LRpk968lStXCgEBAVUWp6145ZVXhLCwMOHGjRsmbVdUVCQ0bNhQmDVrVhVFZjvu3LkjeHt7C//3f/9ncLk1zh/W3P7tjTfewNixY42u8/jjj5u9/06dOmHhwoUoLCw0ONZyUFBQqcb8arUat2/fRlBQkNnHtZSqKJ+LFy+id+/emDBhAmbNmmVyTJ06dYJarcb169fRtGlTk7e3NEuWUVBQUKm7l3V3sZd1Ptj6OVSSOeW1bt061K1bF08//bTJx+vUqRMAsdauYcOGJm9vDZU5p8p7fwQFBaGoqAh3797Vq73NysqyyfPFEFPL53//+x969uyJLl264D//+Y/Jx7PHc6gkPz8/ODs7l+oVw9j/PSgoyKT1HcWkSZOkG3NNrX2Vy+WIjIzE1atXqyg62+Hr64smTZqU+Vqtcf4wuf2bv78//P39q2z/p0+fRu3atQ0mtgAQHR2Nu3fvIjU1FVFRUQCAPXv2QKvVShdUa7J0+Vy4cAG9evXCmDFjsGjRIrP2cfr0aTg5OZX6Kd5aLFlG0dHRWLRoEbKzs6XXl5SUBG9vbzRv3rzMbWz5HCrJ1PISBAHr1q3D6NGjIZfLTT7e6dOnAQDBwcEmb2stlTmnynt/REVFQS6XIzk5GUOGDAEAXLp0CUqlEtHR0WbHXJ1MKZ8///wTPXv2RFRUFNatWwcnJ9NvObHHc6gkV1dXREVFITk5WerBR6vVIjk5GZMmTTK4TXR0NJKTkzFlyhRpXlJSkt2cJ6YSBAGTJ0/G999/j3379qFBgwYm70Oj0eDcuXPo379/FURoW/Ly8pCeno5Ro0YZXG6V86fK6oQd2B9//CGcOnVKmD9/vuDl5SWcOnVKOHXqlHD//n1BEAThp59+Er744gvh3LlzwpUrV4TPP/9cUCgUwpw5c6R9HD16VGjatKlw8+ZNaV7fvn2FyMhI4ejRo8LBgweFxo0bC8OHD6/211dZ5ZXPuXPnBH9/f2HkyJFCRkaG9MjOzpb2UbJ8Dh8+LKxYsUI4ffq0kJ6eLiQmJgr+/v7C6NGjrfIaK6u8MlKr1ULLli2F2NhY4fTp08KOHTsEf39/YebMmdI+HPkcMmT37t1l/gx/8+ZNoWnTpsLRo0cFQRCEq1evCgsWLBBOnDghXLt2Tfjxxx+Fxx9/XOjevXt1h10tKvL+KFlGgiD+5BoaGirs2bNHOHHihBAdHS1ER0db4yVUqZs3bwqNGjUSevfuLdy8eVPvulN8nZpyDm3evFlwc3MTEhIShIsXLwoTJkwQfH19pd5ZRo0aJbz99tvS+ocOHRJcXFyEZcuWCWlpacLcuXMFuVwunDt3zlovoUpNnDhR8PHxEfbt26d3ruTn50vrlCyj+fPnCzt37hTS09OF1NRUYdiwYYK7u7tw4cIFa7yEKvXGG28I+/btE65duyYcOnRIiImJEfz8/KTPcFs4f5jcmmHMmDEG23bt3btXEASx+6W2bdsKXl5egqenp9CmTRth9erVgkajkfaxd+9eAYBw7do1ad6tW7eE4cOHC15eXoK3t7cwbtw4KdmxJ+WVz9y5cw0uDwsLk/ZRsnxSU1OFTp06CT4+PoK7u7sQEREhLF682Gba25qqvDISBEG4fv260K9fP8HDw0Pw8/MT3njjDUGlUknLHfkcMmT48OFCly5dDC67du2aXvkplUqhe/fuQp06dQQ3NzehUaNGwptvvincu3evGiOuPhV5f5QsI0EQhIcPHwqvvvqqULt2bUGhUAiDBw/WS/gcxbp168psk6tT086hTz/9VAgNDRVcXV2Fjh07CkeOHJGWPfnkk8KYMWP01v/mm2+EJk2aCK6urkKLFi2En3/+uZojrj5lnSvr1q2T1ilZRlOmTJHKMzAwUOjfv79w8uTJ6g++Gjz33HNCcHCw4OrqKjz22GPCc889p9cO3RbOH5kgGOmfioiIiIjIjrCfWyIiIiJyGExuiYiIiMhhMLklIiIiIofB5JaIiIiIHAaTWyIiIiJyGExuiYiIiMhhMLklIiIiIofB5JaIqkyPHj30hly0pO7du2PTpk1Vsu/qFh4ejo8++shm9lPcsGHD8OGHH5a7nkwmg0wmg6+vrzRv3rx5aNu2rUXj0Zk9ezYmTJhg8f3KZDL88MMPZS6/ePEi6tevjwcPHujN79Gjh1QGumF6icg6mNwSkd356aefkJWVhWHDhknzwsPDpeTCw8MD4eHhGDp0KPbs2WPFSKtGQkKCXhKpc/z4cYsnfLNmzcKiRYtw7969ctddt24dLl++bNHjG5KZmYmPP/4Y7777rjRv48aNCAkJQe3atTFt2jS99a9fv44mTZogNzfX4P7mz5+PkSNHVujYzZs3R+fOnbF8+XK9+Vu3bsWxY8dMfCVEVBWY3BKR3fnkk08wbtw4ODnpX8IWLFiAjIwMXLp0CV9++SV8fX0RExODRYsWVXlMRUVFVX6M8vj7+0OhUFh0ny1btkTDhg2RmJhY7rq+vr4ICAiw6PEN+b//+z906dIFYWFhAICcnBy8+OKLWLZsGXbt2oXExERs27ZNWv/VV1/F0qVL4e3tbXB/P/74I55++ukKH3/cuHFYtWoV1Gq1NK9OnTrw9/c38xURkSUxuSWianPnzh2MHj0atWvXhkKhQL9+/XDlyhW9db744guEhIRAoVBg8ODBWL58uV4t5V9//YU9e/Zg4MCBpfZfq1YtBAUFITQ0FN27d8d//vMfzJ49G3PmzMGlS5ek9c6fP49+/frBy8sLgYGBGDVqFHJycqTl9+/fx4gRI+Dp6Yng4GCsWLGiVBOL8PBwLFy4EKNHj4a3t7dUY3rw4EF069YNHh4eCAkJwWuvvab3E3Z2djYGDhwIDw8PNGjQABs3biz1OpYvX45WrVrB09MTISEhePXVV5GXlwcA2LdvH8aNG4d79+5JNdXz5s2TYireLEGpVOKZZ56Bl5cXvL29MXToUGRlZUnLdc0GNmzYgPDwcPj4+GDYsGG4f/++XjwDBw7E5s2bS8VpKq1WiwULFqB+/fpwc3ND27ZtsWPHDr11Dh8+jLZt28Ld3R3t27fHDz/8UOqn/s2bN+v9/3///Xf4+PjgueeeQ4cOHdCzZ0+kpaUBAL766ivI5XI8++yzBmO6ceMGLly4gL59+0rzcnJyMHjwYCgUCjRu3Bg//fST3jZ9+vTB7du3sX///soWCRFVASa3RFRtxo4dixMnTuCnn35CSkoKBEFA//79oVKpAACHDh3CK6+8gtdffx2nT59Gnz59StW6Hjx4EAqFAhERERU65uuvvw5BEPDjjz8CAO7evYtevXohMjISJ06cwI4dO5CVlYWhQ4dK20ybNg2HDh3CTz/9hKSkJPz66684efJkqX0vW7YMbdq0walTpzB79mykp6ejb9++GDJkCM6ePYuvv/4aBw8exKRJk/TK4MaNG9i7dy++/fZbfP7558jOztbbr5OTEz755BNcuHAB69evx549e/DWW28BALp06YKPPvoI3t7eyMjIQEZGBqZPn14qNq1Wi2eeeUZKwpKSkvD777/jueee01svPT0dP/zwA7Zt24Zt27Zh//79WLp0qd46HTt2xLFjx1BYWFihMi/Lxx9/jA8//BDLli3D2bNnERcXh6efflr6gpObm4uBAweiVatW+P/27jWmqTOMA/i/BRoRRR14AabTUkSQOwPGyGQfKnEOI+wDhqFOLdTbxMtCWHaBTCUGoyULBC8LcdM4RRMRjVbNjJbLEiBQitJSvBFQQzJFqAhOCs8+EE52VtyAgTr2/BITz3ve23kPHx7Oec5LTU0Ndu3ahfT0dFEfbW1tMBqNePfdd4UyLy8vdHV1Qa/Xo62tDVVVVQgICMCTJ0/w7bffIi8v76VzOnfuHD788EPRU93vvvsOCQkJqKurw9KlS5GUlIS2tjbhvEwmQ1BQEEpLS//VejDGxggxxtgYiY6Opq1btxIRUWNjIwGg8vJy4fyjR4/I0dGRTp06RUREK1asoI8//ljUR1JSEk2ZMkU4zsnJIblcbjPWO++8Qzk5OYPOY+bMmbRx40YiItq1axfFxMSIzre0tBAAMpvNZLFYyMHBgU6fPi2cb29vp4kTJwrXMjBeXFycqB+VSkVqtVpUVlpaSlKplLq7u8lsNhMAqqysFM6bTCYC8NK5ExGdPn2aXFxchOMjR46I1mSwNbhy5QrZ2dlRc3OzcL6+vl40fmZmJk2cOJEsFotQJy0tjSIiIkT9GgwGAkBNTU0vnSMAKioqEpVlZmZSYGCgcOzu7k5ZWVmiOmFhYbRp0yYiIjpw4AC5uLhQd3e3cP6HH34gAKTX64mISK/XEwDRdRERnTlzhvz8/MjT05MyMzOJiGjdunWUk5NDOp2OgoKCaOHChaL7SkS0ePFiysvLE13HN998Ixx3dnYSANJqtaJ28fHxtGbNGlHZvXv3RHNljL0e9q8npGaM/d+YTCbY29sjIiJCKHNxcYG3t7fwCtlsNiM+Pl7ULjw8XJQ/2d3djQkTJgxrbCKCRCIBABgMBly7dg2TJk2yqXfnzh10d3ejp6cH4eHhQvmUKVPg7e1tU//PTw8H+q6rqxOlGhAR+vr6cO/ePTQ2NsLe3h6hoaHC+QULFth8HPbLL79gz549aGhogMVigdVqxfPnz9HV1TXknFqTyYTZs2dj9uzZQpmvry+mTp0Kk8mEsLAwAP2pDJMnTxbquLm52TxJdnR0BAB0dXUNaezBWCwWPHz4EFFRUaLyqKgoGAwGAP33PyAgQHR//3wfgP77D8DmZyA+Pl70s6PT6VBXV4fc3FwoFAqcOHECs2bNQnh4OBYtWoQZM2bAYrFAp9OhoKBA1FdAQIDwfycnJzg7Ow+6Jv9mPRhjY4eDW8bYf4qrqyuePHky5PqPHz/Gb7/9hnnz5gEAOjs7sWzZMmRnZ9vUdXNzw+3bt4fct5OTk+i4s7MT69evR2pqqk3dOXPmDGkngaamJsTGxmLjxo3IysrCW2+9hbKyMqhUKrx48WLUPxhzcHAQHUskEvT19YnKBl7JvwkfTLm6ugLoz99+2Xx+//13bNq0CceOHcPt27dhtVoRHR0NAJg/fz4qKiqwbNkyaLVa+Pr6in4BAIa+Jp6enqN1WYyxUcQ5t4yxV8LHxwdWqxUVFRVC2ePHj2E2m+Hr6wsA8Pb2RlVVlajdX4+Dg4PR2to65AD3+++/h1QqRVxcHAAgJCQE9fX1mDt3LhQKheifk5MT5HI5HBwcRON2dHQMKTANCQmB0Wi06VehUEAmk2HBggWwWq2orq4W2pjNZrS3twvH1dXV6Ovrw/79+/Hee+9h/vz5ePjwoWgcmUyG3t7ev52Lj48PWlpa0NLSIpQZjUa0t7cL6z1UN2/exNtvvy0EliPh7OwMd3d3lJeXi8rLy8tF9//GjRui3N6/3n9PT084OzvDaDS+dKzdu3djyZIlCAkJQW9vr2hXg56eHmHtiouLsXz58hFdz82bNxEcHDyitoyxscXBLWPslfDy8sLy5cuRkpKCsrIyGAwGrFy5Eh4eHkKAsWXLFly8eBEajQa3bt3CoUOHoNVqhZQCoD+4dXV1tQmSgP5dDlpbW9HS0oKSkhKo1Wrs3r0bWVlZUCgUAIDNmzejra0NiYmJqKqqwp07d3D58mWsXbsWvb29mDx5Mj777DOkpaXh2rVrqK+vh0qlglQqFc1jMOnp6fj111/x+eefo7a2Frdu3UJxcbHwQZm3tzeWLFmC9evXo6KiAtXV1UhOThZe+wOAQqFAT08PcnNzcffuXRw7dgwHDx4UjTN37lx0dnbi6tWrePTo0aCvx5VKJfz9/ZGUlISamhpUVlZi9erViI6Otkmn+CelpaWIiYkZVpvBpKWlITs7G4WFhTCbzfjyyy9RW1uLrVu3AgA+/fRT9PX1Qa1Ww2Qy4fLly9i3bx8ACGsvlUqhVCpRVlY26BhGoxGFhYXYuXMngP60D6lUioKCAly4cAENDQ0ICwuD1WqFVqsd1hZgA5qamvDgwQMolcqRLANjbIxxcMsYe2WOHDmC0NBQxMbGIjIyEkSEixcvCq+Bo6KicPDgQWg0GgQGBuLSpUvYvn27KL/Szs4Oa9euHXQLrYyMDLi5uUGhUGDVqlXo6OjA1atXRV/cDzw97O3tRUxMDPz9/bFt2zZMnTpV2DdXo9EgMjISsbGxUCqViIqKgo+Pzz/m+gYEBECn06GxsREffPABgoODkZGRAXd3d9EauLu7Izo6Gp988gnUarVob9jAwEBoNBpkZ2fDz88Px48fx549e0TjvP/++9iwYQNWrFiB6dOnY+/evTZzkUgkKC4uxrRp07Bo0SIolUrI5XIUFhb+7TX81fPnz3H27FmkpKQMq91gUlNTsWPHDnzxxRfw9/fHpUuXcO7cOXh5eQHof7p7/vx51NbWIigoCF9//TUyMjIAiHNsk5OTcfLkSZtUASKCWq2GRqMRUkYcHR3x448/YufOnVCpVMjLy4OHhwd0Oh0mTZqEkJCQYV/HiRMnEBMTI+yzyxh7s0iIiF73JBhj7GVSUlLQ0NAg2naptbUVCxcuRE1NzSsJMJ49ewYPDw/s378fKpVqzMd7kxw4cABFRUW4cuXK39aTSCQoKioS0j9Gy/Hjx4V9fQeecBMRIiIisH37diQmJo6o39TUVFitVuTn5w+r3YsXL+Dl5YWff/7Z5uO4pqYmzJs3D3q9fsz+7DBj7J/xB2WMsTfKvn37sHjxYjg5OUGr1eKnn36yCUBmzZqFgoICNDc3j0lwq9fr0dDQgPDwcHR0dAivuEean/lf5uDggNzc3CHVTUxMhIuLC+7fvz/i8Y4ePQq5XA4PDw8YDAakp6cjISFBlLohkUhw+PBh3LhxY8Tj+Pn5ITIyctjtmpub8dVXX9kEth999BFKSkpGPB/G2OjhJ7eMsTdKQkICrl+/jqdPn0Iul2PLli3YsGHDK52DXq9HcnIyzGYzZDIZQkNDhb8axgY3sMuEnZ2dsDPFSOzduxf5+flobW2Fm5sb4uLikJWVNeq7RIy2Bw8eCNuUzZkzBzKZ7DXPiLH/Lw5uGWOMMcbYuMEflDHGGGOMsXGDg1vGGGOMMTZucHDLGGOMMcbGDQ5uGWOMMcbYuMHBLWOMMcYYGzc4uGWMMcYYY+MGB7eMMcYYY2zc4OCWMcYYY4yNGxzcMsYYY4yxceMPU1+7KyvQZ0sAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "lnDeg = np.log10(results)\n", + "percentile_2p5 = np.percentile(lnDeg, 2.5)\n", + "percentile_97p5 = np.percentile(lnDeg, 97.5)\n", + "bin_edges = np.arange(lnDeg.min(), lnDeg.max() + 0.1, 0.1)\n", + "\n", + "plt.figure(figsize=(8,6))\n", + "plt.hist(lnDeg, bins=bin_edges, edgecolor='blue', histtype='step', linewidth=1, label = 'Degradation in Miami')\n", + "\n", + "plt.axvline(percentile_2p5, color='green', label='2.5th percentile', linewidth=2.0)\n", + "plt.axvline(percentile_97p5, color='purple', label='97.5th percentile', linewidth=2.0)\n", + "plt.axvline(np.mean(lnDeg), color = 'cyan', label = 'mean' )\n", + "plt.axvline(np.median(lnDeg), linestyle='--', label = 'median')\n", + "plt.xlabel('log(Degredation) [log(%/h)]')\n", + "plt.ylabel(f'Counts (out of {n})')\n", + "\n", + "plt.legend()\n", + "plt.grid(True)\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.10.9" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/pvdeg_tutorials/tutorials/Monte Carlo - Standoff.ipynb b/pvdeg_tutorials/tutorials/Monte Carlo - Standoff.ipynb new file mode 100644 index 00000000..d5bc5fee --- /dev/null +++ b/pvdeg_tutorials/tutorials/Monte Carlo - Standoff.ipynb @@ -0,0 +1,359 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Monte Carlo - Standoff Calculation\n", + "\n", + "Author: Tobin Ford | tobin.ford@nrel.gov\n", + "\n", + "2023\n", + "***\n", + "See Monte Carlo - Arrhenius Degredation for a more in depth guide. Steps will be shortened for brevity.\n", + "\n", + "**Objectives**\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# if running on google colab, uncomment the next line and execute this cell to install the dependencies and prevent \"ModuleNotFoundError\" in later cells:\n", + "# !pip install pvdeg==0.2.0" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import pvlib\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pvdeg\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Working on a Windows 10\n", + "Python version 3.10.9 | packaged by Anaconda, Inc. | (main, Mar 8 2023, 10:42:25) [MSC v.1916 64 bit (AMD64)]\n", + "Pandas version 2.1.2\n", + "Pvlib version 0.10.2\n", + "Pvdeg version 0.2.0+12.g9849aa2.dirty\n" + ] + } + ], + "source": [ + "# This information helps with debugging and getting support :)\n", + "import sys, platform\n", + "print(\"Working on a \", platform.system(), platform.release())\n", + "print(\"Python version \", sys.version)\n", + "print(\"Pandas version \", pd.__version__)\n", + "print(\"Pvlib version \", pvlib.__version__)\n", + "print(\"Pvdeg version \", pvdeg.__version__)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Simple Standoff Calculation\n", + "\n", + "This is copied from another tutorial called `4 - Standards.ipynb`, please visit this page for a more in depth explanation of the process for a single standoff calculation.\n", + "\n", + "
\n", + "Please use your own API key: The block below makes an NSRDB API to get weather and meta data. This tutorial will work with the DEMO Key provided, but it will take you less than 3 minutes to obtain your own at https://developer.nrel.gov/signup/ so register now.) \n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " \r" + ] + } + ], + "source": [ + "weather_db = 'PSM3'\n", + "weather_id = (40.633365593159226, -73.9945801019899) # Manhattan, NYC\n", + "weather_arg = {'api_key': 'DEMO_KEY',\n", + " 'email': 'user@mail.com', \n", + " 'names': 'tmy',\n", + " 'attributes': [],\n", + " 'map_variables': True}\n", + "\n", + "WEATHER, META = pvdeg.weather.get(weather_db, weather_id, **weather_arg)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " x T98_0 T98_inf T98\n", + "0 0.0 67.0331 45.660409 70.0\n", + " x T98_0 T98_inf T98\n", + "0 0.0 69.139133 46.496454 70.0\n" + ] + } + ], + "source": [ + "# simple standoff calculation\n", + "height1 = pvdeg.standards.standoff(\n", + " weather_df = WEATHER,\n", + " meta = META\n", + " )\n", + "\n", + "# more arguments standoff calculation\n", + "height2 = pvdeg.standards.standoff(weather_df=WEATHER, meta=META,\n", + " tilt=None,\n", + " azimuth=180,\n", + " sky_model='isotropic',\n", + " temp_model='sapm',\n", + " x_0=6.1,\n", + " wind_factor = 0.33 # default\n", + " )\n", + "\n", + "print(height1)\n", + "print(height2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Defining Correlation Coefficients, Mean and Standard Deviation For Monte Carlo Simulation\n", + "\n", + "We will leave the list of correlations blank because our variables are not correlated. For a correlated use case visit the `Monte Carlo - Arrhenius.ipynb` tutorial.\n", + "\n", + "Mean and standard deviation must always be populated if being used to create a dataset. However, you can feed your own correlated or uncorrelated data into the simulate function but column names must be consistent." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# These numbers may not make sense in the context of the problem but work for demonstraiting the process\n", + "stats = {\n", + " 'X_0' : {'mean' : 5, 'stdev' : 3},\n", + " 'wind_factor' : {'mean' : 0.33, 'stdev' : 0.5}\n", + "}\n", + "\n", + "corr_coeff = []\n", + "\n", + "samples = pvdeg.montecarlo.generateCorrelatedSamples(corr_coeff, stats, 500)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " X_0 wind_factor\n", + "0 0.898718 0.459997\n", + "1 6.956492 0.335493\n", + "2 7.103368 0.598136\n", + "3 7.305353 -0.527580\n", + "4 2.698503 -0.051998\n", + ".. ... ...\n", + "495 -0.534005 0.620055\n", + "496 4.696580 0.518921\n", + "497 1.657083 0.041645\n", + "498 3.767393 0.368091\n", + "499 3.309681 0.483557\n", + "\n", + "[500 rows x 2 columns]\n" + ] + } + ], + "source": [ + "print(samples)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Standoff Monte Carlo Inputs\n", + "\n", + "When using the pvdeg.montecarlo.simulate() function on a target function all of the target function's required arguments must still be given. Our non-changing arguments will be stored in a dictionary. The randomized monte carlo input data will also be passed to the target function via the simulate function. All required target function arguments should be contained between the column names of the randomized input data and fixed argument dictionary, " + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "import pdb\n", + "\n", + "# defining arguments to pass to the target function, standoff() in this case\n", + "function_kwargs = {\n", + " 'weather_df' : WEATHER,\n", + " 'meta' : META,\n", + " 'azimuth' : 180,\n", + " 'tilt' : 0,\n", + " 'temp_model' : 'sapm',\n", + " 'sky_model' : 'isotropic',\n", + " 'conf_0' : \"insulated_back_glass_polymer\", # or is it conf_inf\n", + " 'conf_inf' : \"open_rack_glass_polymer\",\n", + " 'T98' : 70,\n", + "} \n", + "\n", + "# notice how we left off parts we want to use in the monte carlo simulation because they are already contained in the dataframe\n", + "\n", + "results = pvdeg.montecarlo.simulate(\n", + " func=pvdeg.standards.standoff, \n", + " correlated_samples=samples, # in this case correlated_samples is a misnomer, they are not required to be correlated\n", + " **function_kwargs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Dealing With Series \n", + "Notice how our results are contained in a pandas series instead of a dataframe.\n", + "\n", + "This means we have to do an extra step to view our results. Run the block below to confirm that our results are indeed contained in a series. And convert them into a simpler dataframe." + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n" + ] + } + ], + "source": [ + "print(type(results))\n", + "\n", + "# Convert from pandas Series to pandas DataFrame\n", + "results_df = pd.concat(results.tolist()).reset_index(drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " x T98_0 T98_inf T98\n", + "0 0.000000 65.508742 44.608501 70.0\n", + "1 0.000000 66.973564 45.611088 70.0\n", + "2 0.000000 63.445939 43.421469 70.0\n", + "3 1.337077 73.952652 50.320063 70.0\n", + "4 0.120627 70.988196 48.383881 70.0\n", + ".. ... ... ... ...\n", + "495 0.160900 63.027109 43.196684 70.0\n", + "496 0.000000 64.583898 44.101643 70.0\n", + "497 0.025518 70.345243 47.753195 70.0\n", + "498 0.000000 66.674085 45.392219 70.0\n", + "499 0.000000 65.243865 44.415142 70.0\n", + "\n", + "[500 rows x 4 columns]\n" + ] + } + ], + "source": [ + "print(results_df)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Viewing Our Data\n", + "Let's plot the results using a histogram" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArcAAAINCAYAAAAkzFdkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAABkj0lEQVR4nO3deXxTVf7/8XeSphvdKNAFKfu+IwgWEJClLIogzCguCIqoDOhAx41REFAHlxFQBvGro+DGzx0cUYGyFGSXTUQZ1LIqLWUvbWmbJvf3B5KhdiENaVLC6/l45EFy78m9n/tJGj49PTnHZBiGIQAAAMAPmH0dAAAAAOApFLcAAADwGxS3AAAA8BsUtwAAAPAbFLcAAADwGxS3AAAA8BsUtwAAAPAbFLcAAADwGwG+DqAycDgcOnz4sMLDw2UymXwdDgAAAP7AMAydOXNGNWvWlNlcev8sxa2kw4cPKyEhwddhAAAA4CIOHTqkWrVqlbqf4lZSeHi4pHPJioiIqNBz2ewOfbh5v3bv/q+euL2XQoODKvR85ZaTI9Wsee7+4cNSlSq+jcfDbDabli1bpqSkJFmtVl+Hc8Uh/75D7n2L/PsOufctT+Y/KytLCQkJzrqtNBS3knMoQkRERIUXt7kFhXph1W+SwvVMlTBFVAmp0POVm8Xyv/sREX5Z3IaGhioiIoIPOR8g/75D7n2L/PsOufetisj/xYaQ8oUyAAAA+A2KWwAAAPgNilsAAAD4DcbcAgDgx+x2u2w2m6/D8BmbzaaAgADl5eXJbrf7OpwrTnnyb7FYFBAQcMnTslLcAgDgp7Kzs/Xrr7/KMAxfh+IzhmEoLi5Ohw4dYi57Hyhv/kNDQxUfH6/AwEC3z0lxCwCAH7Lb7fr1118VGhqqGjVqXLGFncPhUHZ2tsLCwsqc+B8Vw9X8G4ahgoICHT16VPv27VOjRo3cfr0obr0s0GLW63e205YtWxRo4YcMAFAxbDabDMNQjRo1FBJSyaad9CKHw6GCggIFBwdT3PpAefIfEhIiq9WqAwcOOJ/jDopbLwuwmHV9kxo6m2YogOIWAFDBrtQeW1yePPELCNUVAAAA/AY9t15mszv06bbftDPTpD52h1gsBQDgTQcPSseOee981atLtWt773wAxa2X2ewOPb7wB0kWPWZ3+DocAMAV5OBBqVkzKTfXe+cMDZV27758C9z9+/erXr162r59u9q2bevScwzD0P33369PPvlEJ0+e1Pbt29WmTZti21w9XklGjhypU6dOadGiRW4fw19R3AIAcIU4duxcYfvee+eK3Iq2e7d0553nzutqcXv06FFNnjxZX375pY4cOaKqVauqTZs2mjx5srp06SLp3DjihQsXavDgwRUX/CVYsmSJ5s+fr9TUVNWvX1/Vq1cvcdsfpaam6vrrr5d07hrDw8NVv3599enTRxMmTFB8fLyz7csvv+zyFG9XWiFMcQsAwBWmWTPp6qt9HUXJhg4dqoKCAr399tuqX7++jhw5ohUrVuj48eO+Ds1laWlpio+PV+fOncvcVpo9e/YoIiJCWVlZ2rZtm1544QW9+eabSk1NVatWrSRJkZGRFRb/5Y4vlAEAgErh1KlT+uabb/T888/r+uuvV506ddSxY0dNnDhRN910kySpbt26kqSbb75ZJpPJ+TgtLU2DBg1SbGyswsLCdM0112j58uVFjl+3bl394x//0D333KPw8HDVrl1br7/+epE2mzdvVrt27RQcHKwOHTpo+/btxeJcvXq1OnbsqKCgIMXHx+vxxx9XYWGhpHO9pA8++KAOHjzojK+kbWWJiYlRXFycGjdurGHDhmndunWqUaOGxowZ42wzcuTIIj3Xn3zyiVq1aqWQkBBVq1ZNvXv3Vk5OjqZMmaK3335bn3/+uUwmk0wmk1JTUyVJjz32mBo3bqzQ0FDVr19fkyZNKrKa3ZQpU9S2bVu9++67qlu3riIjIzVs2DCdOXPG2cbhcOiFF15Qw4YNFRQUpNq1a+vZZ5917j906JDuvvtuRUdHKzo6WoMGDdL+/fvLvP5LRXELAAAqhbCwMIWFhWnRokXKz88vsc23334rSZo3b57S09Odj7OzszVgwACtWLFC27dvV79+/TRw4EAdPHiwyPNfeuklZ9H6l7/8RWPGjNGePXucx7jxxhvVvHlzbd26VVOmTNHDDz9c5Pm//fabBgwYoGuuuUbfffed5s6dqzfffFPPPPOMpHPDBaZNm6ZatWo54ytpW3mEhITogQce0Lp165SZmVlsf3p6um677Tbdc8892r17t1JTUzVkyBAZhqGHH35Yt9xyi/r166f09HSlp6c7e4/Dw8M1f/58/fjjj3r55Zf1xhtvaObMmUWOnZaWpkWLFmnx4sVavHixVq9ereeee865f+LEiXruuec0adIk/fjjj1qwYIFiY2MlnZtruX///goLC9Pq1au1bt06hYWFqV+/fiooKChXDsqDYQkAAKBSCAgI0Pz58zV69Gi99tpruvrqq9W9e3cNGzZMrVu3liTVqFFDkhQVFaW4uDjnc9u0aaM2bdo4Hz/99NNauHChvvjiCw0fPty5fcCAAfrLX/4i6VzP5cyZM7Vq1So1adJECxYskMPh0Jtvvqng4GC1aNFCv/76a5Ee01dffVUJCQn617/+JZPJpKZNm+rw4cN67LHHNHnyZEVGRio8PFwWi6VIfCVtK4+mTZtKOvcFt5iYmCL70tPTVVhYqCFDhqhOnTqS5By+IJ0rjvPz84ud+8knn3Ter1u3rh5++GF98MEHevTRR53bHQ6H5s+fr/DwcEnS8OHDtWLFCj377LM6c+aMXn75Zf3rX//SiBEjJEkNGjRQ165dJUkffvihHA6HXnnlFUVGRspsNmvevHmKiopSamqqkpKS3MrFxdBzCwAAKo2hQ4fq8OHD+s9//qN+/fopNTVVV199tebPn1/m87Kzs/Xwww+rWbNmioqKUlhYmHbv3l2s5/Z8kSyd+9JWXFycszd09+7dat26dZGVsRITE4s8f/fu3UpMTCyyOEaXLl2UnZ2tX3/91d3LvqjzXx4raVGONm3aqFevXmrVqpX+/Oc/64033tDJkycveswPP/xQXbp0UVxcnMLCwvTkk08Wy1fdunWdha0kxcfHF8lXfn6+evXqVeLxv/vuO/3yyy9KSEhQRESEwsLCFB0drby8PKWlpbl87eVFcetlgRazXrm1tUY2trP8LgAAJQgODlafPn00adIkrV+/XiNHjtRTTz1V5nMefvhhLVy4UP/4xz/0zTffaMeOHWrVqlWxP39b/zDBvMlkksNR+afm3L17tySVOF7XYrEoJSVFX3/9tZo3b67Zs2erSZMm2rdvX6nH27Bhg+644w4NGDBAixcv1vbt2/XEE0+UK18XW9Y5Oztb7du315o1a7Rt2zbt2LFDO3bs0E8//aTbb7/dlct2C9WVlwVYzOrfMk61HMH6fqdZ27ap1NsffnkCAOCK1Lx5c+Xk5DgfW61W2e32Im3WrVunkSNH6uabb1arVq0UFxdX7i8uNWvWTDt37lReXp5z28aNG4u12bBhQ5FpuNatW6fw8HDVqlWrXOdz1dmzZ/X666+rW7duzmEZf2QymdSlSxdNnTpV27dvV2BgoBYuXChJCgwMLJav9evXq06dOnriiSfUoUMHNWrUSAcOHChXXI0aNVJISIhWrFhR4v6rr75aP//8s6pXr66GDRsWuVXkbA+MufWBgwelceN6Kj+/7PRf7hNfAwAqp987ASvdeY4fP64///nPuueee9S6dWuFh4dry5YteuGFFzRo0CBnu7p162rFihXq0qWLgoKCVLVqVTVq1EifffaZBg4cKJPJpEmTJpW7R/b222/XE088odGjR2vixInav3+//vnPfxZp85e//EWzZs3Sgw8+qHHjxmnPnj166qmnlJycLLPZM32GmZmZysvL05kzZ7R161a98MILOnbsmD777LMS22/atEkrVqxQUlKSYmJitGnTJh09elTNfp/MuG7dulq6dKn27NmjatWqKTIyUo0aNdLBgwf1wQcf6JprrtGXX37pLIZdFRwcrMcee0yPPvqoAgMD1aVLFx09elQ//PCDRo0apTvuuEMvvvii7rjjDj3zzDOqXbu2Dhw4oM8++0yPPvpohf0yQHHrZYV2hxZ/lyFL3QC9+Xg1tW0dWGI7dya+BgCgLNWrn+s4ufNO750zNPTceV0RFhamTp06aebMmUpLS5PNZlNCQoJGjx6tv//97852L730kpKTk/XGG2/oqquu0v79+zVjxgzdc8896ty5s6pXr67HHntMWVlZ5Yo1LCxMX3zxhR544AG1a9dOzZs31/PPP6+hQ4c621x11VX66quv9Mgjj6hNmzaKjo7WqFGjinw561I1adJEJpNJYWFhql+/vpKSkpScnFzql9EiIiK0Zs0azZo1S1lZWapTp45eeukl9e/fX5I0evRopaamqkOHDsrOztaqVat00003acKECRo3bpzy8/N1ww03aNKkSZoyZUq5Yp00aZICAgI0efJkHT58WPHx8XrggQckSaGhoUpNTdXf/vY3/elPf9KZM2d01VVXqVevXoqIiLikHJXFZLi6vIUfy8rKUmRkpE6fPl2hyZak3IJCNZ+8VJL0zk091a1zyeNVtm2T2reXtm718kTbOTlSWNi5+9nZUpUqXjx5xbPZbPrqq680YMCAYuOIUPHIv++Qe9/yRf7z8vK0b98+1atXr8gXpA4ePNdx4i3Vq/u2k8bhcCgrK0sREREe61mF68qb/9Let5Lr9Ro9twAAXEFq1+YvgvBv/AoDAAAAv0FxCwAAAL9BcQsAAAC/QXELAAAAv0FxCwAAAL9BcetlVotZf2nfQse+bK0ApiQBAADwKKorL7NazOpR5yrl7EqguAUAAPAwqisAAAD4DRZx8LJCu0Pb0o8qpL5FdkeUr8MBAADwK/TcelmB3aHnNmxXzJ+3yOZw+DocAAAAv0JxCwAAKo0ePXrowQcf1Pjx41W1alXFxsbqjTfeUE5Oju6++26Fh4erYcOG+vrrr53P2bVrl/r376+wsDDFxsZq+PDhOnbsmHP/8uXL1a1bN0VFRalatWq68cYblZaW5ty/f/9+mUwmffbZZ7r++usVGhqqNm3aaMOGDV69dngGxS0AAFcCw5BycnxzM4xyhfr222+revXq2rx5sx588EGNGTNGf/7zn9W5c2dt27ZNSUlJGj58uHJzc3Xq1Cn17NlT7dq105YtW7RkyRIdOXJEt9xyi/N4ubm5Gj9+vLZs2aIVK1bIbDbr5ptvluMPf0F94okn9PDDD2vHjh1q3LixbrvtNhUWFnok/fAextwCAHAlyM2VwsJ8c+7sbKlKFZebt2nTRk8++aQkaeLEiXruuedUvXp1jR49WpI0efJkzZ07Vzt37tTy5cvVrl07/eMf/3A+/6233lJCQoJ++uknNWzYUDfddJMiIiJk/n2Worfeeks1atTQjz/+qJYtWzqf9/DDD+uGG26QJE2dOlUtWrTQL7/8oqZNm15yCuA99NwCAIBKpXXr1s77FotF1apVU6tWrZzbYmNjJUmZmZn67rvvtGrVKoWFhTlv54vR80MP0tLSdPvtt6t+/fqKiIhQ3bp1JUkHDx4s9bzx8fHOc+DyQs8tAABXgtDQcz2ovjp3OVit1iKPTSZTkW0mk0mS5HA4lJ2drYEDB+r5558vdpzzBeptt92munXr6o033lDNmjXlcDjUsmVLFRQUlHreC8+BywvFLQAAVwKTqVxDAy4XV199tT799FPVrVtXAQHFy5qjR4/q559/1htvvKHu3btLktauXevtMOFFDEvwMqvFrHvaNNXxZS1YoQwAgEs0duxYnThxQrfddpu+/fZbpaWlaenSpbr77rtlt9tVtWpVRUdH64033tAvv/yilStXKjk52ddhowJRXXmZ1WJWvwa1lb29LsUtAACXqGbNmlq3bp3sdruSkpLUqlUrjR8/XlFRUTKbzTKbzXrzzTe1bds2tWzZUhMmTNCLL77o67BRgRiWAAAAKo3U1NRi2/bv319sm3HB9GKNGjXSZ599VuLxDMNQjx49tGvXLudsCX98ft26dYs8lqSoqKhi23B5oLj1MrvD0A9HTygoIUAOI9zX4QAAAPgV/i7uZfmFdk39Zovibt+oArvd1+EAAAD4FYpbAAAA+A2KWwAAAPgNilsAAAD4DYpbAAAA+A2fFrdz585V69atFRERoYiICCUmJurrr7927s/Ly9PYsWNVrVo1hYWFaejQoTpy5EiRYxw8eFA33HCDQkNDFRMTo0ceeUSFhYXevhQAAABUAj4tbmvVqqXnnntOW7du1ZYtW9SzZ08NGjRIP/zwgyRpwoQJ+uKLL/Txxx9r9erVOnz4sIYMGeJ8vt1u1w033KCCggKtX79eb7/9tubPn6/Jkyf76pIAAADgQz4tbgcOHKgBAwaoUaNGaty4sZ599lmFhYVp48aNOn36tN58803NmDFDPXv2VPv27TVv3jytX79eGzdulCQtW7ZMP/74o9577z21bdtW/fv319NPP605c+aooKDAl5dWqgCzWXe2bKSTq5qyQhkAAICHVZrqym6364MPPlBOTo4SExO1detW2Ww29e7d29mmadOmql27tjZs2CBJ2rBhg1q1aqXY2Fhnm759+yorK8vZ+1vZBAaYdVPjesra3IDiFgAAH+nRo4fGjx/vfFy3bl3NmjXLZ/HAc3y+Qtn333+vxMRE5eXlKSwsTAsXLlTz5s21Y8cOBQYGKioqqkj72NhYZWRkSJIyMjKKFLbn95/fV5r8/Hzl5+c7H2dlZUmSbDabbDabJy6rTDZboSSrbLZClXa6c8OGrSostJXapoKCk9V51ybvnrzinX99vfE6ozjy7zvk3rd8kX+bzSbDMORwOORwOLx23srm/BK653Pxx33nt23atElVqlS5onNVEcrKf0kcDocMw5DNZpPFYimyz9WfH58Xt02aNNGOHTt0+vRpffLJJxoxYoRWr15doeecPn26pk6dWmz7smXLFBoaWqHndhjSxj1hCoxrrw0bt+vYsawS26WlRUrqobVr1yk9/XSFxnQhS16ebvz9/tKlS2UPDvbaub0pJSXF1yFc0ci/75B73/Jm/gMCAhQXF6fs7OxKO1TPm86cOVPkcWFhoQoKCpwdXEFBQSosLHQ+hmf9Mf+lKSgo0NmzZ7VmzZpiEwTk5ua6dAyfF7eBgYFq2LChJKl9+/b69ttv9fLLL+vWW29VQUGBTp06VaT39siRI4qLi5MkxcXFafPmzUWOd342hfNtSjJx4kQlJyc7H2dlZSkhIUFJSUmKiIjw1KWVKLegUBM2rlT8iHVq16G7unUOKrHd9u3n/u3atYvatavQkIrKyXHe7du3r1SlihdPXvFsNptSUlLUp08fWa3Wiz8BHkX+fYfc+5Yv8p+Xl6dDhw4pLCxMwZdRR0XPnj3VsmVLWSwWvfPOOwoMDNS0adN0++2368EHH9Snn36q2NhYvfzyy+rfv78kadeuXXr00Ue1du1aValSRX369NGMGTNUvXp1GYahjIwMPf7441q4cKHCw8P1t7/9TQEBAQoMDHT+v1+/fn399a9/1V//+ldJ0syZMzV//nzt3btX0dHRuvHGG/X8888rLCxMkjR//nwlJyfr//2//6fk5GQdOnRIXbp00VtvvaX4+HjfJK8SMgxDZ86cUXh4uEwm00Xb5+XlKSQkRN26dSv2vnX1Fw+fF7d/5HA4lJ+fr/bt28tqtWrFihUaOnSoJGnPnj06ePCgEhMTJUmJiYl69tlnlZmZqZiYGEnnfiuOiIhQ8+bNSz1HUFCQgoKKF5VWq7XCP3Ssxv9eWKvVUur5AgLO/2uVV/8fuuBkVqtV3j2593jjtUbpyL/vkHvf8mb+7Xa7TCaTzGazzBd8xyO3oPTpMs0mk4KtFo+2DQ0sf6nxzjvv6NFHH9XmzZv14YcfauzYsfr88891880364knntDMmTM1YsQIHTx4UAUFBerdu7fuvfdezZo1S2fPntVjjz2mYcOGaeXKlXI4HJo8ebLWrFmjzz//XDExMfr73/+ubdu2qW3btkVycz5fkmSxWPTKK6+oXr162rt3r/7yl7/o8ccf16uvvnru+s1m5ebmasaMGXr33XdlNpt155136tFHH9X7779f7mv2V+eHIlyY27KYzWaZTKYSf1Zc/dnxaXE7ceJE9e/fX7Vr19aZM2e0YMECpaamaunSpYqMjNSoUaOUnJys6OhoRURE6MEHH1RiYqKuvfZaSVJSUpKaN2+u4cOH64UXXlBGRoaefPJJjR07tsTiFQCAK13zyUtL3Xd9kxqad3dH5+P2Ty/XWZu9xLad6kXrw/sTnY+7Pr9KJ3KKD3/Y/9wN5Y6xTZs2evLJJyWdqxWee+45Va9eXaNHj5YkTZ48WXPnztXOnTu1fPlytWvXTv/4xz+cz3/rrbeUkJCgn376SXFxcXrvvff0zjvvqFevXpKkt99+W7Vq1Sozhj9+2eyZZ57RAw884CxupXM98q+99poaNGggSRo3bpymTZtW7uuFZ/m0uM3MzNRdd92l9PR0RUZGqnXr1lq6dKn69Okj6dyfBMxms4YOHar8/Hz17du3yJvKYrFo8eLFGjNmjBITE1WlShWNGDGCNxYAAJex1q1bO+9bLBZVq1ZNrVq1cm47/+XxzMxMfffdd1q1apVzuMCF0tLSlJOTo4KCAnXq1Mm5PTo6Wk2aNCkzhuXLl2v69On673//q6ysLBUWFiovL0+5ubnO7+eEhoY6C1tJio+PV2ZmpnsXDY/xaXH75ptvlrk/ODhYc+bM0Zw5c0ptU6dOHX311VeeDg0AAL/047S+pe4z/2FM5NZJvUtpWbzt2seuv7TALvDHPz+f/zP1hY+lc3/yzs7O1sCBA/X8888XO058fLx++umncp9///79uvHGGzVmzBg9++yzio6O1tq1azVq1CgVFBQ4i9uS4jw/OwB8p9KNuQUAABWnPGNgK6qtJ1199dX69NNPVbduXQUEFI+hQYMGslqt2rRpk+rWrStJOnnypH766Sd17969xGNu3bpVDodDL730knOc6EcffVRh1wDPYhUBAABw2Ro7dqxOnDih2267Td9++63S0tK0dOlS3X333bLb7QoLC9Odd96pxx57TCtXrtSuXbs0cuTIMr/c1LBhQ9lsNs2ePVt79+7Vu+++q9dee82LV4VLQXHrZQFms/7UtL5OrW3ECmUAAFyimjVrat26dbLb7UpKSlKrVq00fvx4RUVFOQvYadOmqWvXrho4cKB69+6trl27qn379qUes02bNpoxY4aef/55tWzZUu+//76mT5/urUvCJWJYgpcFBph1S/OGemmdVQFmVgoCAOBCqampxbbt37+/2LYLx7Y2atRIn332WYnHMwxDYWFheuedd4r01j7yyCNlnmPChAmaMGFCkW3Dhw933h85cqRGjhxZZP/gwYMZc1sJ0HUIAAAAv0HPrZc5HIYOZWXLWj1ADoO5eAEAADyJnlsvyyu062/L16vmqDUqsJc8MTYAAADcQ3ELAAAAv0FxCwAAAL9BcQsAgB/j2/u4nHji/UpxCwCAH7JYLJKkgoICH0cCuC43N1dS8aWNy4PZEgAA8EMBAQEKDQ3V0aNHZbVay1yRy585HA4VFBQoLy/vis2BL7maf8MwlJubq8zMTEVFRTl/OXMHxS0AAH7IZDIpPj5e+/bt04EDB3wdjs8YhqGzZ88qJCREJpPJ1+Fcccqb/6ioKMXFxV3SOSluvSzAbNbARnX03nsWBQzmN0gAQMUJDAxUo0aNruihCTabTWvWrFG3bt0u6U/dcE958m+1Wi+px/Y8ilsvCwwwa3irJvpXKsvvAgAqntlsVnBwsK/D8BmLxaLCwkIFBwdT3PqAL/JP1yEAAAD8Bj23XuZwGMrMOStLhE0Og/QDAAB4EtWVl+UV2jVu6TeqNUYqsPf0dTgAAAB+hWEJAAAA8BsUtwAAAPAbFLcAAADwGxS3AAAA8BsUtwAAAPAbFLcAAADwGxS3XmYxm5RUP0FnttWRxUT6AQAAPInqysuCAiy6t20znUhpKauF9AMAAHgS1RUAAAD8BiuUeZlhGMrKL5A5xCHDMPk6HAAAAL9CcetlZ2123ftlqhIekvJZfhcAAMCjGJYAAAAAv0FxCwAAAL9BcQsAAAC/QXELAAAAv0FxCwAAAL9BcQsAAAC/QXHrZRazSd1r11T297VYfhcAAMDDqK68LCjAorEdWur4V21YfhcAAMDDqK4AAADgN1ihzMsMw1BeYaFMVpMMw/B1OAAAAH6F4tbLztrsuus/K1U7meV3AQAAPI1hCQAAAPAbFLcAAADwGxS3AAAA8BsUtwAAAPAbFLcAAADwGxS3AAAA8BsUt15mNpl07VWxyvlvnMwmk6/DAQAA8CsUt14WbLUouVMbHfu8vQItFl+HAwAA4FcobgEAAOA3KG4BAADgN1h+18tyCwp1y2fLVOcxKa+wpySrr0MCAADwG/TcAgAAwG9Q3AIAAMBvUNwCAADAb1DcAgAAwG9Q3AIAAMBvUNwCAADAb/i0uJ0+fbquueYahYeHKyYmRoMHD9aePXuKtOnRo4dMJlOR2wMPPFCkzcGDB3XDDTcoNDRUMTExeuSRR1RYWOjNS3GZ2WRSu9jqyk2rwfK7AAAAHubT4nb16tUaO3asNm7cqJSUFNlsNiUlJSknJ6dIu9GjRys9Pd15e+GFF5z77Ha7brjhBhUUFGj9+vV6++23NX/+fE2ePNnbl+OSYKtFE7tcraOfdGT5XQAAAA/z6SIOS5YsKfJ4/vz5iomJ0datW9WtWzfn9tDQUMXFxZV4jGXLlunHH3/U8uXLFRsbq7Zt2+rpp5/WY489pilTpigwMLBCrwEAAACVR6Vaoez06dOSpOjo6CLb33//fb333nuKi4vTwIEDNWnSJIWGhkqSNmzYoFatWik2NtbZvm/fvhozZox++OEHtWvXrth58vPzlZ+f73yclZUlSbLZbLLZbB6/rj+y2QolWWWzFaq0050bVWFVYaGt1DYVFJxzzTSbzSbvnrzinX99vfE6ozjy7zvk3rfIv++Qe9/yZP5dPUalKW4dDofGjx+vLl26qGXLls7tt99+u+rUqaOaNWtq586deuyxx7Rnzx599tlnkqSMjIwiha0k5+OMjIwSzzV9+nRNnTq12PZly5Y5i+aKkm+X/r45QAkTLFq7MUjHjmWV2C4tLVJSD61du07p6acrNKYLWfLydOPv95cuXSp7cLDXzu1NKSkpvg7hikb+fYfc+xb59x1y71ueyH9ubq5L7SpNcTt27Fjt2rVLa9euLbL9vvvuc95v1aqV4uPj1atXL6WlpalBgwZunWvixIlKTk52Ps7KylJCQoKSkpIUERHh3gW4KLegUI9uXilzoF0dOlyjbp2DSmy3ffu5f7t27aISOp8rzgXjnfv27StVqeLFk1c8m82mlJQU9enTR1ar9eJPgEeRf98h975F/n2H3PuWJ/N//i/tF1Mpittx48Zp8eLFWrNmjWrVqlVm206dOkmSfvnlFzVo0EBxcXHavHlzkTZHjhyRpFLH6QYFBSkoqHhRabVaK/yNbzX+N0OC1Wop9XwBAef/tcqrP4sXnMxqtcq7J/ceb7zWKB359x1y71vk33fIvW95Iv+uPt+nsyUYhqFx48Zp4cKFWrlyperVq3fR5+zYsUOSFB8fL0lKTEzU999/r8zMTGeblJQURUREqHnz5hUSNwAAAConn/bcjh07VgsWLNDnn3+u8PBw5xjZyMhIhYSEKC0tTQsWLNCAAQNUrVo17dy5UxMmTFC3bt3UunVrSVJSUpKaN2+u4cOH64UXXlBGRoaefPJJjR07tsTeWQAAAPgvn/bczp07V6dPn1aPHj0UHx/vvH344YeSpMDAQC1fvlxJSUlq2rSp/va3v2no0KH64osvnMewWCxavHixLBaLEhMTdeedd+quu+7StGnTfHVZAAAA8BGf9twahlHm/oSEBK1evfqix6lTp46++uorT4UFAACAy5RPe26vRGaTSc2rV1XewWiW3wUAAPAwilsvC7ZaNKXbNTry/xJZfhcAAMDDKG4BAADgNyhuAQAA4DcqxSIOV5LcgkKNWrxKtR40Ka+wqyQmlAYAAPAUilsfOFNgkyXU11EAAAD4H4YlAAAAwG9Q3AIAAMBvUNwCAADAb1DcAgAAwG9Q3AIAAMBvUNx6mdlkUoOoCOWnR7L8LgAAgIdR3HpZsNWi6T2vVcY7XVl+FwAAwMMobgEAAOA3KG4BAADgN1ihzMvOFtg1dskaXfWASfmFncXyuwAAAJ5Dcetlhgwdzc1TQOS5+wAAAPAchiUAAADAb1DcAgAAwG9Q3AIAAMBvUNwCAADAb1DcAgAAwG9Q3HqZSSbVCq+igmNhMonldwEAADyJ4tbLQgItmtGni9Lf7K6gAJbfBQAA8CS357k9ePCgDhw4oNzcXNWoUUMtWrRQUFCQJ2MDAAAAyqVcxe3+/fs1d+5cffDBB/r1119lGP9bhCAwMFDXXXed7rvvPg0dOlRmM53CAAAA8C6XK9CHHnpIbdq00b59+/TMM8/oxx9/1OnTp1VQUKCMjAx99dVX6tq1qyZPnqzWrVvr22+/rci4L1tnC+xKTlmn+FGrlV9o93U4AAAAfsXlntsqVapo7969qlatWrF9MTEx6tmzp3r27KmnnnpKS5Ys0aFDh3TNNdd4NFh/YMjQr2dyFFid5XcBAAA8zeXidvr06S4ftF+/fm4FAwAAAFwKtwbGnj17Vrm5uc7HBw4c0KxZs7R06VKPBQYAAACUl1vF7aBBg/TOO+9Ikk6dOqVOnTrppZde0uDBgzV37lyPBggAAAC4yq3idtu2bbruuuskSZ988oliY2N14MABvfPOO3rllVc8GiAAAADgKreK29zcXIWHh0uSli1bpiFDhshsNuvaa6/VgQMHPBogAAAA4Cq3ituGDRtq0aJFOnTokJYuXaqkpCRJUmZmpiIiIjwaoL8xyaQaocEqPB3C8rsAAAAe5lZxO3nyZD388MOqW7euOnXqpMTEREnnenHbtWvn0QD9TUigRXP6ddNvr/Vk+V0AAAAPc2v53T/96U/q2rWr0tPT1aZNG+f2Xr166eabb/ZYcAAAAEB5uFXcSlJcXJzi4uKKbOvYseMlBwQAAAC4y+XidsiQIS4f9LPPPnMrmCtBns2uiSs3Ku4ukwrsHSRZfR0SAACA33C5uI2MjKzIOK4YDsNQ2qksBcWfuw8AAADPcbm4nTdvXkXGAQAAAFwyt2ZLAAAAACojt79Q9sknn+ijjz7SwYMHVVBQUGTftm3bLjkwAAAAoLzc6rl95ZVXdPfddys2Nlbbt29Xx44dVa1aNe3du1f9+/f3dIwAAACAS9wqbl999VW9/vrrmj17tgIDA/Xoo48qJSVFDz30kE6fPu3pGAEAAACXuFXcHjx4UJ07d5YkhYSE6MyZM5Kk4cOH6//9v//nuej8VHigVfbcQF+HAQAA4HfcKm7j4uJ04sQJSVLt2rW1ceNGSdK+fftkML1VmUIDA/Tmjdfr19l9FBzg9pBnAAAAlMCt4rZnz576z3/+I0m6++67NWHCBPXp00e33nory+8CAADAZ9zqOnz99dflcDgkSWPHjlW1atW0fv163XTTTbr//vs9GiAAAADgKreKW7PZLLP5f52+w4YN07BhwzwWlD/Ls9k1Zc23ir3NpAJ7O7H8LgAAgOe4XNzu3LlTLVu2lNls1s6dO8ts27p160sOzF85DEM/Hjup4NosvwsAAOBpLhe3bdu2VUZGhmJiYtS2bVuZTKYSvzxmMplkt9s9GiQAAADgCpeL23379qlGjRrO+wAAAEBl43JxW6dOHUmSzWbT1KlTNWnSJNWrV6/CAgMAAADKq9xTgVmtVn366acVEQsAAABwSdya53bw4MFatGiRh0MBAAAALo1bU4E1atRI06ZN07p169S+fXtVqVKlyP6HHnrII8H5qyCLWWfPmnwdBgAAgN9xq7h98803FRUVpa1bt2rr1q1F9plMJorbMoQGBujdQb3VqZNVwcNsvg4HAADAr7g1LGHfvn2l3vbu3evycaZPn65rrrlG4eHhiomJ0eDBg7Vnz54ibfLy8pyroIWFhWno0KE6cuRIkTYHDx7UDTfcoNDQUMXExOiRRx5RYWGhO5cGAACAy5hbxe20adOUm5tbbPvZs2c1bdo0l4+zevVqjR07Vhs3blRKSopsNpuSkpKUk5PjbDNhwgR98cUX+vjjj7V69WodPnxYQ4YMce632+264YYbVFBQoPXr1+vtt9/W/PnzNXnyZHcuDQAAAJcxt4rbqVOnKjs7u9j23NxcTZ061eXjLFmyRCNHjlSLFi3Upk0bzZ8/XwcPHnQOdTh9+rTefPNNzZgxQz179lT79u01b948rV+/Xhs3bpQkLVu2TD/++KPee+89tW3bVv3799fTTz+tOXPmqKCgwJ3Lq1B5Nrumr9umGn/arAIWuwAAAPAot8bcGoYhk6n4F6K+++47RUdHux3M6dOnJcl5jK1bt8pms6l3797ONk2bNlXt2rW1YcMGXXvttdqwYYNatWql2NhYZ5u+fftqzJgx+uGHH9SuXbti58nPz1d+fr7zcVZWlqRzc/jabBU7Dja/oFDbjxxTaINz90s737lRFVYVFtpUwSEVZbPJ6rxrk3dPXvHO57uiX2eUjPz7Drn3LfLvO+TetzyZf1ePUa7itmrVqjKZTDKZTGrcuHGRAtdutys7O1sPPPBA+SL9ncPh0Pjx49WlSxe1bNlSkpSRkaHAwEBFRUUVaRsbG6uMjAxnmwsL2/P7z+8ryfTp00vsYV62bJlCQ0Pdit9V+XbpfNq3bPlW2aeySmyXlhYpqYfWrl2n9PTTFRrThSx5ebrx9/tLly6VPTjYa+f2ppSUFF+HcEUj/75D7n2L/PsOufctT+S/pCGxJSlXcTtr1iwZhqF77rlHU6dOVWRkpHNfYGCg6tatq8TExPJF+ruxY8dq165dWrt2rVvPL4+JEycqOTnZ+TgrK0sJCQlKSkpSREREhZ47t6BQj25eKUnq0OEadescVGK77dvP/du1axeV0PlccS4Y79y3b1/pD9O8Xe5sNptSUlLUp08fWa3Wiz8BHkX+fYfc+xb59x1y71uezP/5v7RfTLmK2xEjRkiS6tWrpy5duiggwK1RDcWMGzdOixcv1po1a1SrVi3n9ri4OBUUFOjUqVNFem+PHDmiuLg4Z5vNmzcXOd752RTOt/mjoKAgBQUVLyqtVmuFv/Gtxv96u61WS6nnO5/agACrvPqzeMHJrFarvHty7/HGa43SkX/fIfe+Rf59h9z7lify7+rz3fpCWffu3T1S2BqGoXHjxmnhwoVauXKl6tWrV2R/+/btZbVatWLFCue2PXv26ODBg84e4sTERH3//ffKzMx0tklJSVFERISaN29+yTECAADg8uGZrlc3jR07VgsWLNDnn3+u8PBw5xjZyMhIhYSEKDIyUqNGjVJycrKio6MVERGhBx98UImJibr22mslSUlJSWrevLmGDx+uF154QRkZGXryySc1duzYEntnAQAA4L98WtzOnTtXktSjR48i2+fNm6eRI0dKkmbOnCmz2ayhQ4cqPz9fffv21auvvupsa7FYtHjxYo0ZM0aJiYmqUqWKRowYUa75dgEAAOAffFrcGoZx0TbBwcGaM2eO5syZU2qbOnXq6KuvvvJkaBUmNDBAHw1JOrf87hCmJQEAAPAkt8bcAgAAAJWRWz23OTk5eu6557RixQplZmbK4XAU2b93716PBAcAAACUh1vF7b333qvVq1dr+PDhio+PL3G1MpQsz2bXjE3fqfogkwrsLSQxLQkAAICnuFXcfv311/ryyy/VpUsXT8fj9xyGoY2/HVGVppLDYKoyAAAAT3JrzG3VqlUVHR3t6VgAAACAS+JWcfv0009r8uTJLq/xCwAAAHiDW8MSXnrpJaWlpSk2NlZ169Ytthzatm3bPBIcAAAAUB5uFbeDBw/2cBgAAADApXOruH3qqac8HQcAAABwyVjEAQAAAH7Dp8vvXolCrBa9c1NP9ehhVdDNjos/AQAAAC6j59bLTCaTggMCZNgCWPwCAADAw1wubrOysioyDgAAAOCSuVzcVq1aVZmZmZKknj176tSpUxUVk1/LL7RrzpZdqjbgO9nsDEsAAADwJJeL27CwMB0/flySlJqaKpvNVmFB+TO7w9Dqg4cV1upX2Q2KWwAAAE9y+QtlvXv31vXXX69mzZpJkm6++WYFBgaW2HblypWeiQ4AAAAoB5eL2/fee09vv/220tLStHr1arVo0UKhoaEVGRsAAABQLi4XtyEhIXrggQckSVu2bNHzzz+vqKioiooLAAAAKDe35rldtWqV875hGJLEtFYAAADwObfnuX3nnXfUqlUrhYSEKCQkRK1bt9a7777rydgAAACAcnGr53bGjBmaNGmSxo0bpy5dukiS1q5dqwceeEDHjh3ThAkTPBokAAAA4Aq3itvZs2dr7ty5uuuuu5zbbrrpJrVo0UJTpkyhuC1DiNWif9/QQ337BijoZoZyAAAAeJJbwxLS09PVuXPnYts7d+6s9PT0Sw7Kn5lMJkUEBcpxNohxygAAAB7mVnHbsGFDffTRR8W2f/jhh2rUqNElBwUAAAC4w61hCVOnTtWtt96qNWvWOMfcrlu3TitWrCix6MX/5Bfa9e8duxXdxyybnV8EAAAAPMmtntuhQ4dq06ZNql69uhYtWqRFixapevXq2rx5s26++WZPx+hX7A5Dy/YeUvjVB1h+FwAAwMPc6rmVpPbt2+u9997zZCwAAADAJXF7nlsAAACgsqG4BQAAgN+guAUAAIDfoLgFAACA33CruL3nnnt05syZYttzcnJ0zz33XHJQAAAAgDvcKm7ffvttnT17ttj2s2fP6p133rnkoPxZcIBF/+p7nX6de70CLRZfhwMAAOBXyjUVWFZWlgzDkGEYOnPmjIKDg5377Ha7vvrqK8XExHg8SH9iNpsUUyVE9iyrzCabr8MBAADwK+UqbqOiomQymWQymdS4ceNi+00mk6ZOneqx4AAAAIDyKFdxu2rVKhmGoZ49e+rTTz9VdHS0c19gYKDq1KmjmjVrejxIf1JQ6NC73+9RVA+LCh31fR0OAACAXylXcdu9e3dJ0r59+1S7dm2ZTKYKCcqfFToc+uLnA4rsJBU66vo6HAAAAL/i1vK7Bw4c0IEDB0rd361bN7cDAgAAANzlVnHbo0ePYtsu7MW12+1uBwQAAAC4y62pwE6ePFnklpmZqSVLluiaa67RsmXLPB0jAAAA4BK3em4jIyOLbevTp48CAwOVnJysrVu3XnJgAAAAQHl5dPnd2NhY7dmzx5OHBAAAAFzmVs/tzp07izw2DEPp6el67rnn1LZtW0/EBQAAAJSbW8Vt27ZtZTKZZBhGke3XXnut3nrrLY8E5q+CAyx6qXdn3XZbgAJvZvldAAAAT3KruN23b1+Rx2azWTVq1CiyHC9KZjablBARJtsxlt8FAADwNLeK2zp16ng6DgAAAOCSuf2FstWrV2vgwIFq2LChGjZsqJtuuknffPONJ2PzSwWFDn304y+K7PKTCh0OX4cDAADgV9wqbt977z317t1boaGheuihh/TQQw8pJCREvXr10oIFCzwdo18pdDj0yX/3KqrrzxS3AAAAHubWsIRnn31WL7zwgiZMmODc9tBDD2nGjBl6+umndfvtt3ssQAAAAMBVbvXc7t27VwMHDiy2/aabbir2ZTMAAADAW9wqbhMSErRixYpi25cvX66EhIRLDgoAAABwh1vDEv72t7/poYce0o4dO9S5c2dJ0rp16zR//ny9/PLLHg0QAAAAcJVbxe2YMWMUFxenl156SR999JEkqVmzZvrwww81aNAgjwYIAAAAuMqt4laSbr75Zt18882ejAUAAAC4JC6Puf3jUrtwT1CARf/o0Unpb3dRoIXldwEAADzJ5eK2RYsW+uCDD1RQUFBmu59//lljxozRc889d8nB+SOL2aSG0ZEqyIiS2WTydTgAAAB+xeXidvbs2frnP/+puLg43XrrrXrxxRf1/vvv69NPP9W///1vJScnq2PHjmrbtq0iIiI0ZsyYix5zzZo1GjhwoGrWrCmTyaRFixYV2T9y5EiZTKYit379+hVpc+LECd1xxx2KiIhQVFSURo0apezsbFcvCwAAAH7E5TG3vXr10pYtW7R27Vp9+OGHev/993XgwAGdPXtW1atXV7t27XTXXXfpjjvuUNWqVV06Zk5Ojtq0aaN77rlHQ4YMKbFNv379NG/ePOfjoKCgIvvvuOMOpaenKyUlRTabTXfffbfuu+++SrtSWkGhQ//5aZ8iOlpU6Kjl63AAAAD8Srm/UNa1a1d17drVIyfv37+/+vfvX2aboKAgxcXFlbhv9+7dWrJkib799lt16NBB0rke5gEDBuif//ynatas6ZE4PanQ4dB7u35W1eulQkfliw8AAOBy5tYiDt6UmpqqmJgYNWnSRGPGjNHx48ed+zZs2KCoqChnYStJvXv3ltls1qZNm3wRLgAAAHzI7anAvKFfv34aMmSI6tWrp7S0NP39739X//79tWHDBlksFmVkZCgmJqbIcwICAhQdHa2MjIxSj5ufn6/8/Hzn46ysLEmSzWaTzWarmIv5nc1WeMF9e6nnKyyUJKsKC22q4JCKstlkdd61ybsnr3jn813RrzNKRv59h9z7Fvn3HXLvW57Mv6vHqNTF7bBhw5z3W7VqpdatW6tBgwZKTU1Vr1693D7u9OnTNXXq1GLbly1bptDQULeP64p8u3Q+7Vu2fKvsU1kltktLi5TUQ2vXrlN6+ukKjelClrw83fj7/aVLl8oeHOy1c3tTSkqKr0O4opF/3yH3vkX+fYfc+5Yn8p+bm+tSu0pd3P5R/fr1Vb16df3yyy/q1auX4uLilJmZWaRNYWGhTpw4Ueo4XUmaOHGikpOTnY+zsrKUkJCgpKQkRUREVFj8kpRbUKhHN6+UJHXocI26dQ4qsd327ef+7dq1i9q1q9CQisrJcd7t27evVKWKF09e8Ww2m1JSUtSnTx9ZrdaLPwEeRf59h9z7Fvn3HXLvW57M//m/tF/MZVXc/vrrrzp+/Lji4+MlSYmJiTp16pS2bt2q9u3bS5JWrlwph8OhTp06lXqcoKCgYrMuSJLVaq3wN77V+N/ctlarpdTzBQSc/9cqr/4sXnAyq9Uq757ce7zxWqN05N93yL1vkX/fIfe+5Yn8u/p8t75Qtm3bNn3//ffOx59//rkGDx6sv//97xdd5OFC2dnZ2rFjh3bs2CFJ2rdvn3bs2KGDBw8qOztbjzzyiDZu3Kj9+/drxYoVGjRokBo2bHiuR1FSs2bN1K9fP40ePVqbN2/WunXrNG7cOA0bNqxSzpQAAACAiuVWcXv//ffrp59+kiTt3btXw4YNU2hoqD7++GM9+uijLh9ny5Ytateundr9/nf35ORktWvXTpMnT5bFYtHOnTt10003qXHjxho1apTat2+vb775pkiv6/vvv6+mTZuqV69eGjBggLp27arXX3/dncvyiqAAi566roMyFlzL8rsAAAAe5tawhJ9++klt27aVJH388cfq1q2bFixYoHXr1mnYsGGaNWuWS8fp0aOHDMModf/SpUsveozo6OhKu2BDSSxmk1rUiFb+IavMJr65CQAA4Elu9dwahiGHwyFJWr58uQYMGCBJSkhI0LFjxzwXHQAAAFAObhW3HTp00DPPPKN3331Xq1ev1g033CDp3JjZ2NhYjwbob2x2h5akHVRYu/0q/P0XBAAAAHiGW8XtzJkztW3bNo0bN05PPPGEGjZsKEn65JNP1LlzZ48G6G9sdofe+u6/qpb0A8UtAACAh7k15rZNmzZFZks478UXX1RAwGU1uxgAAAD8iFs9t/Xr19fx48eLbc/Ly1Pjxo0vOSgAAADAHW4Vt/v375fdbi+2PT8/X7/++uslBwUAAAC4o1xjCP7zn/847y9dulSRkZHOx3a7XStWrFC9evU8Fx0AAABQDuUqbgcPHixJMplMGjFiRJF9VqtVdevW1UsvveSx4AAAAIDyKFdxe35u23r16unbb79V9erVKyQoAAAAwB1ujbndt28fha2bAi1mPZ7YTpkfd5DV7Fb6AQAAUAq35+1asWKFVqxYoczMTGeP7nlvvfXWJQfmrwIsZl0dX0Nn91plMbP8LgAAgCe5VdxOnTpV06ZNU4cOHRQfHy+TyeTpuAAAAIByc6u4fe211zR//nwNHz7c0/H4PZvdodQDv6lKS4sKHTG+DgcAAMCvuDXos6CggGV23WSzO/Tq1h9U/YadLL8LAADgYW4Vt/fee68WLFjg6VgAAACAS+LWsIS8vDy9/vrrWr58uVq3bi2r1Vpk/4wZMzwSHAAAAFAebhW3O3fuVNu2bSVJu3btKrKPL5cBAADAV9wqbletWuXpOAAAAIBLxioCAAAA8Btu9dxef/31ZQ4/WLlypdsBAQAAAO5yq7g9P972PJvNph07dmjXrl0aMWKEJ+LyW4EWsyZ0bK2//z1A1sF0nAMAAHiSW8XtzJkzS9w+ZcoUZWdnX1JA/i7AYlZirTjl7mH5XQAAAE/zaNfhnXfeqbfeesuThwQAAABc5lbPbWk2bNig4OBgTx7S7xTaHdrwa4ZCmwTI7oj2dTgAAAB+xa3idsiQIUUeG4ah9PR0bdmyRZMmTfJIYP6qwO7QzM07VWOwZHP09HU4AAAAfsWt4jYyMrLIY7PZrCZNmmjatGlKSkrySGAAAABAeblV3M6bN8/TcQAAAACX7JLG3G7dulW7d++WJLVo0ULt2rXzSFAAAACAO9wqbjMzMzVs2DClpqYqKipKknTq1Cldf/31+uCDD1SjRg1PxggAAAC4xK2pwB588EGdOXNGP/zwg06cOKETJ05o165dysrK0kMPPeTpGAEAAACXuNVzu2TJEi1fvlzNmjVzbmvevLnmzJnDF8oAAADgM2713DocDlmt1mLbrVarHA7HJQflz6wWs/7SvoWOfdlaAWaW3wUAAPAkt6qrnj176q9//asOHz7s3Pbbb79pwoQJ6tWrl8eC80dWi1k96lylnF0JFLcAAAAe5lZ19a9//UtZWVmqW7euGjRooAYNGqhevXrKysrS7NmzPR0jAAAA4BK3xtwmJCRo27ZtWr58uf773/9Kkpo1a6bevXt7NDh/VGh3aFv6UYXUt8juiPJ1OAAAAH7F7XluTSaT+vTpoz59+ngyHr9XYHfouQ3bFfNnlt8FAADwtHINS1i5cqWaN2+urKysYvtOnz6tFi1a6JtvvvFYcAAAAEB5lKu4nTVrlkaPHq2IiIhi+yIjI3X//fdrxowZHgsOAAAAKI9yFbffffed+vXrV+r+pKQkbd269ZKDAgAAANxRruL2yJEjJc5ve15AQICOHj16yUEBAAAA7ihXcXvVVVdp165dpe7fuXOn4uPjLzkoAAAAwB3lKm4HDBigSZMmKS8vr9i+s2fP6qmnntKNN97oseAAAACA8ijXVGBPPvmkPvvsMzVu3Fjjxo1TkyZNJEn//e9/NWfOHNntdj3xxBMVEqi/sFrMuqdNU734okUBg1mhDAAAwJPKVdzGxsZq/fr1GjNmjCZOnCjDMCSdm/O2b9++mjNnjmJjYyskUH9htZjVr0FtPbXdqgCzzdfhAAAA+JVyL+JQp04dffXVVzp58qR++eUXGYahRo0aqWrVqhURHwAAAOAyt1coq1q1qq655hpPxnJFsDsM/XD0hIISAuQwwn0dDgAAgF9h0KeX5RfaNfWbLYq7faMK7HZfhwMAAOBXKG4BAADgNyhuAQAA4DcobgEAAOA3KG4BAADgNyhuAQAA4DcobgEAAOA3KG69LMBs1p0tG+nkqqYKMJN+AAAAT6K68rLAALNualxPWZsbUNwCAAB4GNUVAAAA/Ibby+/CPXaHoV9OnFZgXIAcRqivwwEAAPArPu25XbNmjQYOHKiaNWvKZDJp0aJFRfYbhqHJkycrPj5eISEh6t27t37++ecibU6cOKE77rhDERERioqK0qhRo5Sdne3Fqyif/EK7/p66SfEj1rH8LgAAgIf5tLjNyclRmzZtNGfOnBL3v/DCC3rllVf02muvadOmTapSpYr69u2rvLw8Z5s77rhDP/zwg1JSUrR48WKtWbNG9913n7cuAQAAAJWIT4cl9O/fX/379y9xn2EYmjVrlp588kkNGjRIkvTOO+8oNjZWixYt0rBhw7R7924tWbJE3377rTp06CBJmj17tgYMGKB//vOfqlmzpteuBQAAAL5Xacfc7tu3TxkZGerdu7dzW2RkpDp16qQNGzZo2LBh2rBhg6KiopyFrST17t1bZrNZmzZt0s0331zisfPz85Wfn+98nJWVJUmy2Wyy2WwVdEX6/RyFF9y3l3q+wkJJsqqw0KYKDqkom01W512bvHvyinc+3xX9OqNk5N93yL1vkX/fIfe+5cn8u3qMSlvcZmRkSJJiY2OLbI+NjXXuy8jIUExMTJH9AQEBio6OdrYpyfTp0zV16tRi25ctW6bQ0Ir9kle+XTqf9i1bvlX2qawS26WlRUrqobVr1yk9/XSFxnQhS16ebvz9/tKlS2UPDvbaub0pJSXF1yFc0ci/75B73yL/vkPufcsT+c/NzXWpXaUtbivSxIkTlZyc7HyclZWlhIQEJSUlKSIiokLPnVtQqEc3r5Qkdehwjbp1Diqx3fbt5/7t2rWL2rWr0JCKyslx3u3bt69UpYoXT17xbDabUlJS1KdPH1mt1os/AR5F/n2H3PsW+fcdcu9bnsz/+b+0X0ylLW7j4uIkSUeOHFF8fLxz+5EjR9S2bVtnm8zMzCLPKyws1IkTJ5zPL0lQUJCCgooXlVartcLf+FbDdMH5LKWeLyDg/L9WefVn8YKTWa1Weffk3uON1xqlI/++Q+59i/z7Drn3LU/k39XnV9pFHOrVq6e4uDitWLHCuS0rK0ubNm1SYmKiJCkxMVGnTp3S1q1bnW1Wrlwph8OhTp06eT1mVwSYzfpT0/o6tbYRK5QBAAB4mE97brOzs/XLL784H+/bt087duxQdHS0ateurfHjx+uZZ55Ro0aNVK9ePU2aNEk1a9bU4MGDJUnNmjVTv379NHr0aL322muy2WwaN26chg0bVmlnSggMMOuW5g310jqrAswMbgcAAPAknxa3W7Zs0fXXX+98fH4c7IgRIzR//nw9+uijysnJ0X333adTp06pa9euWrJkiYIv+JLT+++/r3HjxqlXr14ym80aOnSoXnnlFa9fCwAAAHzPp8Vtjx49ZBhGqftNJpOmTZumadOmldomOjpaCxYsqIjwKoTDYehQVras1QPkMEr+MhkAAADcw6BPL8srtOtvy9er5qg1LL8LAADgYRS3AAAA8BsUtwAAAPAbFLcAAADwGxS3AAAA8BsUtwAAAPAbFLcAAADwGxS3XhZgNmtgozo6vak+y+8CAAB4GNWVlwUGmDW8VROdSm1GcQsAAOBhVFcAAADwGz5dfvdK5HAYysw5K0uETQ6D9AMAAHgS1ZWX5RXaNW7pN6o1Riqw9/R1OAAAAH6FYQkAAADwGxS3AAAA8BsUtwAAAPAbFLcAAADwGxS3AAAA8BsUtwAAAPAbFLdeZjGblFQ/QWe21ZHFRPoBAAA8ierKy4ICLLq3bTOdSGkpq4X0AwAAeBLVFQAAAPwGK5R5mWEYysovkDnEIcMw+TocAAAAv0Jx62VnbXbd+2WqEh6S8ll+FwAAwKMYlgAAAAC/QXELAAAAv0FxCwAAAL9BcQsAAAC/QXELAAAAv0FxCwAAAL9BcetlFrNJ3WvXVPb3tVh+FwAAwMOorrwsKMCisR1a6vhXbVh+FwAAwMOorgAAAOA3WKHMywzDUF5hoUxWkwzD8HU4AAAAfoXi1svO2uy66z8rVTuZ5XcBAAA8jWEJAAAA8BsUtwAAAPAbFLcAAADwGxS3AAAA8BsUtwAAAPAbFLcAAADwGxS3XmY2mXTtVbHK+W+czCaTr8MBAADwKxS3XhZstSi5Uxsd+7y9Ai0WX4cDAADgVyhuAQAA4DcobgEAAOA3WH7Xy3ILCnXLZ8tU5zEpr7CnJKuvQwIAAPAb9NwCAADAb1DcAgAAwG9Q3AIAAMBvUNwCAADAb1DcAgAAwG9Q3AIAAMBvUNx6mdlkUrvY6spNq8HyuwAAAB5GcetlwVaLJna5Wkc/6cjyuwAAAB5GcQsAAAC/QXELAAAAv8Hyu16WW1Co4Z8vV8IEk/IKu4vldwEAADyH4tYH8u0OmQN9HQUAAID/qdTDEqZMmSKTyVTk1rRpU+f+vLw8jR07VtWqVVNYWJiGDh2qI0eO+DBiAAAA+FKlLm4lqUWLFkpPT3fe1q5d69w3YcIEffHFF/r444+1evVqHT58WEOGDPFhtAAAAPClSj8sISAgQHFxccW2nz59Wm+++aYWLFignj17SpLmzZunZs2aaePGjbr22mu9HSoAAAB8rNIXtz///LNq1qyp4OBgJSYmavr06apdu7a2bt0qm82m3r17O9s2bdpUtWvX1oYNG8osbvPz85Wfn+98nJWVJUmy2Wyy2WwVdzGSbLbCC+7bSz1fYaEkWVVYaFMFh1SUzeb8ipvNZpN3T17xzue7ol9nlIz8+w659y3y7zvk3rc8mX9Xj1Gpi9tOnTpp/vz5atKkidLT0zV16lRdd9112rVrlzIyMhQYGKioqKgiz4mNjVVGRkaZx50+fbqmTp1abPuyZcsUGhrqyUsoJt8unU/7li3fKvtUVont0tIiJfXQ2rXrlJ5+ukJjupAlL083/n5/6dKlsgcHe+3c3pSSkuLrEK5o5N93yL1vkX/fIfe+5Yn85+bmutSuUhe3/fv3d95v3bq1OnXqpDp16uijjz5SSEiI28edOHGikpOTnY+zsrKUkJCgpKQkRUREXFLMF5Nns2te2lZt325WxwGt1DUxqMR227ef+7dr1y5q165CQyoqJ8d5t2/fvlKVKl48ecWz2WxKSUlRnz59ZLUyDZu3kX/fIfe+Rf59h9z7lifzf/4v7RdTqYvbP4qKilLjxo31yy+/qE+fPiooKNCpU6eK9N4eOXKkxDG6FwoKClJQUPGi0mq1Vvgb32q1amr3jur0qFVVxttKPV9AwPl/rfLqz+IFJ7NarfLuyb3HG681Skf+fYfc+xb59x1y71ueyL+rz6/0syVcKDs7W2lpaYqPj1f79u1ltVq1YsUK5/49e/bo4MGDSkxM9GGUAAAA8JVK3XP78MMPa+DAgapTp44OHz6sp556ShaLRbfddpsiIyM1atQoJScnKzo6WhEREXrwwQeVmJjITAkAAABXqEpd3P7666+67bbbdPz4cdWoUUNdu3bVxo0bVaNGDUnSzJkzZTabNXToUOXn56tv37569dVXfRx12XILCjVq8SrVetCkvMKuYvldAAAAz6nUxe0HH3xQ5v7g4GDNmTNHc+bM8VJEnnGmwCZLxU7KAAAAcEW6rMbcAgAAAGWhuAUAAIDfoLgFAACA36C4BQAAgN+guAUAAIDfoLj1MrPJpAZREcpPj5TZZPJ1OAAAAH6F4tbLgq0WTe95rTLe6apAi8XX4QAAAPgVilsAAAD4DYpbAAAA+I1KvUKZPzpbYNfYJWt01QMm5Rd2FsvvAgAAeA7FrZcZMnQ0N08BkefuAwAAwHMYlgAAAAC/QXELAAAAv0FxCwAAAL9BcQsAAAC/QXELAAAAv0Fx62UmmVQrvIoKjoXJJJbfBQAA8CSKWy8LCbRoRp8uSn+zu4ICWH4XAADAkyhuAQAA4DcobgEAAOA3WKHMy84W2JWcsk7xo0zKL+wklt8FAADwHIpbLzNk6NczOQqszvK7AAAAnsawBAAAAPgNilsAAAD4DYpbAAAA+A2KWwAAAPgNilsAAAD4DYpbLzPJpBqhwSo8HcLyuwAAAB5GcetlIYEWzenXTb+91pPldwEAADyM4hYAAAB+g+IWAAAAfoMVyrwsz2bXxJUbFXeXSQX2DrrY8ru7d1/8mNWrS7VreyY+AACAyxnFrZc5DENpp7IUFH/ufmmqV5dCQ6U777z4MUNDzxXBFLgAAOBKR3FbSdWufa5gPXas7Ha7d58rgI8do7gFAACguK3EatemYAUAACgPitsrxMGDF+8FlqQaoVJCxYcDAABQIShurwAHD0rNmkm5uRdvWz1EOlrxIQEAAFQIitsrwLFj5wrb9947V+SWZvdu6T4XvsAGAABQWVHc+kB4oFWnTnl26d2ypgw7v69ZM+nqqz16WgAAgEqF4tbLQgMD9OaN16tTJ6uC77Rd8vFcnTIsNPRcWwAAAH9GcXuZc3XKMBZ6AAAAVwKKWz/AlGEAAADnUNx6WZ7NrilrvlXsbSYV2NvpYsvv+tKOHZIjpOR93u4JdnUqM3qoAQC4slHcepnDMPTjsZMKrl328ru+UL26FBoi6ey5x126SqXNHubNJX/LM5UZSxEDAHBlo7iFU+3a0rZtkn6fLmzd2pJ7br295G95pjJjKWIAAK5sFLcoIuGC5cnatpVUxVeRFMdUZgAA4GLMvg4AAAAA8BSKWwAAAPgNhiWgwrg6w8HFlLX6GgAAwIUobn0gyGLW2bOeXX63sinPDAeuYIU1AADgCopbLwsNDNC7g3qfW3532KUvv1tZuTrDgavKM39tWT29hYXS0aOlTN4LAAAuexS3qFDenOGgevVzPbx33llWK6uCgnrq+usNNWjgnbgAAID3UNzCb9Sufa7Xtqxxvt9/X6iRIwO0dm2hTp8u+3isdgYAwOWH4tbL8mx2TV+3TTX+ZFKBvbUq8/K7l6PatcsuSCMjDQUFnStwL4bVzvyHK19u5JcZAPAPFLde5jAMbT9yTKENKt/yu1eC2rWlf/1rpVq3vl4BAaX/YnF+tbNvvil7zLCrBZGrM0dczgWWK9foizHPrn650ZO/zFwJrzcAVFYUt7ji1KhxVu3aSdYyOs1dG7/rWkFUnpkjLtfeYtev0ftjnl35cqMnl26+El5vAKjMKG6BErgyftfVgsjVmSNc7S12lSd7BS/WE7l7t2vXeH7M8/HjNq9/oc9bX24s7+vtiYJauvhrxEwh8Fe89/FHflPczpkzRy+++KIyMjLUpk0bzZ49Wx07dvR1WH6trCm3/GHhhYuN3z3vYtd6fv/FiitXe4tdFRoqffaZVKNG2ee82DWW58/6111X9vEKCyv3UBxX3reu/tLgajF9sXN67jWq3DOFMC4a7vCH9z48zy+K2w8//FDJycl67bXX1KlTJ82aNUt9+/bVnj17FBMT4+vw/E55/mTvzwsvlKcYdSUXrvQWu+roUWnIEKlfv4vHdbE/i7vaE1m+uYhNCvDSp4+rBWt5Xsuyfmlw9Re78vwcXeyXFFd6zsvTa+7JMcOuHOv8+7Uyjov2lMJCKS0tUtu366LvfW8X8d7OheS5a3Tl88mXfzGqjK6E7wT4RXE7Y8YMjR49Wnfffbck6bXXXtOXX36pt956S48//riPo/M/rhZhl/MPhivKU4y6mgtXe4td4alhFed54s/61arJ5dkqPOliv1y4+lqW55cGT/wy4+r5zp+zrJ5zV3vNPTlmuLzHWrKk7F8afDEu2nOsknq41NKbY7F9kwvPX2NZn0+V/S9G3nSlfCfgsi9uCwoKtHXrVk2cONG5zWw2q3fv3tqwYUOJz8nPz1d+fr7z8enfJzw9ceKEbLaKXTUst6BQjvxz76pTp07q+HEvf6JcTE6Oc3Iy2/HjUl5eic2qVDl3u5jjxz0XmifYbDbl5ubq+PHjspb1jTIXuZoHyfu5uFhsp05JklWbNxfq1KnSP/x/+skkKUCnTtku+RqqVLHpxRc3qEmTzgrwVtetpOjoc7koK35XXss6daR166QTJy79fK6c09XzuXLOU6cKJYXo22/zlZ1tLvU4P/1kUm5ugObOLVTjxmW/L8aMCdBXX5XeztVjnY8/IaH0/a6+X11Rnrg8pbCwUFu2bFGHDh3KfO+7kldP8kUuPHmNrnw+ufrevxJ48udbkmJjpbi4ss/pyf93z5w5I0kyLjbblHGZ++233wxJxvr164tsf+SRR4yOHTuW+JynnnrKkMSNGzdu3Lhx48btMrsdOnSozNrwsu+5dcfEiROVnJzsfOxwOHTixAlVq1ZNJpOpws+flZWlhIQEHTp0SBERERV+PvwPufct8u875N63yL/vkHvf8mT+DcPQmTNnVLNmzTLbXfbFbfXq1WWxWHTkyJEi248cOaK4UvrKg4KCFBQUVGRbVFRURYVYqoiICH7QfITc+xb59x1y71vk33fIvW95Kv+RkZEXbXPZDz4JDAxU+/bttWLFCuc2h8OhFStWKDEx0YeRAQAAwNsu+55bSUpOTtaIESPUoUMHdezYUbNmzVJOTo5z9gQAAABcGfyiuL311lt19OhRTZ48WRkZGWrbtq2WLFmi2NhYX4dWoqCgID311FPFhkag4pF73yL/vkPufYv8+w659y1f5N9kGBebTwEAAAC4PFz2Y24BAACA8yhuAQAA4DcobgEAAOA3KG4BAADgNyhuK8icOXNUt25dBQcHq1OnTtq8eXOZ7T/++GM1bdpUwcHBatWqlb766isvRep/ypP7+fPny2QyFbkFBwd7MVr/sWbNGg0cOFA1a9aUyWTSokWLLvqc1NRUXX311QoKClLDhg01f/78Co/TX5U3/6mpqcXe+yaTSRkZGd4J2I9Mnz5d11xzjcLDwxUTE6PBgwdrz549F30en/ue4U7++ez3jLlz56p169bOBRoSExP19ddfl/kcb7zvKW4rwIcffqjk5GQ99dRT2rZtm9q0aaO+ffsqMzOzxPbr16/XbbfdplGjRmn79u0aPHiwBg8erF27dnk58stfeXMvnVs1JT093Xk7cOCAFyP2Hzk5OWrTpo3mzJnjUvt9+/bphhtu0PXXX68dO3Zo/Pjxuvfee7V06dIKjtQ/lTf/5+3Zs6fI+z8mJqaCIvRfq1ev1tixY7Vx40alpKTIZrMpKSlJOTk5pT6Hz33PcSf/Ep/9nlCrVi0999xz2rp1q7Zs2aKePXtq0KBB+uGHH0ps77X3vQGP69ixozF27FjnY7vdbtSsWdOYPn16ie1vueUW44YbbiiyrVOnTsb9999foXH6o/Lmft68eUZkZKSXortySDIWLlxYZptHH33UaNGiRZFtt956q9G3b98KjOzK4Er+V61aZUgyTp486ZWYriSZmZmGJGP16tWltuFzv+K4kn8++ytO1apVjX//+98l7vPW+56eWw8rKCjQ1q1b1bt3b+c2s9ms3r17a8OGDSU+Z8OGDUXaS1Lfvn1LbY+SuZN7ScrOzladOnWUkJBQ5m+c8Cze95VD27ZtFR8frz59+mjdunW+DscvnD59WpIUHR1dahve/xXHlfxLfPZ7mt1u1wcffKCcnBwlJiaW2MZb73uKWw87duyY7HZ7sdXRYmNjSx3LlpGRUa72KJk7uW/SpIneeustff7553rvvffkcDjUuXNn/frrr94I+YpW2vs+KytLZ8+e9VFUV474+Hi99tpr+vTTT/Xpp58qISFBPXr00LZt23wd2mXN4XBo/Pjx6tKli1q2bFlqOz73K4ar+eez33O+//57hYWFKSgoSA888IAWLlyo5s2bl9jWW+97v1h+F3BXYmJikd8wO3furGbNmun//u//9PTTT/swMqBiNWnSRE2aNHE+7ty5s9LS0jRz5ky9++67Pozs8jZ27Fjt2rVLa9eu9XUoVyRX889nv+c0adJEO3bs0OnTp/XJJ59oxIgRWr16dakFrjfQc+th1atXl8Vi0ZEjR4psP3LkiOLi4kp8TlxcXLnao2Tu5P6PrFar2rVrp19++aUiQsQFSnvfR0REKCQkxEdRXdk6duzIe/8SjBs3TosXL9aqVatUq1atMtvyue955cn/H/HZ777AwEA1bNhQ7du31/Tp09WmTRu9/PLLJbb11vue4tbDAgMD1b59e61YscK5zeFwaMWKFaWOQUlMTCzSXpJSUlJKbY+SuZP7P7Lb7fr+++8VHx9fUWHid7zvK58dO3bw3neDYRgaN26cFi5cqJUrV6pevXoXfQ7vf89xJ/9/xGe/5zgcDuXn55e4z2vve49+PQ2GYRjGBx98YAQFBRnz5883fvzxR+O+++4zoqKijIyMDMMwDGP48OHG448/7my/bt06IyAgwPjnP/9p7N6923jqqacMq9VqfP/99766hMtWeXM/depUY+nSpUZaWpqxdetWY9iwYUZwcLDxww8/+OoSLltnzpwxtm/fbmzfvt2QZMyYMcPYvn27ceDAAcMwDOPxxx83hg8f7my/d+9eIzQ01HjkkUeM3bt3G3PmzDEsFouxZMkSX13CZa28+Z85c6axaNEi4+effza+//57469//athNpuN5cuX++oSLltjxowxIiMjjdTUVCM9Pd15y83Ndbbhc7/iuJN/Pvs94/HHHzdWr15t7Nu3z9i5c6fx+OOPGyaTyVi2bJlhGL5731PcVpDZs2cbtWvXNgIDA42OHTsaGzdudO7r3r27MWLEiCLtP/roI6Nx48ZGYGCg0aJFC+PLL7/0csT+ozy5Hz9+vLNtbGysMWDAAGPbtm0+iPryd35qqT/ezud7xIgRRvfu3Ys9p23btkZgYKBRv359Y968eV6P21+UN//PP/+80aBBAyM4ONiIjo42evToYaxcudI3wV/mSsq7pCLvZz73K447+eez3zPuueceo06dOkZgYKBRo0YNo1evXs7C1jB89743GYZheLYvGAAAAPANxtwCAADAb1DcAgAAwG9Q3AIAAMBvUNwCAADAb1DcAgAAwG9Q3AIAAMBvUNwCAADAb1DcAsAl6NGjh8aPH1+u5yxatEgNGzaUxWJxPrekba4wmUxatGiRJGn//v0ymUzasWNHueK5HPTo0UMmk8kj1zdy5Ejnsc7nDoD/oLgF4FdGjhypwYMH+zqMMt1///3605/+pEOHDunpp58udVt5JSQkKD09XS1btrxo28uxEB49erTL11eWl19+Wenp6R6KCkBlE+DrAADgSpKdna3MzEz17dtXNWvWLHWbOywWi+Li4jwVaqUTGhrqkeuLjIxUZGSkByICUBnRcwvgsvPJJ5+oVatWCgkJUbVq1dS7d2/l5ORoypQpevvtt/X55587/+ycmpoqSXrsscfUuHFjhYaGqn79+po0aZJsNpvzmFOmTFHbtm317rvvqm7duoqMjNSwYcN05swZZ5ucnBzdddddCgsLU3x8vF566aVisZ08eVJ33XWXqlatqtDQUPXv318///yzJCk1NVXh4eGSpJ49ezrjK2lbSX7++Wd169ZNwcHBat68uVJSUors/2Nv7MmTJ3XHHXeoRo0aCgkJUaNGjTRv3jxJUr169SRJ7dq1k8lkUo8ePSRJ3377rfr06aPq1asrMjJS3bt317Zt24qcx2Qy6d///rduvvlmhYaGqlGjRvrPf/5TpM0PP/ygG2+8UREREQoPD9d1112ntLQ05/5///vfatasmYKDg9W0aVO9+uqrJV5zWVJTU2UymbR06VK1a9dOISEh6tmzpzIzM/X111+rWbNmioiI0O23367c3NxyHx/A5YniFsBlJT09Xbfddpvuuece7d69W6mpqRoyZIgMw9DDDz+sW265Rf369VN6errS09PVuXNnSVJ4eLjmz5+vH3/8US+//LLeeOMNzZw5s8ix09LStGjRIi1evFiLFy/W6tWr9dxzzzn3P/LII1q9erU+//xzLVu2TKmpqcUKv5EjR2rLli36z3/+ow0bNsgwDA0YMEA2m02dO3fWnj17JEmffvqpM76Stv2Rw+HQkCFDFBgYqE2bNum1117TY489VmauJk2apB9//FFff/21du/erblz56p69eqSpM2bN0uSli9frvT0dH322WeSpDNnzmjEiBFau3atNm7cqEaNGmnAgAFFinxJmjp1qm655Rbt3LlTAwYM0B133KETJ05Ikn777Td169ZNQUFBWrlypbZu3ap77rlHhYWFkqT3339fkydP1rPPPqvdu3frH//4hyZNmqS33367zOspzZQpU/Svf/1L69ev16FDh3TLLbdo1qxZWrBggb788kstW7ZMs2fPduvYAC5DBgBcRrZu3WpIMvbv31/i/hEjRhiDBg266HFefPFFo3379s7HTz31lBEaGmpkZWU5tz3yyCNGp06dDMMwjDNnzhiBgYHGRx995Nx//PhxIyQkxPjrX/9qGIZh/PTTT4YkY926dc42x44dM0JCQpzPO3nypCHJWLVqlbNNSdv+aOnSpUZAQIDx22+/Obd9/fXXhiRj4cKFhmEYxr59+wxJxvbt2w3DMIyBAwcad999d4nH+2Pb0tjtdiM8PNz44osvnNskGU8++aTzcXZ2tiHJ+Prrrw3DMIyJEyca9erVMwoKCko8ZoMGDYwFCxYU2fb0008biYmJpcbRvXt3Z57PW7VqlSHJWL58uXPb9OnTDUlGWlqac9v9999v9O3bt9gxL8wdAP/BmFsAl5U2bdqoV69eatWqlfr27aukpCT96U9/UtWqVct83ocffqhXXnlFaWlpys7OVmFhoSIiIoq0qVu3rnOIgCTFx8crMzNT0rle3YKCAnXq1Mm5Pzo6Wk2aNHE+3r17twICAoq0qVatmpo0aaLdu3df0nXv3r1bCQkJRcbkJiYmlvmcMWPGaOjQodq2bZuSkpI0ePDgEnuFL3TkyBE9+eSTSk1NVWZmpux2u3Jzc3Xw4MEi7Vq3bu28X6VKFUVERDhztWPHDl133XWyWq3Fjp+Tk6O0tDSNGjVKo0ePdm4vLCx0exzshbHExsY6h55cuO18TzUA/8ewBACXFYvFopSUFH399ddq3ry5Zs+erSZNmmjfvn2lPmfDhg264447NGDAAC1evFjbt2/XE088oYKCgiLt/liMmUwmORyOCrkOb+jfv78OHDigCRMm6PDhw+rVq5cefvjhMp8zYsQI7dixQy+//LLWr1+vHTt2qFq1auXKVUhISKnHz87OliS98cYb2rFjh/O2a9cubdy40Z3LLBKLyWTyu9cRQPlQ3AK47JhMJnXp0kVTp07V9u3bFRgYqIULF0qSAgMDZbfbi7Rfv3696tSpoyeeeEIdOnRQo0aNdODAgXKds0GDBrJardq0aZNz28mTJ/XTTz85Hzdr1kyFhYVF2hw/flx79uxR8+bN3bnUIsc+dOhQkSmsXCkGa9SooREjRui9997TrFmz9Prrr0s6lydJxXK1bt06PfTQQxowYIBatGihoKAgHTt2rFyxtm7dWt98802RL+ydFxsbq5o1a2rv3r1q2LBhkdv5L7kBwKVgWAKAy8qmTZu0YsUKJSUlKSYmRps2bdLRo0fVrFkzSeeGFixdulR79uxRtWrVFBkZqUaNGungwYP64IMPdM011+jLL790FsOuCgsL06hRo/TII4+oWrVqiomJ0RNPPCGz+X99BI0aNdKgQYM0evRo/d///Z/Cw8P1+OOP66qrrtKgQYMu6bp79+6txo0ba8SIEXrxxReVlZWlJ554osznTJ48We3bt1eLFi2Un5+vxYsXO/MUExOjkJAQLVmyRLVq1VJwcLAzV++++646dOigrKwsPfLII2X2xJZk3Lhxmj17toYNG6aJEycqMjJSGzduVMeOHdWkSRNNnTpVDz30kCIjI9WvXz/l5+dry5YtOnnypJKTk93OEQBI9NwCuMxERERozZo1GjBggBo3bqwnn3xSL730kvr37y/p3ET/TZo0UYcOHVSjRg2tW7dON910kyZMmKBx48apbdu2Wr9+vSZNmlTuc7/44ou67rrrNHDgQPXu3Vtdu3ZV+/bti7SZN2+e2rdvrxtvvFGJiYkyDENfffVVieNPy8NsNmvhwoU6e/asOnbsqHvvvVfPPvtsmc8JDAzUxIkT1bp1a3Xr1k0Wi0UffPCBJCkgIECvvPKK/u///k81a9Z0Ft9vvvmmTp48qauvvlrDhw/XQw89pJiYmHLFWq1aNa1cuVLZ2dnq3r272rdvrzfeeMOZg3vvvVf//ve/NW/ePLVq1Urdu3fX/Pnz6bkF4BEmwzAMXwcBAEBZevToobZt22rWrFkeO6bJZNLChQsr/Yp2AMqHnlsAwGXh1VdfVVhYmL7//vtLOs4DDzygsLAwD0UFoLKh5xYAUOn99ttvOnv2rCSpdu3azi/EuSMzM1NZWVmSzk33VqVKFY/ECKByoLgFAACA32BYAgAAAPwGxS0AAAD8BsUtAAAA/AbFLQAAAPwGxS0AAAD8BsUtAAAA/AbFLQAAAPwGxS0AAAD8BsUtAAAA/Mb/B0rlaU+/pz1QAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "bin_edges = np.arange(results_df['x'].min(), results_df['x'].max() + 0.1, 0.05)\n", + "plt.figure(figsize=(8,6))\n", + "plt.hist(results_df['x'], bins=bin_edges, edgecolor='blue', histtype='step', linewidth=1, label = 'Standoff Distance')\n", + "plt.ylabel('Counts (out of n trials)')\n", + "plt.xlabel('standoff distance [m]')\n", + "plt.axvline(np.mean(results_df['x']), color = 'red', label = 'mean' )\n", + "plt.axvline(np.median(results_df['x']), linestyle='--', label = 'median')\n", + "\n", + "plt.legend()\n", + "plt.grid(True)\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "monte-carlo", + "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.10.9" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tests/data/correlated_samples_arrhenius.csv b/tests/data/correlated_samples_arrhenius.csv new file mode 100644 index 00000000..e9b7312e --- /dev/null +++ b/tests/data/correlated_samples_arrhenius.csv @@ -0,0 +1,51 @@ +Ea,LnR0,X +74.07708998694397,9.73022090701182,-0.01429488294411331 +57.56168948006327,15.2070907191989,0.15753724881478248 +58.179029072132565,14.938657735666816,0.11616413357559288 +54.155268350478956,16.347485864196095,0.09921572043947152 +68.4717276686662,11.566588203430717,-0.053853851412172214 +45.0812954925816,19.45784411812305,0.019095068155474897 +74.96683072815007,9.474552112619692,0.07108497196921254 +56.457878071368945,15.67691261652651,-0.09238424398306455 +64.43635895565852,12.955942496207507,-0.000844905990239779 +60.23820028079894,14.408001544809919,-6.505623837697178e-05 +72.87883680142677,10.04881802131346,-0.054920984313326725 +46.86421274779222,18.91318473804258,0.005471338600242955 +59.69869101459703,14.559472872160235,0.08852465969029857 +59.24345134729001,14.648421257566195,-0.032262746260674956 +70.45379434720107,10.957708679512502,0.048666688908438546 +53.95642307787203,16.435451715160227,-0.0842192298475958 +60.80647974467399,14.23706841161549,-0.01644726793751712 +55.59631329691633,16.01132817058643,-0.176717207411224 +62.391782290492024,13.788872823630117,0.045484819833095054 +66.38455660546232,12.172312391059052,0.13007652506635406 +53.95104688094081,16.330225251277138,0.08839670557759077 +70.53470037613342,10.852972849724695,-0.012066261476659824 +68.73896874415426,11.505991614722141,0.1449330954339707 +65.79132268826142,12.548606810790648,0.17698799838276 +68.73354187007709,11.519976481160565,-0.14425761121317024 +57.03012277771021,15.25441467161328,0.22715474085393308 +61.17235737236437,14.002160199713693,0.19334898345535437 +55.16859411244937,16.10038556969974,0.02803166329161076 +60.10143222149817,14.402551888847054,-0.08437780913910284 +65.9970994062349,12.47080337755871,0.08209823684289543 +56.97153201990721,15.414783557746738,0.02507054561407481 +59.14965780134712,14.687427023045222,-0.01322509570440275 +57.00467987145666,15.435665533961899,-0.0864853224211003 +55.837480173018754,15.843818423107802,0.06467088077195328 +57.12231032686542,15.397211274452863,0.09572980809899696 +61.98646180530478,13.762923688080797,-0.027135893374793676 +53.827769227049565,16.4319749499913,0.10573526787300028 +63.811347460937476,13.172333892654557,-0.08414852905999287 +74.33896691969808,9.628620056268247,0.10592132197778342 +67.56058976119208,11.976223383154856,-0.0062179580528424205 +60.66314097736758,14.290307775579663,-0.03262197763903248 +55.52414999706222,15.931105803345623,0.014918887960663107 +56.56163827401506,15.540017757474615,0.12601300646588615 +74.58013119227073,9.488421192542235,0.1326771818322532 +62.4552559152248,13.62982508677269,0.06460749508418803 +57.37527755356807,15.303080703617074,0.04164837733383297 +63.490063586656966,13.223453269989704,0.05579429488419243 +77.59206438700544,8.53365394103253,0.09375079055971544 +62.967469991238815,13.376775366357666,0.08094152780854805 +66.63853872767706,12.251102912815893,0.06863948562286684 diff --git a/tests/data/h5_pytest.h5 b/tests/data/h5_pytest.h5 index ddc427af..86eee1ac 100644 Binary files a/tests/data/h5_pytest.h5 and b/tests/data/h5_pytest.h5 differ diff --git a/tests/data/monte_carlo_arrhenius.csv b/tests/data/monte_carlo_arrhenius.csv new file mode 100644 index 00000000..87c83430 --- /dev/null +++ b/tests/data/monte_carlo_arrhenius.csv @@ -0,0 +1,51 @@ +0 +9.49800132810379e-10 +0.0001801309368080377 +0.00010856798451580064 +0.0023110060076571622 +5.9213709986490384e-08 +2.19655591255509 +4.99458768992465e-10 +0.0004955624232773357 +1.2105164789706175e-06 +2.8689618592820938e-05 +2.1543921203454037e-09 +0.6168242487142045 +4.0343672743271745e-05 +5.542332995072857e-05 +1.3903376050724424e-08 +0.002937269884460099 +1.928499281769252e-05 +0.0010197958817409255 +6.312444928090319e-06 +2.396018051136527e-07 +0.0024789537465739887 +1.2344632582145316e-08 +4.697551145164206e-08 +4.381579926732448e-07 +5.236549815638251e-08 +0.00022923180311753554 +1.2240710290375742e-05 +0.0012246041838437723 +3.111592735634318e-05 +3.838461126222035e-07 +0.0002955579124060609 +5.946475918096224e-05 +0.00031054673377995426 +0.0007111660380219962 +0.00026623920665945483 +7.441460023868276e-06 +0.0028680558449895527 +1.9965595690662264e-06 +7.447524502891329e-10 +1.273198256175642e-07 +2.1691344179317994e-05 +0.0008984041627532367 +0.00038205161762326563 +5.824818415191378e-10 +5.213079393742761e-06 +0.0002227757975378916 +2.283375516333262e-06 +6.661537463037927e-11 +3.266806849587072e-06 +2.3827314137386823e-07 diff --git a/tests/data/noCorrEmpty_samples_arrhenius.csv b/tests/data/noCorrEmpty_samples_arrhenius.csv new file mode 100644 index 00000000..8d5d71a7 --- /dev/null +++ b/tests/data/noCorrEmpty_samples_arrhenius.csv @@ -0,0 +1,51 @@ +Ea,LnR0,X +74.07708998694397,14.464733976474415,-0.010289001259125324 +57.56168948006327,12.851072045304905,0.15566385954995926 +58.179029072132565,10.89646361986343,0.07415691518293091 +54.155268350478956,12.858262374021516,0.09302792339932608 +68.4717276686662,13.205640324161344,-0.07459813992423395 +45.0812954925816,15.173231532409613,0.05091555960960241 +74.96683072815007,15.797406113824323,0.10761926007209843 +56.457878071368945,16.025247609879408,-0.060579294834903964 +64.43635895565852,14.428665159778518,0.007670971463536215 +60.23820028079894,15.911570280527151,0.037337831953165854 +72.87883680142677,11.856419972663234,-0.10221718314962083 +46.86421274779222,16.82108699520686,0.06538766925657294 +59.69869101459703,14.990962201850891,0.11810319060145749 +59.24345134729001,12.98502116595036,-0.051229046687254025 +70.45379434720107,14.930583643917144,0.06890069741406432 +53.95642307787203,13.535393275902877,-0.0961778542489839 +60.80647974467399,16.521221365337983,0.03025847623075629 +55.59631329691633,17.481343857715505,-0.12630693155381786 +62.391782290492024,19.127996348636866,0.1454295279774913 +66.38455660546232,10.26828737263206,0.07469388713405513 +53.95104688094081,10.150512811920201,0.03165613447356299 +70.53470037613342,12.474588908063561,-0.04285471235615173 +68.73896874415426,14.118135720834273,0.16055301159566365 +65.79132268826142,15.889378803377674,0.22938540312848824 +68.73354187007709,14.502983377112344,-0.15035245018025642 +57.03012277771021,8.720701633460473,0.15682104943802505 +61.17235737236437,12.964959403510765,0.19568616787438486 +55.16859411244937,15.770177594510447,0.06765634778561869 +60.10143222149817,14.2914126890975,-0.08495817339769823 +65.9970994062349,15.607027015439254,0.119809210740905 +56.97153201990721,13.172413595382832,0.016139010366980157 +59.14965780134712,13.225763887940387,-0.02585464306861507 +57.00467987145666,14.1837395910409,-0.08801488245903258 +55.837480173018754,14.736508706704784,0.08875499528928596 +57.12231032686542,14.212772560652184,0.11280645657894939 +61.98646180530478,14.016658162768165,-0.027801449660472086 +53.827769227049565,12.06352736331704,0.08578058026639483 +63.811347460937476,14.656154930051413,-0.07950529249381022 +74.33896691969808,14.023614762854171,0.11370531524175004 +67.56058976119208,16.515914848410635,0.038723001148165945 +60.66314097736758,16.687649204721627,0.015578155295057614 +55.52414999706222,14.18026460302728,0.023999107282558337 +56.56163827401506,12.794098224344301,0.12035928145706742 +74.58013119227073,12.142506003028426,0.10859764083574658 +62.4552559152248,14.769757195057345,0.08666304070084224 +57.37527755356807,13.913597281714553,0.0477703840133375 +63.490063586656966,12.87183869551266,0.04182574694832342 +77.59206438700544,13.83013858645015,0.09549013337608692 +62.967469991238815,12.188830726222818,0.05718106010866063 +66.63853872767706,15.448784339959484,0.10186076870245656 diff --git a/tests/data/noCorrR0_arrhenius.csv b/tests/data/noCorrR0_arrhenius.csv new file mode 100644 index 00000000..8d5d71a7 --- /dev/null +++ b/tests/data/noCorrR0_arrhenius.csv @@ -0,0 +1,51 @@ +Ea,LnR0,X +74.07708998694397,14.464733976474415,-0.010289001259125324 +57.56168948006327,12.851072045304905,0.15566385954995926 +58.179029072132565,10.89646361986343,0.07415691518293091 +54.155268350478956,12.858262374021516,0.09302792339932608 +68.4717276686662,13.205640324161344,-0.07459813992423395 +45.0812954925816,15.173231532409613,0.05091555960960241 +74.96683072815007,15.797406113824323,0.10761926007209843 +56.457878071368945,16.025247609879408,-0.060579294834903964 +64.43635895565852,14.428665159778518,0.007670971463536215 +60.23820028079894,15.911570280527151,0.037337831953165854 +72.87883680142677,11.856419972663234,-0.10221718314962083 +46.86421274779222,16.82108699520686,0.06538766925657294 +59.69869101459703,14.990962201850891,0.11810319060145749 +59.24345134729001,12.98502116595036,-0.051229046687254025 +70.45379434720107,14.930583643917144,0.06890069741406432 +53.95642307787203,13.535393275902877,-0.0961778542489839 +60.80647974467399,16.521221365337983,0.03025847623075629 +55.59631329691633,17.481343857715505,-0.12630693155381786 +62.391782290492024,19.127996348636866,0.1454295279774913 +66.38455660546232,10.26828737263206,0.07469388713405513 +53.95104688094081,10.150512811920201,0.03165613447356299 +70.53470037613342,12.474588908063561,-0.04285471235615173 +68.73896874415426,14.118135720834273,0.16055301159566365 +65.79132268826142,15.889378803377674,0.22938540312848824 +68.73354187007709,14.502983377112344,-0.15035245018025642 +57.03012277771021,8.720701633460473,0.15682104943802505 +61.17235737236437,12.964959403510765,0.19568616787438486 +55.16859411244937,15.770177594510447,0.06765634778561869 +60.10143222149817,14.2914126890975,-0.08495817339769823 +65.9970994062349,15.607027015439254,0.119809210740905 +56.97153201990721,13.172413595382832,0.016139010366980157 +59.14965780134712,13.225763887940387,-0.02585464306861507 +57.00467987145666,14.1837395910409,-0.08801488245903258 +55.837480173018754,14.736508706704784,0.08875499528928596 +57.12231032686542,14.212772560652184,0.11280645657894939 +61.98646180530478,14.016658162768165,-0.027801449660472086 +53.827769227049565,12.06352736331704,0.08578058026639483 +63.811347460937476,14.656154930051413,-0.07950529249381022 +74.33896691969808,14.023614762854171,0.11370531524175004 +67.56058976119208,16.515914848410635,0.038723001148165945 +60.66314097736758,16.687649204721627,0.015578155295057614 +55.52414999706222,14.18026460302728,0.023999107282558337 +56.56163827401506,12.794098224344301,0.12035928145706742 +74.58013119227073,12.142506003028426,0.10859764083574658 +62.4552559152248,14.769757195057345,0.08666304070084224 +57.37527755356807,13.913597281714553,0.0477703840133375 +63.490063586656966,12.87183869551266,0.04182574694832342 +77.59206438700544,13.83013858645015,0.09549013337608692 +62.967469991238815,12.188830726222818,0.05718106010866063 +66.63853872767706,15.448784339959484,0.10186076870245656 diff --git a/tests/test_montecarlo.py b/tests/test_montecarlo.py new file mode 100644 index 00000000..414dd579 --- /dev/null +++ b/tests/test_montecarlo.py @@ -0,0 +1,99 @@ +# TODO: +# correlation list is empty AND correlation list is populated with r = 0's + +import pytest +import os +import json +import pandas as pd +import pvdeg +from pvdeg import TEST_DATA_DIR + +WEATHER = pd.read_csv( + os.path.join(TEST_DATA_DIR, r"weather_day_pytest.csv"), + index_col=0, + parse_dates=True +) + +with open(os.path.join(TEST_DATA_DIR, "meta.json"), "r") as file: + META = json.load(file) + + +CORRELATED_SAMPLES_1 = pd.read_csv( + os.path.join(TEST_DATA_DIR, r"correlated_samples_arrhenius.csv"), +) + +CORRELATED_SAMPLES_2 = pd.read_csv( + os.path.join(TEST_DATA_DIR, r"noCorrEmpty_samples_arrhenius.csv") +) + +CORRELATED_SAMPLES_3 = pd.read_csv( + os.path.join(TEST_DATA_DIR, r"noCorrR0_arrhenius.csv") +) + +ARRHENIUS_RESULT = pd.read_csv( + os.path.join(TEST_DATA_DIR, r"monte_carlo_arrhenius.csv"), +) + +def test_generateCorrelatedSamples(): + """ + test pvdeg.montecarlo.generateCorrelatedSamples + + Requires: + --------- + list of correlations, stats dictionary (mean and standard deviation for each variable), number of iterations, seed, DataFrame to check against + """ + # standard case + result_1 = pvdeg.montecarlo.generateCorrelatedSamples( + corr=[pvdeg.montecarlo.Corr('Ea', 'X', 0.0269), pvdeg.montecarlo.Corr('Ea', 'LnR0', -0.9995), pvdeg.montecarlo.Corr('X', 'LnR0', -0.0400)], + stats={'Ea' : {'mean' : 62.08, 'stdev' : 7.3858 }, 'LnR0' : {'mean' : 13.7223084 , 'stdev' : 2.47334772}, 'X' : {'mean' : 0.0341 , 'stdev' : 0.0992757}}, + n = 50, + seed = 1 + ) + + # EMPTY CORRELATION LIST + result_2 = pvdeg.montecarlo.generateCorrelatedSamples( + corr=[], + stats = {'Ea' : {'mean' : 62.08, 'stdev' : 7.3858 }, 'LnR0' : {'mean' : 13.7223084 , 'stdev' : 2.47334772}, 'X' : {'mean' : 0.0341 , 'stdev' : 0.0992757 }}, + n = 50, + seed = 1 + ) + + # populated correlation list, ALL R = 0 + result_3 = pvdeg.montecarlo.generateCorrelatedSamples( + corr=[pvdeg.montecarlo.Corr('Ea', 'X', 0), pvdeg.montecarlo.Corr('Ea', 'LnR0', 0), pvdeg.montecarlo.Corr('X', 'LnR0', 0)], + stats = {'Ea' : {'mean' : 62.08, 'stdev' : 7.3858 }, 'LnR0' : {'mean' : 13.7223084 , 'stdev' : 2.47334772}, 'X' : {'mean' : 0.0341 , 'stdev' : 0.0992757 }}, + n = 50, + seed = 1 + ) + + pd.testing.assert_frame_equal(result_1, CORRELATED_SAMPLES_1) + pd.testing.assert_frame_equal(result_2, CORRELATED_SAMPLES_2) + pd.testing.assert_frame_equal(result_3, CORRELATED_SAMPLES_3) + +def test_simulate(): + """ + test pvdeg.montecarlo.simulate + + Requires: + --------- + target function, correlated samples dataframe, weather dataframe, meta dictionary + """ + + poa_irradiance = pvdeg.spectral.poa_irradiance(WEATHER, META) + temp_mod = pvdeg.temperature.module(weather_df=WEATHER, meta=META, poa=poa_irradiance, conf='open_rack_glass_polymer') + + poa_global = poa_irradiance['poa_global'].to_numpy() + cell_temperature = temp_mod.to_numpy() + + function_kwargs = {'poa_global': poa_global, 'module_temp': cell_temperature} + + new_results = pvdeg.montecarlo.simulate( + func=pvdeg.degradation.vecArrhenius, + correlated_samples=CORRELATED_SAMPLES_1, + **function_kwargs) + + new_results_df = pd.DataFrame(new_results) + + new_results_df.columns = ARRHENIUS_RESULT.columns + + pd.testing.assert_frame_equal(new_results_df, ARRHENIUS_RESULT) \ No newline at end of file