lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
24234ed5a5ef29270d9f1bc09542d3d6d488d7b4
FYS-STK4155/doc/Programs/DimRed/Bayes.ipynb
T
mhjensen 918b257206 update on intro
2020-08-19 11:44:46 +02:00

1521 lines
101 KiB
Plaintext
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import GPy\n",
"import pandas as pd\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.stats as stats \n",
"import matplotlib.mlab as mlab\n",
"import seaborn as sns\n",
"from sklearn.decomposition import PCA\n",
"from sklearn.preprocessing import StandardScaler, MinMaxScaler, Normalizer, RobustScaler\n",
"from IPython.display import display\n",
"from sklearn.gaussian_process import GaussianProcessRegressor\n",
"from sklearn.gaussian_process.kernels import RBF\n",
"from scipy.stats import norm\n",
"from scipy.stats import uniform"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"#Read csv file and get the outputs(features) and inputs\n",
"features = ['SN112EK42.5', 'SN112EK47.5', 'SN112EK52.5', 'SN112EK57.5', 'SN112EK62.5', 'SN112EK67.5', 'SN112EK72.5', 'SN112EK77.5', 'SN112EK82.5', 'SN112EK87.5', 'SN124EK42.5', 'SN124EK47.5', 'SN124EK52.5', 'SN124EK57.5', 'SN124EK62.5', 'SN124EK67.5', 'SN124EK72.5', 'SN124EK77.5', 'SN124EK82.5', 'SN124EK87.5', 'DREK12.5', 'DREK17.5', 'DREK22.5', 'DREK27.5', 'DREK32.5', 'DREK37.5', 'DREK42.5', 'DREK47.5', 'DREK52.5', 'DREK60.0', 'DREK70.0', 'DREK80.0', 'DREK87.5']\n",
"inputs = ['S0', 'L', 'ms', 'mv']\n",
"error_features = ['SN112EK42.5_Error', 'SN112EK47.5_Error', 'SN112EK52.5_Error', 'SN112EK57.5_Error', 'SN112EK62.5_Error', 'SN112EK67.5_Error', 'SN112EK72.5_Error', 'SN112EK77.5_Error', 'SN112EK82.5_Error', 'SN112EK87.5_Error', 'SN124EK42.5_Error', 'SN124EK47.5_Error', 'SN124EK52.5_Error', 'SN124EK57.5_Error', 'SN124EK62.5_Error', 'SN124EK67.5_Error', 'SN124EK72.5_Error', 'SN124EK77.5_Error', 'SN124EK82.5_Error', 'SN124EK87.5_Error', 'DREK12.5_Error', 'DREK17.5_Error', 'DREK22.5_Error', 'DREK27.5_Error', 'DREK32.5_Error', 'DREK37.5_Error', 'DREK42.5_Error', 'DREK47.5_Error', 'DREK52.5_Error', 'DREK60.0_Error', 'DREK70.0_Error', 'DREK80.0_Error', 'DREK87.5_Error'\n",
"]\n",
"df = pd.read_csv(\"C:/Users/danny/OneDrive/Desktop/bayesian_example/e120_bugfix_model_new_mv.csv\", usecols=features)\n",
"df2 = pd.read_csv(\"C:/Users/danny/OneDrive/Desktop/bayesian_example/e120_bugfix_model_new_mv.csv\", usecols=inputs)\n",
"df3 = pd.read_csv(\"C:/Users/danny/OneDrive/Desktop/bayesian_example/e120_exp_result.csv\", usecols=features)\n",
"df4 = pd.read_csv(\"C:/Users/danny/OneDrive/Desktop/bayesian_example/e120_exp_result.csv\", usecols=error_features)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>SN112EK42.5</th>\n",
" <th>SN112EK47.5</th>\n",
" <th>SN112EK52.5</th>\n",
" <th>SN112EK57.5</th>\n",
" <th>SN112EK62.5</th>\n",
" <th>SN112EK67.5</th>\n",
" <th>SN112EK72.5</th>\n",
" <th>SN112EK77.5</th>\n",
" <th>SN112EK82.5</th>\n",
" <th>SN112EK87.5</th>\n",
" <th>...</th>\n",
" <th>DREK27.5</th>\n",
" <th>DREK32.5</th>\n",
" <th>DREK37.5</th>\n",
" <th>DREK42.5</th>\n",
" <th>DREK47.5</th>\n",
" <th>DREK52.5</th>\n",
" <th>DREK60.0</th>\n",
" <th>DREK70.0</th>\n",
" <th>DREK80.0</th>\n",
" <th>DREK87.5</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1.049740</td>\n",
" <td>1.009370</td>\n",
" <td>0.971587</td>\n",
" <td>0.933543</td>\n",
" <td>0.895236</td>\n",
" <td>0.860123</td>\n",
" <td>0.828202</td>\n",
" <td>0.797841</td>\n",
" <td>0.769038</td>\n",
" <td>0.742806</td>\n",
" <td>...</td>\n",
" <td>1.25778</td>\n",
" <td>1.28492</td>\n",
" <td>1.31206</td>\n",
" <td>1.33919</td>\n",
" <td>1.36544</td>\n",
" <td>1.39081</td>\n",
" <td>1.42918</td>\n",
" <td>1.47076</td>\n",
" <td>1.50127</td>\n",
" <td>1.51709</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1.101460</td>\n",
" <td>1.076550</td>\n",
" <td>1.053560</td>\n",
" <td>1.032650</td>\n",
" <td>1.013830</td>\n",
" <td>0.995491</td>\n",
" <td>0.977648</td>\n",
" <td>0.959740</td>\n",
" <td>0.941767</td>\n",
" <td>0.928342</td>\n",
" <td>...</td>\n",
" <td>1.23401</td>\n",
" <td>1.26444</td>\n",
" <td>1.29513</td>\n",
" <td>1.32608</td>\n",
" <td>1.36158</td>\n",
" <td>1.40161</td>\n",
" <td>1.45320</td>\n",
" <td>1.51812</td>\n",
" <td>1.59659</td>\n",
" <td>1.65511</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1.033870</td>\n",
" <td>0.999617</td>\n",
" <td>0.967117</td>\n",
" <td>0.937034</td>\n",
" <td>0.909368</td>\n",
" <td>0.884840</td>\n",
" <td>0.863449</td>\n",
" <td>0.842654</td>\n",
" <td>0.822455</td>\n",
" <td>0.805631</td>\n",
" <td>...</td>\n",
" <td>1.25178</td>\n",
" <td>1.27908</td>\n",
" <td>1.30598</td>\n",
" <td>1.33249</td>\n",
" <td>1.36163</td>\n",
" <td>1.39341</td>\n",
" <td>1.43365</td>\n",
" <td>1.48005</td>\n",
" <td>1.51807</td>\n",
" <td>1.54192</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>0.979704</td>\n",
" <td>0.948603</td>\n",
" <td>0.921646</td>\n",
" <td>0.900123</td>\n",
" <td>0.884032</td>\n",
" <td>0.867298</td>\n",
" <td>0.849921</td>\n",
" <td>0.837539</td>\n",
" <td>0.830154</td>\n",
" <td>0.821927</td>\n",
" <td>...</td>\n",
" <td>1.23993</td>\n",
" <td>1.26528</td>\n",
" <td>1.29176</td>\n",
" <td>1.31937</td>\n",
" <td>1.34603</td>\n",
" <td>1.37173</td>\n",
" <td>1.40500</td>\n",
" <td>1.44910</td>\n",
" <td>1.48531</td>\n",
" <td>1.50366</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1.026250</td>\n",
" <td>0.998031</td>\n",
" <td>0.974097</td>\n",
" <td>0.951106</td>\n",
" <td>0.929057</td>\n",
" <td>0.908939</td>\n",
" <td>0.890750</td>\n",
" <td>0.876159</td>\n",
" <td>0.865163</td>\n",
" <td>0.855486</td>\n",
" <td>...</td>\n",
" <td>1.21375</td>\n",
" <td>1.23669</td>\n",
" <td>1.26262</td>\n",
" <td>1.29155</td>\n",
" <td>1.32245</td>\n",
" <td>1.35531</td>\n",
" <td>1.40611</td>\n",
" <td>1.47931</td>\n",
" <td>1.54623</td>\n",
" <td>1.59240</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>5 rows × 33 columns</p>\n",
"</div>"
],
"text/plain": [
" SN112EK42.5 SN112EK47.5 SN112EK52.5 SN112EK57.5 SN112EK62.5 \\\n",
"0 1.049740 1.009370 0.971587 0.933543 0.895236 \n",
"1 1.101460 1.076550 1.053560 1.032650 1.013830 \n",
"2 1.033870 0.999617 0.967117 0.937034 0.909368 \n",
"3 0.979704 0.948603 0.921646 0.900123 0.884032 \n",
"4 1.026250 0.998031 0.974097 0.951106 0.929057 \n",
"\n",
" SN112EK67.5 SN112EK72.5 SN112EK77.5 SN112EK82.5 SN112EK87.5 ... \\\n",
"0 0.860123 0.828202 0.797841 0.769038 0.742806 ... \n",
"1 0.995491 0.977648 0.959740 0.941767 0.928342 ... \n",
"2 0.884840 0.863449 0.842654 0.822455 0.805631 ... \n",
"3 0.867298 0.849921 0.837539 0.830154 0.821927 ... \n",
"4 0.908939 0.890750 0.876159 0.865163 0.855486 ... \n",
"\n",
" DREK27.5 DREK32.5 DREK37.5 DREK42.5 DREK47.5 DREK52.5 DREK60.0 \\\n",
"0 1.25778 1.28492 1.31206 1.33919 1.36544 1.39081 1.42918 \n",
"1 1.23401 1.26444 1.29513 1.32608 1.36158 1.40161 1.45320 \n",
"2 1.25178 1.27908 1.30598 1.33249 1.36163 1.39341 1.43365 \n",
"3 1.23993 1.26528 1.29176 1.31937 1.34603 1.37173 1.40500 \n",
"4 1.21375 1.23669 1.26262 1.29155 1.32245 1.35531 1.40611 \n",
"\n",
" DREK70.0 DREK80.0 DREK87.5 \n",
"0 1.47076 1.50127 1.51709 \n",
"1 1.51812 1.59659 1.65511 \n",
"2 1.48005 1.51807 1.54192 \n",
"3 1.44910 1.48531 1.50366 \n",
"4 1.47931 1.54623 1.59240 \n",
"\n",
"[5 rows x 33 columns]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y = df.loc[:, features].values\n",
"#scaling the data\n",
"#y_scaler = MinMaxScaler().fit_transform(y)\n",
"#y_scaler = Normalizer().fit_transform(y)\n",
"#y_scaler = RobustScaler().fit_transform(y)\n",
"#y_scaler = StandardScaler().fit_transform(y)\n",
"pd.DataFrame(data = y, columns = features).head()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>PCA 1</th>\n",
" <th>PCA 2</th>\n",
" <th>PCA 3</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0.004585</td>\n",
" <td>0.189818</td>\n",
" <td>-0.017025</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>-0.889633</td>\n",
" <td>0.023971</td>\n",
" <td>0.013887</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>-0.174100</td>\n",
" <td>0.061279</td>\n",
" <td>-0.026724</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-0.011284</td>\n",
" <td>-0.096832</td>\n",
" <td>-0.013462</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>-0.320008</td>\n",
" <td>-0.111562</td>\n",
" <td>0.062695</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" PCA 1 PCA 2 PCA 3\n",
"0 0.004585 0.189818 -0.017025\n",
"1 -0.889633 0.023971 0.013887\n",
"2 -0.174100 0.061279 -0.026724\n",
"3 -0.011284 -0.096832 -0.013462\n",
"4 -0.320008 -0.111562 0.062695"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#applying PCA on the outputs and reducing it from 33D to 3D\n",
"pca = PCA(n_components=3)\n",
"yPCA = pca.fit_transform(y)\n",
"PCADf = pd.DataFrame(data = yPCA, columns = ['PCA 1', 'PCA 2', 'PCA 3'])\n",
"PCADf.head()\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.95664316 0.03929512 0.00237127]\n"
]
}
],
"source": [
"#Seeing how much each PCA takes into account\n",
"print(pca.explained_variance_ratio_) #it seems like the 3PCA account for 98% of the data"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1.05226638 1.00984326 0.96921742 ... 1.46182261 1.48709511 1.50213919]\n",
" [1.10304362 1.0782522 1.05524568 ... 1.52386133 1.59330332 1.64266425]\n",
" [1.03278436 0.99840936 0.96670162 ... 1.47629668 1.51377543 1.53713317]\n",
" ...\n",
" [0.95791579 0.90589328 0.85745636 ... 1.38459667 1.38846038 1.38896563]\n",
" [1.01619208 0.95436387 0.89447893 ... 1.3714765 1.37153888 1.37266782]\n",
" [0.99489351 0.96558251 0.94040069 ... 1.45655175 1.51113185 1.55003265]]\n"
]
}
],
"source": [
"#testing if the reverse pca works\n",
"G = np.dot(yPCA, pca.components_) + pca.mean_\n",
"print(G)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>S0</th>\n",
" <th>L</th>\n",
" <th>ms</th>\n",
" <th>mv</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>31.159</td>\n",
" <td>43.5</td>\n",
" <td>0.980</td>\n",
" <td>0.865000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>35.073</td>\n",
" <td>56.1</td>\n",
" <td>0.700</td>\n",
" <td>0.875000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>32.395</td>\n",
" <td>66.9</td>\n",
" <td>0.860</td>\n",
" <td>0.855000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>28.275</td>\n",
" <td>86.7</td>\n",
" <td>0.972</td>\n",
" <td>1.095001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>27.657</td>\n",
" <td>79.5</td>\n",
" <td>0.612</td>\n",
" <td>0.975000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" S0 L ms mv\n",
"0 31.159 43.5 0.980 0.865000\n",
"1 35.073 56.1 0.700 0.875000\n",
"2 32.395 66.9 0.860 0.855000\n",
"3 28.275 86.7 0.972 1.095001\n",
"4 27.657 79.5 0.612 0.975000"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#extracting and looking at the inputs\n",
"x = df2.loc[:, inputs].values\n",
"xdf = pd.DataFrame(data = x, columns = ['S0', 'L', 'ms', 'mv'])\n",
"xdf.head()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,\n",
" kernel=RBF(length_scale=1), n_restarts_optimizer=10,\n",
" normalize_y=False, optimizer='fmin_l_bfgs_b',\n",
" random_state=None)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#using sckitlearn to do the gaussian process\n",
"kernel = RBF()\n",
"gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10)\n",
"gp.fit(x, yPCA)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>S0</th>\n",
" <th>L</th>\n",
" <th>ms</th>\n",
" <th>mv</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>67.435407</td>\n",
" <td>96.000057</td>\n",
" <td>1.060497</td>\n",
" <td>0.897444</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>52.220715</td>\n",
" <td>40.647893</td>\n",
" <td>0.619497</td>\n",
" <td>0.998561</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>112.420615</td>\n",
" <td>80.647919</td>\n",
" <td>0.754445</td>\n",
" <td>0.877165</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>111.576292</td>\n",
" <td>93.796187</td>\n",
" <td>0.954161</td>\n",
" <td>0.648290</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>102.520674</td>\n",
" <td>56.963868</td>\n",
" <td>0.880245</td>\n",
" <td>1.001980</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" S0 L ms mv\n",
"0 67.435407 96.000057 1.060497 0.897444\n",
"1 52.220715 40.647893 0.619497 0.998561\n",
"2 112.420615 80.647919 0.754445 0.877165\n",
"3 111.576292 93.796187 0.954161 0.648290\n",
"4 102.520674 56.963868 0.880245 1.001980"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#creating new randomized inputs x_ similar to x\n",
"msv_min = min(x[:,2].min(),x[:,3].min())\n",
"msv_max = max(x[:,2].max(),x[:,3].max())\n",
"sl_min = min(x[:,0].min(),x[:,1].min())\n",
"sl_max = max(x[:,0].max(),x[:,1].max())\n",
"x_sl = np.random.uniform(sl_min,sl_max,94)\n",
"x_msv = np.random.uniform(msv_min,msv_max,94)\n",
"x_sl = x_sl.reshape(47,2)\n",
"x_msv = x_msv.reshape(47,2)\n",
"x_ = np.concatenate((x_sl,x_msv),axis=1)\n",
"x_df = pd.DataFrame(data = x_, columns = ['S0', 'L', 'ms', 'mv'])\n",
"x_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(47,)\n",
"(47, 3)\n"
]
}
],
"source": [
"#using the gp model trained on the inputs(x) and the outputs(y) on x_ to predict y_pred and y_std\n",
"y_pred, y_std = gp.predict(x_, return_std=True)\n",
"#y_pred, y_cov= gp.predict(x_, return_cov=True) #should I inverse pca on y_cov?\n",
"y_preddf = pd.DataFrame(data = y_pred, columns = ['PCA 1', 'PCA 2', 'PCA 3'])\n",
"y_preddf.head()\n",
"print(y_std.shape)\n",
"print(y_pred.shape)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>SN112EK42.5</th>\n",
" <th>SN112EK47.5</th>\n",
" <th>SN112EK52.5</th>\n",
" <th>SN112EK57.5</th>\n",
" <th>SN112EK62.5</th>\n",
" <th>SN112EK67.5</th>\n",
" <th>SN112EK72.5</th>\n",
" <th>SN112EK77.5</th>\n",
" <th>SN112EK82.5</th>\n",
" <th>SN112EK87.5</th>\n",
" <th>...</th>\n",
" <th>DREK27.5</th>\n",
" <th>DREK32.5</th>\n",
" <th>DREK37.5</th>\n",
" <th>DREK42.5</th>\n",
" <th>DREK47.5</th>\n",
" <th>DREK52.5</th>\n",
" <th>DREK60.0</th>\n",
" <th>DREK70.0</th>\n",
" <th>DREK80.0</th>\n",
" <th>DREK87.5</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1.007239</td>\n",
" <td>0.972153</td>\n",
" <td>0.940416</td>\n",
" <td>0.911383</td>\n",
" <td>0.885053</td>\n",
" <td>0.860981</td>\n",
" <td>0.839165</td>\n",
" <td>0.819277</td>\n",
" <td>0.801317</td>\n",
" <td>0.784839</td>\n",
" <td>...</td>\n",
" <td>1.241723</td>\n",
" <td>1.267337</td>\n",
" <td>1.293529</td>\n",
" <td>1.320298</td>\n",
" <td>1.346925</td>\n",
" <td>1.373411</td>\n",
" <td>1.410665</td>\n",
" <td>1.455663</td>\n",
" <td>1.494245</td>\n",
" <td>1.519706</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1.007240</td>\n",
" <td>0.972153</td>\n",
" <td>0.940416</td>\n",
" <td>0.911383</td>\n",
" <td>0.885053</td>\n",
" <td>0.860981</td>\n",
" <td>0.839165</td>\n",
" <td>0.819277</td>\n",
" <td>0.801317</td>\n",
" <td>0.784839</td>\n",
" <td>...</td>\n",
" <td>1.241723</td>\n",
" <td>1.267337</td>\n",
" <td>1.293529</td>\n",
" <td>1.320298</td>\n",
" <td>1.346925</td>\n",
" <td>1.373411</td>\n",
" <td>1.410665</td>\n",
" <td>1.455663</td>\n",
" <td>1.494245</td>\n",
" <td>1.519706</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1.007239</td>\n",
" <td>0.972153</td>\n",
" <td>0.940416</td>\n",
" <td>0.911383</td>\n",
" <td>0.885053</td>\n",
" <td>0.860981</td>\n",
" <td>0.839165</td>\n",
" <td>0.819277</td>\n",
" <td>0.801317</td>\n",
" <td>0.784839</td>\n",
" <td>...</td>\n",
" <td>1.241723</td>\n",
" <td>1.267337</td>\n",
" <td>1.293529</td>\n",
" <td>1.320298</td>\n",
" <td>1.346925</td>\n",
" <td>1.373411</td>\n",
" <td>1.410665</td>\n",
" <td>1.455663</td>\n",
" <td>1.494245</td>\n",
" <td>1.519706</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1.007239</td>\n",
" <td>0.972153</td>\n",
" <td>0.940416</td>\n",
" <td>0.911383</td>\n",
" <td>0.885053</td>\n",
" <td>0.860981</td>\n",
" <td>0.839165</td>\n",
" <td>0.819277</td>\n",
" <td>0.801317</td>\n",
" <td>0.784839</td>\n",
" <td>...</td>\n",
" <td>1.241723</td>\n",
" <td>1.267337</td>\n",
" <td>1.293529</td>\n",
" <td>1.320298</td>\n",
" <td>1.346925</td>\n",
" <td>1.373411</td>\n",
" <td>1.410665</td>\n",
" <td>1.455663</td>\n",
" <td>1.494245</td>\n",
" <td>1.519706</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1.007239</td>\n",
" <td>0.972153</td>\n",
" <td>0.940416</td>\n",
" <td>0.911383</td>\n",
" <td>0.885053</td>\n",
" <td>0.860981</td>\n",
" <td>0.839165</td>\n",
" <td>0.819277</td>\n",
" <td>0.801317</td>\n",
" <td>0.784839</td>\n",
" <td>...</td>\n",
" <td>1.241723</td>\n",
" <td>1.267337</td>\n",
" <td>1.293529</td>\n",
" <td>1.320298</td>\n",
" <td>1.346925</td>\n",
" <td>1.373411</td>\n",
" <td>1.410665</td>\n",
" <td>1.455663</td>\n",
" <td>1.494245</td>\n",
" <td>1.519706</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>5 rows × 33 columns</p>\n",
"</div>"
],
"text/plain": [
" SN112EK42.5 SN112EK47.5 SN112EK52.5 SN112EK57.5 SN112EK62.5 \\\n",
"0 1.007239 0.972153 0.940416 0.911383 0.885053 \n",
"1 1.007240 0.972153 0.940416 0.911383 0.885053 \n",
"2 1.007239 0.972153 0.940416 0.911383 0.885053 \n",
"3 1.007239 0.972153 0.940416 0.911383 0.885053 \n",
"4 1.007239 0.972153 0.940416 0.911383 0.885053 \n",
"\n",
" SN112EK67.5 SN112EK72.5 SN112EK77.5 SN112EK82.5 SN112EK87.5 ... \\\n",
"0 0.860981 0.839165 0.819277 0.801317 0.784839 ... \n",
"1 0.860981 0.839165 0.819277 0.801317 0.784839 ... \n",
"2 0.860981 0.839165 0.819277 0.801317 0.784839 ... \n",
"3 0.860981 0.839165 0.819277 0.801317 0.784839 ... \n",
"4 0.860981 0.839165 0.819277 0.801317 0.784839 ... \n",
"\n",
" DREK27.5 DREK32.5 DREK37.5 DREK42.5 DREK47.5 DREK52.5 DREK60.0 \\\n",
"0 1.241723 1.267337 1.293529 1.320298 1.346925 1.373411 1.410665 \n",
"1 1.241723 1.267337 1.293529 1.320298 1.346925 1.373411 1.410665 \n",
"2 1.241723 1.267337 1.293529 1.320298 1.346925 1.373411 1.410665 \n",
"3 1.241723 1.267337 1.293529 1.320298 1.346925 1.373411 1.410665 \n",
"4 1.241723 1.267337 1.293529 1.320298 1.346925 1.373411 1.410665 \n",
"\n",
" DREK70.0 DREK80.0 DREK87.5 \n",
"0 1.455663 1.494245 1.519706 \n",
"1 1.455663 1.494245 1.519706 \n",
"2 1.455663 1.494245 1.519706 \n",
"3 1.455663 1.494245 1.519706 \n",
"4 1.455663 1.494245 1.519706 \n",
"\n",
"[5 rows x 33 columns]"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#Reverse PCA on y_pred to look like y\n",
"Y = np.dot(y_pred, pca.components_) + pca.mean_\n",
"Ydf= pd.DataFrame(data = Y, columns = features)\n",
"Ydf.head()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>SN112EK42.5</th>\n",
" <th>SN112EK47.5</th>\n",
" <th>SN112EK52.5</th>\n",
" <th>SN112EK57.5</th>\n",
" <th>SN112EK62.5</th>\n",
" <th>SN112EK67.5</th>\n",
" <th>SN112EK72.5</th>\n",
" <th>SN112EK77.5</th>\n",
" <th>SN112EK82.5</th>\n",
" <th>SN112EK87.5</th>\n",
" <th>...</th>\n",
" <th>DREK27.5</th>\n",
" <th>DREK32.5</th>\n",
" <th>DREK37.5</th>\n",
" <th>DREK42.5</th>\n",
" <th>DREK47.5</th>\n",
" <th>DREK52.5</th>\n",
" <th>DREK60.0</th>\n",
" <th>DREK70.0</th>\n",
" <th>DREK80.0</th>\n",
" <th>DREK87.5</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>SN112EK42.5</th>\n",
" <td>0.000313</td>\n",
" <td>0.000341</td>\n",
" <td>0.000364</td>\n",
" <td>0.000384</td>\n",
" <td>0.000400</td>\n",
" <td>0.000414</td>\n",
" <td>0.000428</td>\n",
" <td>0.000439</td>\n",
" <td>0.000448</td>\n",
" <td>0.000458</td>\n",
" <td>...</td>\n",
" <td>-0.000031</td>\n",
" <td>-0.000027</td>\n",
" <td>-0.000015</td>\n",
" <td>0.000005</td>\n",
" <td>0.000031</td>\n",
" <td>0.000063</td>\n",
" <td>0.000118</td>\n",
" <td>0.000197</td>\n",
" <td>0.000295</td>\n",
" <td>0.000373</td>\n",
" </tr>\n",
" <tr>\n",
" <th>SN112EK47.5</th>\n",
" <td>0.000341</td>\n",
" <td>0.000381</td>\n",
" <td>0.000415</td>\n",
" <td>0.000446</td>\n",
" <td>0.000472</td>\n",
" <td>0.000497</td>\n",
" <td>0.000520</td>\n",
" <td>0.000541</td>\n",
" <td>0.000559</td>\n",
" <td>0.000578</td>\n",
" <td>...</td>\n",
" <td>-0.000037</td>\n",
" <td>-0.000030</td>\n",
" <td>-0.000014</td>\n",
" <td>0.000010</td>\n",
" <td>0.000042</td>\n",
" <td>0.000082</td>\n",
" <td>0.000147</td>\n",
" <td>0.000245</td>\n",
" <td>0.000364</td>\n",
" <td>0.000458</td>\n",
" </tr>\n",
" <tr>\n",
" <th>SN112EK52.5</th>\n",
" <td>0.000364</td>\n",
" <td>0.000415</td>\n",
" <td>0.000461</td>\n",
" <td>0.000502</td>\n",
" <td>0.000539</td>\n",
" <td>0.000575</td>\n",
" <td>0.000608</td>\n",
" <td>0.000639</td>\n",
" <td>0.000666</td>\n",
" <td>0.000693</td>\n",
" <td>...</td>\n",
" <td>-0.000042</td>\n",
" <td>-0.000033</td>\n",
" <td>-0.000013</td>\n",
" <td>0.000015</td>\n",
" <td>0.000053</td>\n",
" <td>0.000100</td>\n",
" <td>0.000175</td>\n",
" <td>0.000291</td>\n",
" <td>0.000431</td>\n",
" <td>0.000538</td>\n",
" </tr>\n",
" <tr>\n",
" <th>SN112EK57.5</th>\n",
" <td>0.000384</td>\n",
" <td>0.000446</td>\n",
" <td>0.000502</td>\n",
" <td>0.000554</td>\n",
" <td>0.000602</td>\n",
" <td>0.000647</td>\n",
" <td>0.000691</td>\n",
" <td>0.000731</td>\n",
" <td>0.000766</td>\n",
" <td>0.000802</td>\n",
" <td>...</td>\n",
" <td>-0.000048</td>\n",
" <td>-0.000035</td>\n",
" <td>-0.000012</td>\n",
" <td>0.000020</td>\n",
" <td>0.000063</td>\n",
" <td>0.000117</td>\n",
" <td>0.000202</td>\n",
" <td>0.000334</td>\n",
" <td>0.000493</td>\n",
" <td>0.000613</td>\n",
" </tr>\n",
" <tr>\n",
" <th>SN112EK62.5</th>\n",
" <td>0.000400</td>\n",
" <td>0.000472</td>\n",
" <td>0.000539</td>\n",
" <td>0.000602</td>\n",
" <td>0.000659</td>\n",
" <td>0.000715</td>\n",
" <td>0.000768</td>\n",
" <td>0.000817</td>\n",
" <td>0.000861</td>\n",
" <td>0.000905</td>\n",
" <td>...</td>\n",
" <td>-0.000052</td>\n",
" <td>-0.000037</td>\n",
" <td>-0.000011</td>\n",
" <td>0.000025</td>\n",
" <td>0.000073</td>\n",
" <td>0.000133</td>\n",
" <td>0.000227</td>\n",
" <td>0.000374</td>\n",
" <td>0.000552</td>\n",
" <td>0.000684</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>5 rows × 33 columns</p>\n",
"</div>"
],
"text/plain": [
" SN112EK42.5 SN112EK47.5 SN112EK52.5 SN112EK57.5 SN112EK62.5 \\\n",
"SN112EK42.5 0.000313 0.000341 0.000364 0.000384 0.000400 \n",
"SN112EK47.5 0.000341 0.000381 0.000415 0.000446 0.000472 \n",
"SN112EK52.5 0.000364 0.000415 0.000461 0.000502 0.000539 \n",
"SN112EK57.5 0.000384 0.000446 0.000502 0.000554 0.000602 \n",
"SN112EK62.5 0.000400 0.000472 0.000539 0.000602 0.000659 \n",
"\n",
" SN112EK67.5 SN112EK72.5 SN112EK77.5 SN112EK82.5 SN112EK87.5 \\\n",
"SN112EK42.5 0.000414 0.000428 0.000439 0.000448 0.000458 \n",
"SN112EK47.5 0.000497 0.000520 0.000541 0.000559 0.000578 \n",
"SN112EK52.5 0.000575 0.000608 0.000639 0.000666 0.000693 \n",
"SN112EK57.5 0.000647 0.000691 0.000731 0.000766 0.000802 \n",
"SN112EK62.5 0.000715 0.000768 0.000817 0.000861 0.000905 \n",
"\n",
" ... DREK27.5 DREK32.5 DREK37.5 DREK42.5 DREK47.5 DREK52.5 \\\n",
"SN112EK42.5 ... -0.000031 -0.000027 -0.000015 0.000005 0.000031 0.000063 \n",
"SN112EK47.5 ... -0.000037 -0.000030 -0.000014 0.000010 0.000042 0.000082 \n",
"SN112EK52.5 ... -0.000042 -0.000033 -0.000013 0.000015 0.000053 0.000100 \n",
"SN112EK57.5 ... -0.000048 -0.000035 -0.000012 0.000020 0.000063 0.000117 \n",
"SN112EK62.5 ... -0.000052 -0.000037 -0.000011 0.000025 0.000073 0.000133 \n",
"\n",
" DREK60.0 DREK70.0 DREK80.0 DREK87.5 \n",
"SN112EK42.5 0.000118 0.000197 0.000295 0.000373 \n",
"SN112EK47.5 0.000147 0.000245 0.000364 0.000458 \n",
"SN112EK52.5 0.000175 0.000291 0.000431 0.000538 \n",
"SN112EK57.5 0.000202 0.000334 0.000493 0.000613 \n",
"SN112EK62.5 0.000227 0.000374 0.000552 0.000684 \n",
"\n",
"[5 rows x 33 columns]"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#obtaiting the covariance from y_pred, however this is incorrect\n",
"Ycov = Ydf.cov()\n",
"Ycovdf = pd.DataFrame(data= Ycov)\n",
"Ycovdf.head()"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
" C:\\Users\\danny\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:8: MatplotlibDeprecationWarning:Support for passing a (n, 1)-shaped error array to errorbar() is deprecated since Matplotlib 3.1 and will be removed in 3.3; pass a 1D array instead.\n"
]
},
{
"ename": "ValueError",
"evalue": "operands could not be broadcast together with shapes (33,47) (33,) ",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-13-a4695423b236>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 7\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfeatures\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mYy\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'r'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 8\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0merrorbar\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfeatures\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mY_ex\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0myerr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0my_ex_std\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 9\u001b[1;33m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfill\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfeatures\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mYy\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mYcov\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mYy\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mYcov\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0malpha\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0.5\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'k'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 10\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'$E_{c.m.}$ (MeV)'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 11\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mylabel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'$R_{n/p}$ & $DR_{n/p}$'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mValueError\u001b[0m: operands could not be broadcast together with shapes (33,47) (33,) "
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd1QVVxfF99CbgAgIIiCIvaCCNXaNLcbeNWrsMUaNxhiNnxqNJZoYTey9d03svcaCBo0oVgSlidIUC0XK/f7YPB8ISK/e31qzHvDmzcyDx753zj1nH0UIAYlEIpF8HGjk9wVIJBKJJO+Qoi+RSCQfEVL0JRKJ5CNCir5EIpF8REjRl0gkko8Irfw6sbm5uShTpkx+nV4ikUgKJdeuXQsVQlhk9fX5JvplypSBu7t7fp1eIpFICiWKovhm5/UyvCORSCQfEVL0JRKJ5CNCir5EIpF8REjRl0gkko8IKfoSiUTyESFFXyKRSD4ipOhLJBLJR4QUfYlEIskrhABmzAA8PPLtEvKtOEsikUg+KoQAJk0CfvkFiIwEnJ3z5TLkTF8ikUhyGyGAKVMo+CNGALNn59ulSNGXSCSS3GbaNAr90KHAkiWARv5JrxR9iUQiyU1++gmYORMYPBhYvhyIj8/Xy5GiL5FIJLnFzJnA9OnAwIHAypVARATQoAGweHG+XVK6oq8oylpFUYIVRfH8wD5NFUW5oSjKbUVRzuXsJUokEkkhZPZsYOpUoH9/YPVqCn6rVsDNm0A+2spnZKa/HkCbtJ5UFMUUwFIAHYQQVQB0z5lLk0gkkkLK3LnAjz8C/foBa9cCL1+qBb9DB8DKKt8uLV3RF0KcBxD+gV36ANgrhPBL3D84h65NIpFICh/z5zM1s08fYP164NUrCr6HB9ClC7B7N7B/f75dXk7E9MsDKK4oyllFUa4pitI/rR0VRRmmKIq7oijuISEhOXBqiUQiKUAsWAB8/z3QqxewYQMF/9NPKfgDBwLbt3NB96ef8u0Sc0L0tQC4APgMQGsA/1MUpXxqOwohVgohXIUQrhYWWe72JZFIJAWPDRuA8eOB7t2BTZuA16/Vgj96NLBqFdCpEzN4FCXfLjMnRD8AwFEhxBshRCiA8wDyp9RMIpFI8oMnT4AxY4AmTYAtW5IL/uTJwKJFQNOmwLZtgFb+GiHkhOjvA9BIURQtRVEMANQFcDcHjiuRSCQFHyGAr78GYmKYpfPmjTqG//PPwLx5QPXqwL59gJ5efl9t+t47iqJsA9AUgLmiKAEApgHQBgAhxHIhxF1FUY4CuAkgAcBqIUSa6Z0SiURSpNizB/j7b1osmJtT8G/cAH77jSmbpUsDR44Axsb5faUAAEUIkS8ndnV1Fe7u7vlybolEIskRnj8HKlUCbGyAY8eAdu0o+EuW0HpBUYCLF3M0L19RlGtCCNesvl66bEokEklW+e47IDQUOHyYvjo3bjDEM3MmEB0NnD+fr4VYqSFFXyKRSLLCqVMsvJo4kWL/998U+99/58LuyZNA1ar5fZUpkKIvkUgkmSUyEhg2DChXjjYLdesCjRsDJ04Anp7AgQNA/fqpvzY2lhk8+ZS2KQ3XJBKJJLNMmwb4+DDnfuhQQFOTC7b//ANs3Ai0ScO55tEjoGFDmq/lE1L0JRKJJDO4u7PydtgwwM0NuHQJ+OorYOtWYNw4oHfv1F+3Zw9QsyZw/z5gaZm315wEKfoSSX5y7x7w5Zcs7JEUfGJjaaNQsiS9daZNo5/Oli1A+fKM6b9PdDQwahTQrRtQoQLw339A5855f+2JyJi+RJIf3LoFzJoF7NzJ4h4DA2Dhwnwtz5dkgPnz6ZS5fTvbHpYsCZiYAAEBwIULgL5+8v29vICePSn048fTbvn1a9osm5jky1uQM32JJC+5fp2zvOrVgUOHmPkxeTIXBsM/ZGYryXfu3wdmzOCM/cIF3qWNGQOsWwd8+y2boyRl2zagVi3A15cLu1OnAnPmAA4OtF7OJ6ToSyR5gZsb8NlngIsLcOYMBcDXlyJQqxb38ffP32uUpE1CAhdsDQwYzlm8GBg5kkVY5cvTbkFFZCT37dMHcHYGLl8G7tyh2E+fDrRoAfTtm29vRYZ3JJLc5Px5xnlPngRKlKA4jBqV/Nbezo6P/v5AjRr5c52SD7NyJTNz/viDYZrKlYG4OMDPjz9XhXXu3AF69ODj998D1tZM5Xz2jBk9M2cCrlkups0RpOhLJDmNEMDZs/RMP3eOcd/58xkDNjJKub+tLR/9/PL0MiUZ5MkThuFatODfNTSUf9thwxjW+eQT7rdvH2f3BgYM+2zbxoG8SRM2TmnYkPvdv8+tQ4d8eTtS9CWSnEIIzuhnzGDM19qa1ZnDhlEI0sLSEtDRkeGdgsq339JBs1Uriv9PPzEs5+SkDuucOUMffVtbfg4WLgTq1AHWrAFatlQv0CddqM8n3zMp+hJJdhECOHqUYu/mxiKdxYuZ2pcRK10NDb5GzvQLHkePMsNq3DgKfKNGDNU8fsy7OAMDZuZ06ADo6rJgS2Wj/PnnapHX0Mg3kX8fuZArkWQVIYCDB1mC364dwwDLlgEPH9JfPTPe6ba2cqZf0IiK4t+xQgXgyhX+bORIYOlSdsJq1Ih/61atgLdvuYC7aJF6EFAUFmopSnLBNzLiwnA+IUVfIsksQjDd0tWVs7mQELbC8/Ji3F5XN/PHtLOTM/2CxuzZnLk3b0575F9/BX78EShbljUWQUF8Ljwc0NDA/GGz0VO3Nmf1AMV++/bkx3zyhOGiM2fy/v0kIsM7EklmuH6ddrpnzgCOjnRZ7NcP0NbO3nFtbYHAQCA+nj4ukvzl7l02RenQgXn4n39OIzUfH4Z1YmOZlePvz0yso0fh7pE4m0+twO7WLWDzZqZ3vn7NfZo3z9v3lIgUfYkkI/j6cpa3ZQu7I/35JzB8ePbFXoWdHQU/KIjxfUn+IQS9dIyMWGlrYECrjC5dgG++AWrX5vbwIRfrz54FypfH9PY2qBj2JPmxmjRhHUadOrRj6NmTn6P8tFwWQuTL5uLiIiSSAs/z50J8/70QurpC6OkJ8cMPQrx4kfPnOXRICECIS5dy/tiSzLFhA/8WXbrwcd06IZychHB05N/e2Zk/d3QU4tkzvgYQCRwu1NvIkfzcaGoK0b+/EPfu5cjlAXAX2dBeGdOXSFLj7Vsuyjk5Mce+Z0/gwQOm6uWGZ0rSAi1J/hEezuIrZ2daJ3TtygbnDx9y3aZJE35fpQo9eCwt34VzFLBJOAYN4h3gqlXAF18wJ3/DBi4IFwBkeEciSYoQwN69zMf29mZBzvz5tMTNTWSBVsHghx/Y99bCgoP7l18ynj9sGDN2bt/mZ+Hq1RSNUBIACA1NhgCHDuVnSDWYFyCk6EskKu7fZ4reqVOMuR4+zNL5vHC+NDEBjI3lTD8/uXiRs/NPPuHXW7dy0d7GBjh9mrN9V1emb2popPhcaAA4W7c1mu5eBZQqlT/vIQPI8I5EEhkJTJkCVKvGBhmLF7Pnadu2eWt1bGsrZ/r5RWws022trCjqvXoxfHPvHn11Hj7kgqybW6qCD319LBkwBcsGTCnQgg9I0Zd87Bw6xPjsrFn8R1fN9vMjbVIWaOUfCxcyJVNPj8Z4gwcD8+Yxk8rLC7C3p7XGoUMpBb9cOcDNDefrt8vQqYbOP4Rh8w/lwpvIGFL0JR8nfn70tW/fnjO5M2fY27Rkyfy7JlmglT/4+tLyuHx52issWcICKl1dpmwaG7MlYuvWQMeOyV/brRvvDqtX//A5wsPpw9OqFZb/0BHtTu3IrXeTLjKmL/m4ePuWJmgzZvD7X34Bxo6l4VluosrB9/NLfatZk5WeISHM586MhYMke4weTVuEhw+B/v1ZSOXpyee0tIDTp/G8dBmYili8m+NrarJCd8yYtEOAL1/Sg2fHDuD4cYaQHB2xv1VfXKjTCp3y4r2lghR9ycfDqVP8B79zB+jUibf09vY5c+yYGIq3ry9ni76+ybfAQPqvJ8XUlLP7N2+Y8bFiBX8eEMBUUUnus2sXsH8/7/A0NZl107ix+vnNmwFXV5gCasEvVYqve79TFgDdmCgatG3fzkSAmBiG7caMYdqviwu2r3TLi3eWJlL0JUWf69eZinfiBFCmDPOv27fP+OsTEoDgYMbbAwK4+furRd7Xl7P4pKicM+3tacxlZ5d8s7Vl2AAA1q9naqBW4r+jn58U/bzgxg3+3q2t+ffbv59++CpztClTUna4atmSA7SlZbIfx798BfOLZ/A6Og4xe2dC18KcaZ69egH16qn9eAoAUvQlRZeHD/mPu2MHYGYGLFjA8vqkoZPXr/kPHxREMyzV1yphDwjgLD02NvmxdXQo6mXKMMvH3p5bmTJ8tLHJuEVDxYp8jIrio1zMzX2Cgph/b2gIPH3KhdstW9S/+y5dGL6JjwfAWf7etgPQ5cCa5Iv8QiB4x1/49ugjXCzlDADYs+4w+vRpXmA9lKToS4oWcXFMs/v5Z3Yr0tKiKLu4cAY9aFBygX/1KuUxdHUp6KVLM2fb1lb9veprC4ucS+dUVWqqGqPLxdzcJTKSC7Lh4ZwA2Nuz/eH48Xze2ZmhwOjody9Z1WcCTjbujC5JhdzbGxcm/YKxlo3x2sIR9RPCcVPfAsuDtNADSoEV14J6XZKPESEYA42I4CJY0kfV9uJF8u+T/jwkhI9JiY8HjhzhVqwYb8utrfmP3bYtv7a2ZpxW9XXx4rmfn//mDTNCzp3jpqfHOxNLSznTz00SEoABA5hxU7o0P19Tp3KmDzBdMzCQnykVW7bg5CsH9ffR0Yj7ZR4WnfXB4jpdUVYrFltGNsXUg/dQ6s1beAW/xsGbQehU0yZv31sGkaIvyThCAJcvc3asKNxUhSrvfx0by3+otLZXr1IX9/fDKKlhZMQK1qRbbKxaLGvVYhZG1ap0xLSw4D9zVnzuc4qXL1nlqRJ5d3felWhqMvQUHU07X5m2mbtMm8Y7wHLlgEePuJj/1Vf8bGtp8TPyJIlT5v79DAOtuMzvjx3D0+9+xOhq3XC1bnd0q2SGGb1rw0CHUlrcQBsVShbDkjMP0cG5FDQ08rC4L4NI0ZdknBs31E2gs4K2ttpuwNiYM28bG95aGxurn0vtUbUZG6sXPENCWD27ZAkQFgZ8+ikN0Vxccub9ZocXL4B//gHOn6f17vXrnGVqa9OWd8IEmnc1aEBh6dePWUXNmrFATJLzbN7MsF/58jTPGz2aNgsxMXze2jr5XdaZM0DTpgAAs+fBGLBrEc4+j8C4jt8jSt8Iv3V1RleX5DbYiqJgZLOyGLP9Bk7cfYbWVazy6M1lHCn6koyjmoGuW8eZkspENiEh5dcqgS9WTC3yOTXT9vbmouzatZwhd+hAEW3YMGeOnxVCQynw589zJu/hwd+Djg7bKU6eTAGpXz9lk/TKlfn44gXvSk6c4Gvz0gKiqHPxIkM4jo4U/DZtgMWLEaeKvdvYJBf8q1c5OAPAwYOYO7M/FtfthlWfdkHFkkZY3NcFTpZGqZ7qs2rWWHDiAZaceYhWlUtCKWB/Ryn6kowTHMzHFi3UrpB5ibs7HS9VC7RffMHFt0qV8vY6hGDc9+JFtcjfvs3n9PUp7NOncyZfpw5/9iEqVlT3UdXWZkZRRATz+CXZx8eHdRlmZvy6YkU2PDcwgFZkJPwUPdgFBqr3v3WLoUEhgDlz8HDBcozuNRt3LMqgdx1bTPu8CvS0087M0dLUwIgmZTFp7y384xWKxuUt8uBNZhwp+pKMExLCR4s8/BALARw7Rh+UM2d4xzBhAm/N88rY6vlz4N9/uV29ykdVXr6REUNeffuyqKd27cxX9+rrM5bv6/suRRD+/lL0c4KICMbkY2J4J2VuzuwuMzMgPByPFUPYizfq/R8+ZGX0mzeI+XIwlj3RwNLBi5GgoYmy5oaY0yUdu4VEutSywaKTXlhy5mHhE31FUdYCaA8gWAiRZo8vRVFqA3AD0FMIsTvnLlFSYAgORpSuAQZu+A87htfP3XO9fcuqxl9/5czLxoaz/GHD1EVNucHLlzyfuzsF/upVCoGKChV4p1OnDsM2tWqp1xiyQ/XqDJ+9fMnv/fzo+inJOnFx6uY3enr8O4WGMkXT1xcAYC/eqCttAwL4OXv0CFe/HItJ5drC29EWHZytEfA8CtqaGS+w0tXSxLDGjphx8A7cH4fDtYxZzr+/LJKRT+t6AIsBbExrB0VRNAH8AuBYzlyWpEASEoKIYrk8+4yIAFauZNeqwEDeZq9bx0rJnPTHiY7mjM/Tk5vKbyVp5oyNDcV90CA+urjk3uy7alVWCj99yu9l2mb2EILWB8eOcVb/4gVFv2pVta8OWHQlACghIYC5OSKOn8bcJYexrd4w2OgB63rXRrMKluipyt7JBL3q2GLxmYdYcuYh1n1ZJ+feWzZJV/SFEOcVRSmTzm7fANgDoHYOXJOkoBISglfFiufOsf39mT63ahXTOVu0AFavprNhdhbCwsM50/Py4nb3Lv/pvbzUoRRtba4LNGxIUahalQKfl77oqsXchw8pTjJtM+tERwNDhrDC1sKCYUkTE/6N3ZL73iQA6P/nGWwyM8OReWsxzU8bYRWbYGjV4vi2R513qZhZwUBHC4MbOmD+sfvwDIxAVZtcaLOZBbJ9X6ooig2AzgCaIx3RVxRlGIBhAGBXANuIST5MRPhLTG4yGG9j43PuoP/9B/z2G60ShODt+PjxDJtkBCGYrunrqxb2pCKvqnIFWEfg4MCwSbduFPdq1ZiJlFHLhNyiShU+BgSwaEjO9LNGcDAtsy9dUtc8WFpyUD17NsXuvZddRHx0DIaOX4OTuqVQRSsUawfVQrUKOVNY1a+ePZaf9cays95Y0jeDn+lcJicWchcCmCiEiE8vNUkIsRLASgBwdXUVOXBuSR5yE0b4t2R5mIVHZu9AqsXZX39lubuREfDNN7Q4fn8yEBXFYpmkNsS+vsm/V3nWqChdmrnY3btT0MuV4/cODvlboPUhVFYMCQmcncqZfua5fZtGekFBDMP5+SG0uCXMa9SgtfH7CIHIGQfwMCIWQsMMkzV9MWj+MGhp51x+i4m+Nvo3sMfSs954GPw6zTTPvCQn3p0rgO2Jgm8OoJ2iKHFCiL9z4NiSgoIQCH/LcTo8MhY3fcNR3T4Ti1Px8VxE27iRFsLe3qyS7dWLi5ivXtHjPjiYt+PBwdxev055LCsrDg7VqgGffaZ2rnRyYubF+3nwhQEDAw5WAQH8Ws7000QVX0+WTHDsGNCjBycUsbFATAyuV6kPiASYvy/4dnaI9fbBrF9241akPuo89cJvn5WDbZ+R2bqutJIbBn3igDUXHmHZWW/81sM5W+fICbIt+kKId6YUiqKsB3BQCn4RJCICobqcpRjFRGL+90uwae8MZkXo6/Mx6ZaQQCFX2S5EpnJ3EBbGDJ3t2xnHtrDgrbilJcXb0pI/s7JixoXKkrigztazS7VqFH1F4WNCQoGy5C2wLFnCO0V9fX7OtLWBP/6A5pKNcL57Nfm+27cj3LU+vh67BpeNbNDJ6yLiq9eAbZ8uuXZ5JYx00buOHTZe9sXYluVy7TwZJSMpm9sANAVgrihKAIBpALQBQAixPFevTlJwCA5GuL4JNBPi0erZbex1qI1L3/2MBjHPuHAWFcVH1aYobEzh48MZPsCimLZtmQmjqtJVCb2pqaxAdXamMVxUFGerz57RGkCSOnFxbGu4eDEHx8hIWltPnQr88AOqq4oJVcTH4/aqrRi24DRC9C3wm8kz7Gr5GUQeDKzDGjtis5svVp73yfVzpUdGsnd6Z/RgQoiB2boaScElJAThBiYwi4zA7JkD4LbTB78YN8TfIxskLzN/+xY4fRpYtowpiNradDX89lv1YqUkdVQZPCqx8veXop8G+lGvgVatWLAH8K6oe3feZQ4aBCBJOiYABAZi/xfj8L1VIxTXeYvdXZxQvUE17MxCKmZWsDbRR9dapbHD3R9VrI2ho5V/d3Dy3lGSMUJCEGpgghKREdCzscbYluXh4f8Cx24/o03w3r00DbO05Gz+4kU2MPHzY+qlFPz0Uf2Onj3jo1zMTZWoJ0/x7NZ9zEmwx5XSVRCnrUPjtEuXgE2b3u0nAMQCiF+/AXMGzcBo209RzUBg//SOqN4gc4VvO4bXz3ZB4ogmZREXn4CnL6PT3zkXkTYMkowRHMyZftRLoH17dGnUGCs0quDXNSfx6coR0Ix8wyKYzp3ZdejTT2Vz78yi6qClat4hF3NTcHftDtx9owtzXSOsrd0RK+p1g3F8DBo/uILmpuXRJPwVSkSxqnlGr2/R9P5NrDn0GOdrfI6+FU0wrV/bfJtllzE3RPvqpXDw5hOUMsm//w0p+pKMERKCcAMLVH36ELh4EVqnTuG78g3wVefJ2PtJZ3SvYw+MGsVFV8k7Us00SQsDA/7+nj7lgCln+mri4xHduy/GFGsME31jHFw/BjqVKuJitB5O29fAGUdXHKzUGIpIQI2gB2he2hB1/tqCqY0H4ElxK8zuWAV96pfJ73eBr5s5Yb/HEzx9GZNv1yBFX5IxgoMRauCEEpERXHxt3x5tzEqg+ptQLCzXEh3mDoHu7NmsZG3dGmjenOZjxYrl95UXLqpUoeibmsqZvgpPT6BFC8yt1gkPHO2xfudUvDC3QtnrV9EGQJs755EABbdLOuL0wHE4ramFBfE2EB1/gC4SsG1EgwLjfVPBqhhKF9eHsV7+FQNK0ZdkiLchYXhlZ8Twzu7dQJ06UAB87xWKfmuuYMvGExjk/Q/zpefOBWbNYjZO5cpAvXo0J6tbl6JWQBtGFwhcXFiwpqMjZ/rx8bSonjULZxxcsN61A7787yDqPbqO94MjGs7VUe3hQ1SbNwpjundHaPem6H3+BYx0tT4o+LluHJgKNqbpWG3nMlL0JRki/PkrwA4o8eaFOvYMoGE5czQoWwJLvF+hx/eTYPS//9Hcys0NuHKF219/AWvW8AWGhoCrq3oQcHFh/v3Hnq6pQrWYGxPzcc/0794FunYF7t5FiJEZJrQbg4phfph4ak0KwYe+PnDzJgv9pkwBKleGOQCza3mTmVPYkKIvyRBhiRkHJaIiUlgbf9+mIjotuYg1/zzCmJblGJpo04YbwCpJb2/1IHDlCvD77+p+uCVK0Gsn6Va27Mc5EKhE/9Ur5p3HxBTdYrTUiI1l74Rp04D4eAhdXXzf5hu81DXElu1ToBefpIeylhZTNbt2BX78MdlkRJI2UvQlGSI8kv9sZjFvUjxXw9YUrauUxKp/fPBFfXuYGb5ngawotEhwcmKzEYBiduMGe8devw5cu8YWiKqBwMQEqFmTA0DNmqxWrVix6AugSrhUFcyBgWzx9zHw33/Al1+y1SQAaGtjY5VPcaZsbfx0YjkqhPqq99XQ4Gfpxx/prSTJMFL0JekjBMLiOOsu8SqMIn7hAitrE90pv2tVASfunMfSMw9xKzACQDrxUl1ddYhHRUwMTbNUg8D16yyxVzWu1tKiMVm1avTrqVaNW1EKDxkastYhaYFWURf96Ghg5kzgl1+S/R3vm5TCrGaD0Mz7X/S/flC9f+fOvBtwcsqHiy38SNGXpM+LFwjTo+9OiUgK+rsm5B07Ap9+inKffooutWyw0c0Xla2LQVcrC4u1urrq8M6QIfxZbCytkm/d4nbzJnD5Mv16VJiY0Ca5UiUOChUr8tHBIWe6WuU1FSqoRb+oL+ZevswK2nv3OIFI7HEQramNMZ9/B+OYN5h3eBGraitV4t+9esZaFkpSpxD+R0jynOBghBmYQis+DsbR74V39u3jBmBsxRrY3+EnxHp5w9LGPGfOra3NOHeVKlyoUxERoe54dfMmv963T93HV/VaJyeKqGowKFeOM2crq4J7d+DqCvzzD78uqou5b94wNPPHH7TWBoDY2He2Cb80HYh7lg5Yt2s6LOKj6Mw6bFh+XnGRQYq+JH0SfXeKR72EBtJug1D63g30tTmIDbXa44+5w4FFUM/cXVz4WKpUzoitiQkbkn/ySfKfh4cD9+9zu3dP/fWhQ+r1AoDFTw4O6s3RMfnXudmHNz2cE+13i2ra5qlTwNChwKNHNOVT2U6Agn/WoRbWuXbEgGsH0MxKB7gSwIbmkhxBir4kfYKDEaZvrA7tVKkCNG7M1oZxccl2/fryTuys9ilmNh2M9QY+jM0fPMgMHoDx6qRZOjVrUmhzatZtZgbUr88tKXFxFJmHD/no46N+vHiRdw5JMTFR+/Tb2qb82sYm97ptJfUpKiIz/Z4rLkM/6jXW39nFz02pUpzhJxF8AAjTN8Z37b5F+VBfTGpXCRi7LMufjfzIwS8MSNGXpE9ICMIMTdWif+EC0zLHjWNe9I4d/McUAuaRERjptgvzmwzAqd1H0ML7LhuheHgkX6A9cULdo9bEBKhRQz0I1KzJUExOxuO1tNRdtFLj+fPkA4G/v7ozl5sbvf+ToigMEdnYUMBKlaIjpurrxE1JSMi8dW+lSnyMjS0yol/r5kUM2ToPeBnGxfsrV5I9H2xoiptW5bGmdke81DPCps5O0GtRL5+utmgjRV+SPsHBCNe3QrUIL35vktjg2cmJC2sTJgCTJ7MlnY0NBl/9C/srN8GUViNRd81IGKlittbWbH0IMGPj1i2m6am25cvVrQ/19JiZU6MGB4EaNfi9US61mytenCEoF5fUn3/zhgKcdDDw92drvseP6fCo6huQhLG2VeHwOgRYUYqhjLS2pB2/DA1Zu6Dq/VuYCQsDxo7FxM2bEWhlD1gWx8v/buGWvTM8rMrBw7o8blqXQ5CxBQBAQySgrIkOKknBzzWk6EvSJyQEYQbl1TP992+3XVxov3D6NPDDD9ALDMS4qzsxot13+LXRF5h+aiX3CwpSv/bePXrz1K6tPk5cHDN1kg4Eu3czHKA6b7lyjHnXqMHN2Tnn1gk+hKEh7z4+VAD09i19c548AZ48gfujMPQKKYUePpcx781/DGV4evIx6foCwB6+d++qO2WVK0fBVHUey881hqyyZw8wciQQHo6NTfvA3cganlbl4NO59Ltd7N+EwjXgHpzvnIDzoO745Y0lNDUK6AJ7EUGKvtPEehgAACAASURBVCRdYkLC8MreSC36adG8OXDlCn4bMQfdDq1D/+uHsMGlPTpUsUCtP2Yl31clnkZGrD4FGIKpXJmbqohLCM6ob9zg5uHBENGuXepjmZtT/CtUUIdwnJy4VqDzXqFYbqKjo479A9i0/T8g5An2V/gEv8yYqW42IwStKp4943bsGDBnDsNmjRtzn1q1GFYC+P4LUz+CZ8/ouLp7N1CzJk6Uq4dZtbrCIDYargF30OX2aVQP8kL1ctYwPXGY/Re2rQGsraGZR01NPmak6EvS5Xn4K8AeMEtP9AFAUXC1VjNcrdkUa6zCcfzsS0x6bY0DpW2hM2Y0m5+rRB5gvF8lhiNGsOPWe8d7J6QdOqh/HhHBVE0PD/VgsGVL8gVZDQ321lUNAqqBQBV7t7LKtcXY0NcxOHwrCDpaGoiOTcAN/xeoaVdc/Z6KF+dWsSJTNP/4gw1AVKKf9A6osIi+EMDWrcDo0QyH9emDHZ4hmNRgMKo99caa3dNhHvWS6ZeTTgIX7rH47quvCm76bBFEir4kXUIjaAmQ7kw/kaRZEzPLP8WQjdewqnEffD1hAmflEydysbJr1+QvXL6cG8DwT4UKaZ/ExARo1IibCiEYEvHyYpZO0sf3BwQV5ubqQcDaWv11mTKMszs4ZKkZzI5//REbL1ClVDHcDXqJXdcC1KL/PoaGbDyzaxfw5588X1KRLwxpm4GBHLQPHgRcXCCuXcOSR3H4te0YNPa5hqV/z8GRyvXRvYYdMHw4B7rNmz/8N5bkClL0Jemi8t3JqOgnpWVlK3xWzRqLtJqg7fH2cFw0lxk/xsZc/P36a+ZsHz6c/IVJY+ci7dqAZCgKRdzcPGXKpmpA8PXl2kLS7ckTPt6+zZj8e2mosLHhAODomPyxXDmmiL5HfILA1it++MSpBOLiBYob6OCAxxNMbV8ZetppVCp/8QVn+gcOsNerKoMHKNgZPEIAa9cykys2FqhbF/FX/8WMlsOxweVzdPY8jV+O/IE5Yxbiq42zgFvngP/9j1sqd1kyzTL3kaIv+TAJCQhP4OKiWVTmRR8ApnWojPNeIZj8SBPbDhyA4uEBzJ7NOPaCBWycfv8+M1ZKllSncqpQ3foXK8ZFzayQdED4EAkJrOp99IjOoD4+6sfjx9XZRyqsrWkBkWQ7o2WJwBdRmPJZJay/9BgWxXRx7+krHLv9FB1r2KR+3ubNeazNmyn6RkZMi33xgucviPj6csA+cYJrKh4eiHG/hnEdJuBQxUYYenUvJm2YDo2ZrzDt91EcKA9cSDkgS/IUKfqSD/P8OUL1maJp/uZFlg5hWUwPk9tVwqS9t7DT3R89a9cAdu5kps6CBcD69cDKlfTxOXeOVbY3b6orU1W8eqUeAOzsciedUUNDnUZZL5W0wago9YBw/z6zcTw9k6Wbbuo+HSWtnNDyx68QHWuCACt7vHKoj91XHqct+pqaQJ8+wKJFTP00N+cdxfXrDE8VJBISGJf//nv1XZiHB17p6GN45ym4VMYZPz69hKHftGdO/rNnbFz+00/qtFRJvpE/HYIlhYdECwat+DgYp2KrnFF6utqijoMZZh26i+BX9ObveSYMPWsOYMx6yhTg/HkauTVoQFGNi6PD5tKlKQ/o58cBQFGSL3rmNvr6zC76/HMK2fr1gLs7ByQvL/hu+wvnHV3QO8YX2j4P8fnxLfhm/Ux0PbEZF3zC8cSpCmf1X31FgT92jINXQgJDPHFxLHYD1INeQQrv+PgALVowFVNPjwu2YHFVzz5zcdW2Cn5vao2hr+9zncLSErh6FZg/Xwp+AUGKvuTDJPHdyU5+hYaGgjldqiE6LgE/HbiT/ElLS2b1+PlxIfPpUwpGpUqMFw8cyBnl69cU2vdxd1cPANWqZeMqs4GmJuDkhK3GFaChoYHe88cBnp7o/+cZjJ2+Dd1G9YBQNPBX6/4sTNu+HRg7lo1mypTh3U2VKrz+zZt5zAYN+BgamvF1jdwiIYEZRtWqUcRV1wXgkak1uvb7FY+tHbHa9iU6927BBd05c4B//0274E2SL0jRl3yY4GCEGpioF3Gz0cSkrIURvmnmhEM3g3Dq7rOUOxgaMr/by4vhH1NTzohtbbnoGx7OGaMQFJzWrVMew9NTPQDkcUFTdGw8drj7o3WVkihpzIyfeE0tBFnZw65nR9R1MMOuMnUgLl7ke3n2jHc3//sfc/I3b+Zs382NvwOVhXB8fHL30LzGywto0gQYMwYJkVHw1iuOfZWaYFazQejdaxbafvkngk0tsfXWVjQd05/XffMm8MMPuedPJMkyMqYv+TCJM/13op9oo5xVhjcpi4M3g/C/vz1RylQ/9epLTU0uZnbrRlFctIgNNubN4x3AmDGcBR89yv39/bmvagaqIukagIZGygXiHObQzSC8iIxFv3r2qT7fzaU0Juy+iWu+z9ms29KSW8OGwJEjjHmfOsWU1s2bk9/V+Plx3zyk59J/UPvCYTjc/Q+3SlTG7T6f446lA97oMkyjE/cWlYrroLq/L37ePA/lY54z1j9kiLqyWFLgkH8ZyYcJDkaYgQlKRCYu4qqap2QRHS0NzOlaDUEvo+H/PPLDOysKZ5h79zLG/+23zBRp2JB53hs2MOZva0sDLyEYH1c1YElKQgKEoiBBdReQC2xy80VZC0PUdyyR6vPtqlnDQEcTu9wDkv2850o3zK7fhx4+R48y5r95MzN4VHHwPPbg8T/rhtfXPLDYriHGt/4GO6q3glCA7rdOYv6h33HU/2/cbqyNfZu/w84VX6P8JzVpIzFsmBT8Ao7860g+TOJM3ywyMVUyBxbjatkVR/969nj2MgavY+LSfwHAuPf8+UBAADNloqMZ67ezA6ZOVadS2tnRq0cI7jtp0rtDKEjygVeJv6IwhJRJeq64jJ5JLAM8AyNww/8F+tWzV9stvIehrhbaVbPGoVtBiHyb/H17VK7LQrOffwZ69OCC6eXLrCgGWHWcB4jwcOwd/TPa7vOHn7ElZhxfhpOrR8BzYQ/s3jIR00+tRPepw1HxTQi0O37OheeDB4G//mJRm6TAI0Vf8kFigp7hla5h2mZrWWRCm4rQ1lTwKPQN4hMysUhpaMiKTk9Pzvrr1qVQ2ttTLM+cUS962tiwHkAIICgIa3qNR0Jqx1y+XD0AZDEGvdnNF/ramuhSq/QH9+vuUhqvY+Jw1PNp8icUhe8jKIjxe319FmupMng8PbN0XRlGCESs24TRwxZgnEFNVA72wZG1o9D/v0NwCguApkhg57JJk7jucPIkMHcur+uzz3L32ooYO4bXz9ciNCn6kg8S/pQ+8u/COzmEka4W7M0MEPk2HluuZCF0oShAy5bA/v2AlxcONumK14cSQyOVK3Md4EWSa7aywvGmXdF7+SX+fPt2Cuv7xMUlvwtYsiTdS4mIisXfNwLRqWYpmOgnHzTe/wev42AGOzMD7L4W8P5h6LvTqhXw++8U0h071Omo9+5l6NeSJTw94daqO9r+G4cjDrUx4dwGbNs2GaVfJi4em5oC06fTEG7OHK63PHjAtYdsLOxL8gcp+pIPEpbou/MuvJODmBnqwFhPC/OP3UfIq5isH6hsWWzqPgYj5u5n3ryJCdMhS5UCBg9mSmdSTEyAnj2ByEjaIZ88mfZsddSo5INAKuy5FoDo2IQ0F3CToigKurmUxiXvMPiHp7Km8fPPtIvQ02NjF5UlRFBQusfONK9f4+3IrzHv6/noXWsAdONisWfzBHzttoszewBo146VxtOnAxYWFP5Nm2QopxAjRV/yQcITbd/fhXeSOmRmE0VRUMbcEDGxCZhz+G62jxero0tLBzc3VrL268cZfaJvf9OLB6HzNjr5i7S1WWx08CBz0T08gJkzP3TRgKJg64gGGLJgHIQQ2Ozmi5p2pqhSyiRD19mllg0UBdh7PTDlk7VrA5068Q7GwoKtHAHWKOQUQgCbN8O7TCV0fVMOS+v3QC+P4zi4fgycnyZW/zo6Ap07c2H53j1WTP/7b8qexJJChxR9SdokJCBMk/nm78I7pqa0OPbxyZFT6GtrYngTR+z9LxBuPmHpvyCj1KxJoXryhAVfkZH4atNsLJ/YgSZv16+nfI2iMMd8yhQKY0gI3TnbtUuxqyaAlg/ccMmhJnxC3+CLpf/LUCgIAEoXN0CDsiWw+7o/ElJbz5g5k4OrvT1FV0uL6abvN17JCufOIb5kSWz5dQva91+IABNLLN87C3OOLYZhbDTvMOrXZybRkSO0SX7wgB47mmmYxUkKFVL0JWkTHo4wA85e3830ExLoiFm5MvPKVe0Ns4Aq3j2yqRNKF9fH//72RGx8qkutWcfEhCEaT09MH78E16p/wipfFxc2KlmyhGGU1DA3px/OoUMMs1y+zEKqRPFTAGyq2Q7FIyPQ7t6FlKGg+/fTvKxuLqXhHx6Fq4/DUz5ZtSrQuzcXSd++Va893LmTct+M4ukJYWmJo0Mnoe3n0/Fjm1FwDbyDo+u+QRuvxCykMmWYAnvjBmshfHy4vlA8DUtoSaEkXdFXFGWtoijBiqKkmj6gKEpfRVFuJm6XFEVxTm0/SSEkJARhBqYpfXcuXOCt//TpFKhDh7J1Gn0dTfzUoQq8gl9j7YVH2bvmNAiPjMXFUlXwS58fGR9XzcpHjWJ8+osvaPaWlt2BpiYN2GbMAOLiMOyXA5gzaAZOlK+PHjdPQC8+lVl4xYrJB4H9+9891aaKNYx0tVLk7L/jp584s08quBcuZP6N+/pCGBvjbIcB6NB2MkZ0+RFxGhr4c99crN85DSVfh9O9FACCg4Hx42kot2ABXT8lRY6MzPTXA2jzgecfAWgihKgOYCaAlTlwXZKCQHAwwg2MYRYZkdx3p0YNYNs2Vo/q6gLt29Mh81HWBbtFpZL4tHJJLDzphScvsn73oCIoIgr7bgTix79u4dMF51Br5gl4Bb/G7ScvcSwwmoZh16+z9eKgQfSxb9oUKF8eWzt/hRG/7P/g8SNMSuBguQZI0NBA3z2LaTvw22/00kmLjh3fDQD6ulpoH3oXRzyDUk9ZdXICvvySVtKqdZTMiP7Tp4CODtwatUePz6dgYI8ZeK5fDL8eWoDja77G5w8uqauh4+OZifP4MWshSpbM+HkkhQ5FZMDISVGUMgAOCiGqprNfcQCeQog0/GPVuLq6Cvf3syokBYvduzFkmwcCTEri6Lpv+DMTk+SpkG/fMj3yp58oHpMm0XI3C92mAp5HouWCc2hWwRLL+mXOpKvj4gt4GR2HWnbFcfVxGPzDOXAY6WrBxb446jiY4aDHE/g9j0RsnMCqAa5oUt5CfYDISFb+rl4NnDuHBEUDGq1bsQCsY8cU76f78ku44f8CDZ3Mse7LOskvJiYGuHSJdQR//ZVmuuU1m4ro2u9XzDu8EF1unYTW+93C/PzYqOXtW35fpUr6+fqBgUDp0rhhXR6/NeqHfxxqoeSrMIy6tB09b56AjpaG+niGhgzjfPtt+n0GJAUGRVGuCSFcs/r6nI7pDwZwJK0nFUUZpiiKu6Io7iH5aSAlyRhBQQgzMIW5ahG3ePGUs0AdHWDCBApbx47AtGkM+ezfn2lnyNLFDfBN83I44vkUZ+4HA0hZ+fo+z9+8xcTdN+EREIFHoW9w5n4wKlsb43/tK+PgNw1xY+qn2DCoDr5u5gRjfW1ULFkMZS2NMHyTO64kXTg2MGC2z9mzGD1jJ/5q25+dtHr1Yi/dESOYFZT4np5HxiI2XqSepqmrCzRrxsKwu3c5SB46xN9THfUAUSvwHhzDArC7aguaYKnCQaoKYTs7nldFYCrZPgDXWQ4dAhQF3tXqYEiXKejUfwFulyyLKadX49zmMfji7mnoJMRR8EuVYmGVnx8wa5YU/I+MHBN9RVGagaI/Ma19hBArhRCuQghXCwuLtHaTFBR8fN6FdwBQlNK69S9dmumRJ08yDbJjR1bLHjmSKfEf2sgRjhaGmL7/NqJj0zZIE0Jgz7UAtFhwDruvB8DKWA/VbUxwbUpLrPjCFYMbOqCqjQm0NJN/xLU0NbBpcB3YmOpj8AZ3ePinLDp7ZlkaOzsMY7jj5El652/cyKyWSpXw90/L8Cj4FXS1NNC0QgZM0ExMmAE0bx49gl6+BI4cgTJxIrpGPMBVu2rwNbVS76+qEDYw4J2Tjg5//n7XsKdPudagqQm0b4/Tjq7oOGAhrthWxXfnN+L8zgkYcuso9F6+4N3HJ5/Qv9/fn+GcVFo9Soo+OSL6iqJUB7AaQEchRA7m3UnyFV9fhBmYqguzwsPTd3ps0YLx7TVrmPLYrh3F8vjxDIm/jpYGZnasCt+wSCw/l3qbQO+Q1+i7+grG7/KAfQkDHBjVEPYlDKCvo5mm701SzI10sWVIPRQ31Eb/tVdxNyiNwjMNDb6fTZuAp0/xetVajGswEGOj7FDT/w6WHPkNmps3ZT6Hvlgxxv7nzkWXVT8DAL7r/7PaP19FVBQXU1XhmIQEiveUKRwUrK1pJAdgeZ2uGNxtKuyfP8GOnT9i1L97YfTsCbOOunfnXcKFC6z4lYZoHzXZ/usrimIHYC+AL4QQD7J/SZKCQkzAE7zWNVCHd0JCMrbIp63NxdH795krHxRE7/uGDbn4m474f+Jkjg7OpbD0rHey2X50bDx+P/EAbRf+g1uBEfi5U1XsGdEAlUtl3jffykQPW4fUg762Jr5YcwU+IR8W7lsvBdo/L4O/S1bD2Fol0OPFPVT3vcNisJIlgb59eVfzflP1dLA20YeJvjauGVhh9dxNeBvxEti1K+0X6OkxJJNItKY2xn82DnObfYl2Xm7YtXUiKgd5U9hHjODdwc6dsoJW8o6MpGxuA3AZQAVFUQIURRmsKMoIRVFUwcapAEoAWKooyg1FUeTqbBEh/J0FQ2J4Jzw8c5kdOjos6nnwAFi2jDHkli2ZJXP27AdfOuWzStDR1MDjsEgIIXDpYSjaLvoHi055oW01K5wa3wT96tlDIzU//gxia2aAzUPqQgig7+orqdoiJCQIrDrvgy7LLuJtXAK2D6uPsT3qYU+Hofjm592cPffvT8Fv147iOmYMq1czGNZyKGGAYnra+PnQXbRZdQ1nqjVWdwr7AM+MzNCzz1zsrdoc4//ZhMV/zYaBkQHXDkJC+DtPzV9I8lGTrugLIXoLIayFENpCiNJCiDVCiOVCiOWJzw8RQhQXQtRI3LK8qiwpWIQlpp4n893JSiMPXV3OOr28WB378CEXOps1Y4ZLKuJoaayHcZ+WR0RULO4+fYU+q69ACIFNg+tgUa+asCyW+eyg1HCyNMKmwXXxJiYO/dZcwbOXapuGkFcxGLj+X8w6fBctKpbE4TGNUMchSRxcURgnX7aMdzN//03//xUruGBbsSKrax98+AZYV1sTFa2KYe1AVwgBfLnuXwxadxU+R84yYycVPKzKocOAhfAyt8Pyv+fgmxKRUA4cYK79vHl53jVMUniQwT1JmoRp0EHRPKnDZnZyuPX0WAzl7c00z/v3GWOuUYNxc1XsOpH+9e1hoKOJ19FxGN3cCUfHNkajcjmfAFC5lDE2DKqD0Fcx6Lf6CmLjE/Ai8i3aLjqPKz5hmNW5Kpb1qwVTA520D6Kry8XrXbu4wLp6NWf9U6cyDbNGDWbzeHmleYjmFUvi2NjGmOxsjKt3AtD6Sixml2uFV9t3sXgMAMzNsa9ma/ToMxdaisAek8doc2oHs6Xat6dlg0TyAaToS1InPh7hukYAkoR3gJwp3NHTo6fLo0fAunXM7+/fnyZf8+cDETyflqYGKlkVg7OtKca1qgA97dzzfqlpVxxrBtaGX3gkPANf4v6z1yhhqIsD3zRE37ppN0ZJFVNTunueOcNMmd9/Z078jz8C5cvTFyi1AcDPDzqDBmJY3yY4vX0COheLwqqKLdDsgQl2Po5CnKKBX6q2x5hW38DZQg/7Z3VHpWnfsW+ARJJBpOhLUic8HGGGpgCS+O4AOdunVVeXxU+3btHPp0IFFnbZ2rI/rL8/tDQ1oKuVMx/T9JpX1HMsgRVfuCAuIQGWxXSxb9QnKF+yWPZOWro0bZ4vXuSaxu+/M86eZADodGQDev+1lN/v3Al8/z0sb13DvKm9se/rT2CnROP7BgPhPG4XltXtht517LD5uzYoYSzj9ZLMI0VfkjpPniBM3wTa8bHJfXdyo0RfUYC2bZnZc+0a8+IXLgQcHTFq7XQ4+KbfQCSnuhE1rWAJV/vicDA3zPk7C1tbDgCXLnEAWLAA0NND730r0OnYZnb+evCAhVOmHHCrlzbFnn8WY+HFtRD6+ihTwgCzO1dlZa1EkgXkJ0eSOvfuIczQBMUjX6p9d/T12aw7N6lVi3bGPj7A6NFw9biAuXMG0Wd+1aqc9ZVPg0yFcrKKrS3tDy5fxsjZf2H0zJ0sALOzS76flxeU48fRqXVN3JnZFmcnNMub65MUWaToS1LHywvh+iYpQzt5JTh2dsBvv+GruX9jbc9xbIQ+bBgXR7/6Cvjvv7y5jjwgzKwknlmk0Vt3+XIuzg4ZkrcXJSmySNGXpI6PD8IMTFAiKocXcTNJlL4RjjXrxirfixdp6bx+Pe8I6tRh5e+bN+kep1ASFcWF7i5d6P8jkeQAUvQlqRMQgDADE3XmjoFBvoj+u1i9otCmYMMGdsNatIhiP2QI7QhGjqSvTSZN3rJ9Xbl6kh1s8DJyZO6eR/JRIUVfkjpPniDcIEl4x8goZzN3skPx4kz59PQE/vmHPWXXrmWTEycnetNkp8tUQWHpUnYoa9w4v69EUoSQoi9JlejnL/Ba10At+qGhBa+5hqLQz2fjRuDZM4ZCypYF5sxhJWuNGqxO9fPL7yvNPO7utHL46qu8W0eRfBRI0ZekSngcPxrJeuMWNNFPiokJc/6PH6ej5KJFLAKbOJENxhs1ol1CYenjsGwZQ2qqSlyJJIeQoi9JlfBEC4Ycr8bNC6ysGP5xc6PPz88/A2FhjI2XLEm/nLlz2SQlj9YAMsXz52xH2a8fBzOJJAeRRh2SlMTGIlSfhl25Vo2bV5Qty+rXyZOZAbR3L3DwIJuTTJoEODjQs+bzz2mWpqOT+wu075HifBs2MHNH1UFLIslB5ExfkpLAQIQbcIZZojDO9FNDUQBnZ/byvXYNCAhgDnzlyiz6atWKbQO7d+caQWho/lynEAzt1K/PNQmJJIeRoi9Jye3b70S/UIZ3MoKNDTB8OGf9YWHAvn3sh3vxoroxSqNGNIC7fz/vruv0aVoxyDRNSS4hRV+Skrt3EWrwnu+OpiZTJYsiBgZAhw7s8hUQAFy9ypDQ69c0gKtYkWZo330HnDuX6e5YmWLpUqBECaBbt9w7h+SjRoq+JCU+Pgg3MIFZUt8dS8uPo7eqhgZ9fmbMoNWDry+weDFj/3/8wa5flpbMqtm1K2Wz8uwQGMg7jsGDmXkkkeQCH8F/sSTTPH6McH2TohvayQx2dsDXXwPHjjHOv2sXF34PH6Yrprk5m5wvXcq7hOywahVTY4cPz5lrl0hSQYq+JCVBQQg1NEneMaswZu7kNMbGDLuoisHOn2dqqLc3BwZbW8DFhXcJHh6ZSweNjWV4qU0bNpORSHKJwin6+ZVZ8bEQHJw4008SuvhYZ/ppoaXFhd5ff+XC6507zP3X1QWmT2fmjYMDB4VDh4BXrz58vH372GdXLuBKcpnCJ/oHDgBlygC//cbZkSTnef6cZmsqh00dHSn6H0JRgEqVWP176RLFe80apoiuXs1wkJkZB4mffuI+7392ly1j5XDbtvnzHiQfDYVO9EfdfItrjs7MpKhVi4ZbkpwjLg7Rb+PwRtcA5m8Swztv38rwTmYoWRIYNIiz9/BwdgSbMAGIiaHof/IJM3Q6dODi8MGDTNUcPpxZUhJJLlLoRN/PpCS+7vUTYvf+zcyJxo3puRIcnN+XVjQIClLn6EfJ8E620dMDmjdnI/SrVxma3LMH6NsXuHsXGDOG1cDa2szakUhymUIn+tGx8Xj6MgYnneowjjppErB1K5tqL1sGxMfn9yUWbhJ99IEiXJiVn5iZsSnKsmWAlxfw6BFDQNu3y7spSZ5Q6ETfVF8bOpoKtl71AwwNOYO6eZOhnpEj6an+77/5fZmFlySinyx7R4p+7lCmDGf4Xbrk95VIPhIKnehX9L6FQW578Y9XKB6HJDbJrlgROHmSM/6AAKBuXQ4Az5/n78UWRh4/TmLBkCS8I2ehEkmRoNCJ/ozetfFlpBc0E+KxbcJvwL17fEJRgN69+f3o0cCKFeyitGgRFyIlGcPbG2EGpgCAEm+SzPQtLPLpgiQSSU5S6EQftWqh5IXTaGkYjd0m5RBTsxbj+qrm2CYmwMKFwPXrQM2awNix7KL0118F0zu9oOHjgzADY2jHx6LY20j+rEQJLjRKJJJCT+ETfQDQ1ESf3s0Qpm+CYwPGsyimUiV6pauE3dkZOHGChTHa2oyZNmki4/3pERCAMAPTlL47EomkSFA4RV8INHIyh62ZPrbWasdc/eLFga5dgXbtmBUBMOTTrh0XepctY+inTh2my/n65u97KKg8e5bYEF0u4kokRZHCJ/q3bwP160Pj3l30rmMHN59wPCxfg40xFi6kH3rVqsDUqUBkYnhCSwsYMYKt8yZP5h1BhQrADz8AEREfPt/HRFxcYjWucdFpniKRSJJR+EQ/OJgGV66u6H7nDLQ0FGy76kdhHzOGDS969ABmzqSwr12r9j83NgZmzeI+3bsDv/zCxd4FC9ie7mPn6VNAiHfhnXfI8I5EUmQofKLfqBG9Sxo2hMXXw9A6wht73P0RHZtYlGVtDWzaRAdEGxvmQDs7syReFe+3s+M+7u58bvx4iv/SpSyV/1hJtAYO1zeW4R2JpIhS+ET/+HEKtZkZMGgQ+h5djxfRcTi883Ty/Ro1Ai5fZsl7XBzQqRPQsCFw4YJ6HxcX5vefOUM726+/Tnl3b+tkdgAAIABJREFU8DHh749oLR280TWQ4R2JpIhS+ETfwQH48kuK9dq1qB/mDccXT7D1yA22uEvqXqgozNq5fZte5Y8ecTDo0AHw9FTv17Qp7wyOHmUoY/BgZgNt2fJx2ToEBCBMPxXfHRnekUiKDOmKvqIoaxVFCVYUxTON5xVFUf5QFOWhoig3FUWplfOXmYRx4xiGOXMG2L8fSuvW6H3zONxLV8b9lVs4Yz9zJvlrtLSAoUO5kDtnDgW+enUatfn5qd4I0Lo1cOUKQ0EGBkC/ftxvzx52NCrqBAQg3MAYwHuFWXKmL5EUGTIy018PoM0Hnm8LoFziNgzAsuxf1gc4epSP1aszHLNjB7ruWQYdRWBr8z6MSzdvDpQrR8/9mzfVsXwDA2bseHszjr99O/cbOVKdwqkovBP47z9gxw6KfbduLPTatq1oh338/RFqyObnJaJkeEciKYqkK/pCiPMAwj+wS0cAGwVxA2CqKIp1Tl1gClasUH/t7w8oCsyszdHO2QZ7KzVB5PmLbEbx8CE9952dubjbrx+wYQPw5AkrTOfPZ8ejgQPpcujkBAwZwtcBbJDdowfDQBs3MmzUpw9j/suXA9HRufYW8w1/f4Trq2b6SURfhnckkiJDTsT0bQD4J/k+IPFnKVAUZZiiKO6KoriHhIRk7Wzx8ZyNJz8w+tS1x6voOBzUs6VwT5vGsI6xMWfzx49T4G1saMswdixw6xbvBry9mce/eTNF/Ysv1J4+mpr83tOTVg7m5sBXX3FtYd48evoXFXx9k5itJYq+kRHvkCQSSZEgJ0RfSeVnqZrcCCFWCiFchRCuFlk18GrVCqhdO8WPazuWQDlTHWxR5exPn86UTCcnZuw0bEhbhnnzKPwrVrCNXfHizOxJSOBzAwYwhl+5MtCzJwcGgDP/Tp0ANzd2Qqpale3x7Oy4gFzYm7jExQHBwQg1MEnuuyNDOxJJkSInRD8AgG2S70sDeJIDx02dsmWZpz97NoU4EQVAn91/wsP/BTwDE2epzs5cmJ07Fzh8mOEaS0vg2DG2sTtxgm3sTEwYwhkzBli3js2tHRyAv//m2sFnn7HrEcC7jObN+dp//wVatuTicJkywDffqC0gChvPngHx8YkWDBHqkVyKvkRSpMgJ0d8PoH9iFk89ABFCiKAcOG7aaGrSWdPDAyhV6t2Pu3iehl5sNLYO+lG9r5YWZ+QeHpy9DxxIP56QEAr27NnsT/riBRd9V62ih4+BgTr98/BhevSbmQEdO7Lp9dWrDAXt3s0OXr168e6hfHke/+jRwpXx488IXZiBSfKOWTKeL5EUKTKSsrkNwGUAFRRFCVAUZbCiKCMURRmRuMthAD4AHgJYBWBkrl3t+1StCjx+zBk8AJOYN2h/7x/sq9wEr3Xfi0NXqMBUzT/+oEFblSo0YVMJs6YmUK0aF3NXr2ZY58ULhnKmTmX2zps3wP793KduXa4X2NlxfcDUlBYP/fpxQGjbls1d/vijcMT9E6txwxJn+u+QM32JpGghhMiXzcXFReQomzYJoaUlrpWqIOwnHhSbndsIAQixb1/KfR89EqJlSz7fpIkQd+5k7BxxcULs2SNE7dp8rba2EBUqCFGlihD6+vxZ0k1Li4+amkI4Owvxww9C7N4txMWLQty9K0RwsBCxsTn5W8g6v/8uBCAaDl8txrQfr34PU6fm95VJJJIkAHAX2dBerfwedHKMfv2A+vVRs0EDVHrmg6012qCPx1EoHTvy+aQNVMqUYTbPunXM169enbP1qVOBYsXSPoemJit8u3RhHv+iRczdj41lSKdPH6B0aRZ8qbZbt1gR7OHBLTVMTBg6KlGCm5kZNxMTXo+xMbfUvjY0BPT1GcbKDv7+gKamDO9IJEWcoiP6AFC2LBQfH/QZ9CP+5/ApPKzLo0bQAz6nKMmFX1GAQYOAzz9nwdavv7LH7m+/MWvn/bTQ96lZE1i/novEy5ZxO3SI6aFDhzIElFQwnz4F/vyTawYhIRTs2rUZAtLUBMLCuLgcFsYU0vBw2j5n1AZCS4vir68P6Ompv1Z9r63NTUsr+ab62blziNY1QKSOvgzvSCRFGEXkUwtBV1dX4e7univHfhUZA5dpR1DX7xZW75kB3fgkVbRt23Jh9n3c3Fjhe/060KwZBbpKlYyfNDoa2LmTon7hAsW0UycOAC1aqDONYmOBI0d4l3HwIFMl69XjANSzJwcDFULwuC9fcnv1KvnXERFcZ4iOpjX0+1vSn8fFpbsFaujjk/6LMffIH+h18ziv4fx5+hVJJJICgaIo14QQrlk+QHZiQ9nZcjym/x5N5p0W9hMPiv69Z4koLZ2U8fbUiIsTYtkyIYoXZzx+/HghXr7M/Mlv3xbi22+FMDPjuRwchJg1S4gnT5Lv9/SpEL/+KkTlytxPX1+IL74Q4swZIeLjM3/ebOLhVFPYTzwojjnVVf+e7t3L8+uQSCRpg2zG9Aufy2YGKWmsBwdzQ5y3c8aAQQvwWkc/+Q6KwmYqSdHUZGXugwd08lywgFk/W7dmrql65cp8bWAgX2tvzwIuW1ugc2fO8GNjGToZP57VvleuAP370+ytWTMWlf3vf+rK4DwgTNEBABnekUiKMEVW9AHAspguFvaqAfcSDug3ehUidA2T71CxIhdO38fcnFbMbm6s3u3bF2jcmP78mUFPD+jdm66f9+/TIfTiRa4jWFvTzuHCBQ4oderQ0ycoiA1eHB2ZAlqpEn3/f/uNg0hu8fYtwhJ/P+9EX0eHi8kSiaTokJ3bhOxsuR3eScoxzyBRbvJh0Wb6PhFiUSplqAcQIiEh9RfHxQmxcqUQVlbcr0uX7IU8YmKEOHBAiF691Gme9vZM57x1K/m+T54wlVKVIqooQjRrJsSqVUKEh2f9GlIjMFCsqNNZ2E88KCJ0DHi+0qVz9hwSiSTbQIZ30qdVFSusGeiKR3Ha6PntOjx1aZByJw0NzurfR1NT7cU/YwZTPatU4Sz96dPMX4yODj1/tm2j9cGmTZzNz5/P4jBnZ/bu9fPj3cDYsSz2un+fJnIBAbweKyuGinbs4KJudgkJQZi+CXTipO+ORFKUKXyi/+oVvW6SdsjKAI3KWWDjoLp4Fi3QvdsM+H/5FZ/Q1FTv9OQJY/2p5dMbGjLG7u1NwVfZMU+blnXRLVaM9QVHjvDcf/7J8/zwA9cB6tZlSuj9+7R3mDaNX//7LzONrlyh/YOFBe0hNm4Enj/P2rUEBSHM0ARmUdJ3RyIpyhQ+0d+7F5g8GWjThjntmaCOgxm2Dq2LV2/j0d2pK2YP/AlvNTRTxvVr1KD3zqNHKQ9iaUlxvnuXRmwzZtAEbsmSTA9EKY47ahTN5Ly96QkkBD2GKlbk3cDkyXQOdXHhQrG/P3DuHDB8OFNNBwzgcdq0UdcDZBRvb4Trm8iOWRJJEafwif6AASyKuniRxU03b2bq5dVLm2L7sHqISxBYa+2KET9s5Kwa4CxbRVQUF1NHj07dNtnJiaGVq1cZ7hk1ilk7Gzdmv7uWoyPF/upVCvuffzL0NG8eF3xtbXm+s2eB+vVZGezry4Xnb7+l0+ewYQwBNW3K16c2gCXl8WNW48reuBJJkabwiT5A4T93jsVH9evT/z4TVLQyxs7h9aChKDgXVwx+h05RZKOiaKmclD//5Ix3ypTUZ861a9Ol89AhDhoDBnBmvm5d9mb+KkqXpsCfPMnBZ8MGnnPtWrqEmpvTFmLlSor0vHlcf/jvP94ZBAdz4HJ05HWNG0db6JiY5OcJCJBmaxLJR0DhFH2A8W53dy5+dutG35xMWBk7WhihknUxCCEw7cgDiFmzOJAIwUXdpJWxANMn7e3pue/vn/w5RaH3zn//0YPfxIQVthUqMPb/9m0OvGHQj6d/f3bwCg3lY8+ewLVrXGdwdGTsf/RoLvhOnEjb53v3gN9/5/UvXcpGNGZm7AW8bBmdSoOC3nnpv0OKvkRS5Ci8og/QS//sWRZSzZzJbJZM2BjraWuidHF9nLkfguN3nrG7locHhfXlS8bOkxIVRatkR0eK+vvFXYrCBVV3d+DAAc7Chw6lEK9YkXJ2nR0MDGjzsHIlRfvePYZ5ypXjXcDnn1PYmzXj4FCvHovCwsL4OHAgi8JGjgQcHPD60hVE6uhLszWJpIhTuEUfYAHUmjUU40OHGO5RNTfPACWN9VDRqhhmHLiDyLdxnOGvW8fmKI8e0bCsdOnkL4qL4z6VKvEu49q15M8rCtMyr1xhZo61NSt9y5XjTDunm6orCu8qRo/m7yA8nH0Avv2W2TyTJvH3YmrK2b27O6/71q13dwHepcoCkNW4EkmRJztJ/tnZcqU469Qp+t2Ymgpx7FiGX3b1UZiwn3hQzD1yN/kTgYFCtG7NQqVWrVIv6jIx4eOnnwpx+nTqRV4JCUIcPy7EJ59w35IlhZg5U4iQkGy+4Qzy7Bl9/MeMEaJmTRZ5qfoB1K8v/t/emcdXUWR7/NtJIGyySYBENhcUUBhA1EFcEFBEFFBZZNT3nPGNy5MRHUXxjfoEEZWZp6I4IzouuLGpCIojKoiCoMIoi4DsIbImLCFkIzf39vvjnKIrzQ1LQIWkfp9Pf+69fU//urqq+vSpU6dP+fff77/f4Qq/6f0f+p+cdm5wbdu2/TLlc3BwOGTgXs6y0KWLxLA3bizZNJ944pBSE5/TrC59z27EP+esY02mFXOfliaW+vPPSyhltWoSomljt1rGixfL+du2FZdLXl4g43lw6aWyYtfMmdC+vcT8N24sI4Ay5tcZMHY+A8YeQmqI+vVlCchnnpHQzp07ZURwzz1StqeeouZOiVCqm58TlDleigoHB4fjGuVL6YP42+fNEyX3wAOiiDdsOOhhQ3u0oGqlRB56fxm+nVzN88TvvXy5pEh++GHx9TdvXpIgM1NSEHuexM2fdFIQPmlzdekiqZ2XLYMbb5Tw05YtJeZ/5szDS+xWVtSuLRPPjz8uoa/Z2bx89R0A1MvXOP2UlJIvrjk4OJQLlD+lD1CjhsTQv/qqRNS0aSPpDg6gUOvVSOa+y1swf90Opi3evL9A48aSAXPSJImMWbdOLGUbc+aIxf/mmzLSGDNGJnG7d5e1de1RR6tWMiLIyIBhw8TP3q2bjBTGjTuqk74HHRFUq8aGlCYAwUSu8+c7OJRLlE+lD2JV33STKOHWrSUiZ8CAA77FO/DcJrRpVIsR01eQUxgnxt7zoF8/eRv397+XzJcnnwwTJ1IiWPSGG2DCBAntHD5crPreveXN3SefLBnvX7++jB42bJAJ6WhUyt2smbiAMjKOTn0cBMUxn0rFEWoUFQTlcnBwKHcov0rf4OSTJf7+8cclhr51a0maFgeJCR4j+pzF9ty9PP3pqtI569SRNAezZ8sKWQMGMPe8y5l7auuScqmpotDXr4fJk6UsQ4fKXEGvXhIhZCJ5qlSRMNClS2HGDPH7P/aYHNOrl7iEDnXpxDIgEo1RLz/b5d1xcCjnKP9KH8Q3PXSohFDWri3uljvvlLj7ENo0qs315zVh3Lx0lm3eHYfMwsUXy0jiwQc5f8GntNmaITnxq4YWbKlcWeL1P/9cYuPvvlvcOf36yYPh1luDvPqeJy9PTZ8uLqShQyUdQ8+eMlIYObJs2T1LQXZ+EW9/k0FOQaRkjL5T+g4O5RIVQ+kbtGsnMfWDB0t6hbPPlmiWEIZc1oI61Srz0Ps/EIsdZGK1ShV49FGG/uU1NqU2k2icFi1EwdvIzgbPY8r1gxlw6tXi+vnkE4nnf/NNmQQ2WTvNewbNmom1n5EhcwmnnhqswNW/P41/XER2fhHbcw/P/59fVMy0xZv5r3ELOOexz/ifKUtJwufWb98LhJx7x8GhXKJcLox+SPj0U/GdZ2bCvfeK/7xatX1/v/Pvjdw7eTFPXtuaAec0OShd/xfmUVQcY8qJG/HuGyKTvQMHSs4ea4F1X7eEqVPFbQOQmyvZQ994I4jg6dhR0iZfe23JPP8rV8KLL/Lx50sZfMlt7K2UDED9akmc2bgOZ6bV4sy0mrRKq0mTutXwPI8BY+cT831u73wqUxdt5tPl28gvitKwZhV6tU2j12/SeHnsdJ5+9PrgPK+8IvMWDg4OxxSOdGH0iqv0QeLV77lHwiabNZN4/CuuAOSltQFjv2Z15h5OrledSokJTLy1Y4nDCyNR5q/dwSfLt/HOv38iEvVpWLMKXU6rQ9clszn/mUeo6kcl8Vn79vu4S8BW/iAPi7ffFut/6VJx93TqBP37w7XX4qem8vLc9Tw2fQUNi/Zw+4L3KCosYnn9U1h+SmtWV08hqp75E5KTaJlWk7WZuewuiFAc86ldrRI9zkqld9s0zm1Wl4QEkR0xeDQPPntXUI7p0+OX18HB4VeFU/pHA198IQnLVqwIXmJq1Igft+bQ89m51K1emVPqVWfirR3Jzi9i1o+ZfLp8G1+syiK/KEr1yokkJyVQo0oSZ6bV4stVWeQVRamS5NEpewNd5k6ja2QrDUc8zAvPjOfWLycHE6YGY8bIwiiwL7xyYue6MgE8aRIsXUpxQiLDbvhf3khtT49Ta7O92CMhwWPi5Wnw1lswbhyFa9ez+qTmLLvsGpa16ciyxFos3bSbmlUqMapvGy5snkLlpP29ei/+bgi3jP9bsMPk7XdwcDim4JT+0UJRkYRgDh8OSUkwYgTccQePzVjFS3PWk1arCk1OrMaC9F1EYz71T0imW6sGXNqqAeefeiL/8fK3AEy8tSN7i6N8u34nM1dk8tmKbWzcJRPGZ25dQ+vsjew4qx0vFX4n5wtjzBgGJLXfx2WQt3Q5gyYs4vNoLW795l3u/2IcK09tzTftOnPTk4Mlysf3ZdL39ddlOcZduyA1lWlnXcK8Dt144rGbZOQQB1Muv5GrZ7wZ7MjIkLkDBweHYwpO6R9trFsnFvfHH0O7duQ+93faf7SLomiM0xvU4NJWDbi0VUPanFRrn2sELOs85ALyfZ81mbl8tmwLs2Yv4d+FycQSErhg71Zu7nMuF38wjoSQ8o8Bj/S+k+HvjwZg6+5C/vDaAlZu28Pw3mdyfa0CmDyZjBdeo8nmdXLQWWdJZs1evWShlUhEXDSvv07xBx+SFIvK28r9+snWvn2JB8Cc87pz4bdWKGthISQnH8WKdXBwOBpwSv/ngO/LwiyDB8OWLUzoej3vX3YDE4Z0L/WQ0pR+GDf+9V/UWvIdC2o1ZluNEzkllsvvu7Xi2hmvU+2pwL3iAx6w/MbbuLllX3IKIoy5vj2XnBFE1QwYO58GWRt5tsYmeeP3yy8llr9+fYkK6tULunXjD2Pncs6iL7k96ztZjKW4eL8HwJKW59JmpbZHzZpBTiEHB4djCkeq9CtWyOahwvMk9fCKFTBoEP1nvs2LwwaIr/8I0yMU1axN1gVdmDvkEkbnLqRG5hYempVBx0oX8OT7i9gy8m/7FP7sk9vT78Qu+Bs3MvmNe0sofINtKY3grrtk9a6sLJkE7tJFHlp9+sCJJzLo1WFUKSqQfPtbt8qbv6efLu6lDh3gtNM4+SfrZTQXo+/gUG7hLP1DwNC/vMb17z1P6x8XSpTPo4/C734nK2wdJsIjAn/lShY+NoaX95zAJ807kpDgUadqEqmb01lWK40zstJ55Z3hNMy10kccSptFIpILaNo0tr79Dg2zNsn+Zs3k5bTu3eW9hVmzYNIkYjM+wcOXCeYLL5RRg4ODwzEH5975JfHpp7IE4fffw29+I6mbu3cvdXI0Hkp1A333HT8Ne4LXcmsx4Tfdyatclc7N6zGmU11qtDw9Ptkhtt0+N1DdLEnxMGuWvBuQmCjvA3TvLtE/JsVz374SNeTg4HDMwbl3fklceqmEMo4fD3v2SCbNrl0lh/+Ron17Gk+dxEMP38Ab05/g1cn/yz9HDKTG7JmSLiLeouyeJ9vzzx+UfltKI0kRPXWqJJ2bPRvuu0+4H3qoZE5/9zaug0O5hbP0y4qiIkmNPHy4KOR+/SRlQjjPfhkw4IV5tP1hPg989y7Mny8J2oYMkfw9CQmS8C3e3EJCQtmSsmVllVT0w4ZJojgHB4djDr+Ipe953uWe5630PG+N53lD4/zfxPO8zz3P+97zvCWe55X/VzkrV4ZBg2DtWsmX89FHshjKTTftv2D64cLzWNT6fFngZOZMWf/27rslFn/0aEkdUVwsidlsxGKB9X84ZUhJKfnbTeQ6OJRbHFTpe56XCDwP9ABaAQM9z2sVEnsQmOT7fjvgOuDvR7ugxyxOOAEeeUSU/5/+JG/PtmwpufuXLDkybrPS1qxZkoWzQwdZDaxpU5lMHj9e/PrLlu1/bIsWcnxYoR8KnHvHwaHc4lAs/XOBNb7vr/N9vwiYAPQOyfhATf1eC4iz9FQ5R4MG8PTTkJ4u6ZD/9S+Z7O3dW96SPVJ06iSjiYULoXNnccE0bSpunxo1RPnv2QP16pU8bvv2w7f+naXv4FBucShK/yTgJ+v3Rt1n4xHgBs/zNgIfAX86KqU7HlG/vuS837BBFPPcuXDeeeKKORphkGefDVOmyCjiyivlQXPKKTKy+OEH8c/HYhKvH4ax/g8WbeSUvoNDucWhKP14GiI8+zsQeM33/UbAFcAbnuftx+153i2e5y30PG9hVrxolPKEOnVkMjQ9HUaNksVWLr5YYuCnTxfFXAom3trxoG/20rq1uHfWrRN//4wZEn7529+Ki6lfP7H+MzP3X9QFAuU/cOD+/zml7+BQbnEoSn8jYGfeasT+7pubgUkAvu/PB6oAIT8D+L7/ou/7HXzf75BSFl/z8YgTThAXTHo6PPusfF55pUzOjh595OkOmjSBv/5VUjI/95yEY153nVj/o0ZJ8rj8fHnIfP31/sdPmFDS8k9OhurVj6xMDg4OxywORekvAJp7nney53mVkYnaaSGZDKArgOd5LRGlX85N+cNE1aoy0bt2rVjoKSmSPqFRI4kCsuPky4IaNYRn5UrJw3PaafIiWaNGkkBu+XJxM/k+5OWVzOFvo169w3rZzMHB4fjCQZW+7/vFwCBgBrACidJZ5nnecM/zjOa4B/ij53mLgfHATf6v9QLAsY7KlcUSnzdPXuq65hpZZL1lS3kz9sMPD+j6OSgSEiTb5qxZsGiRLL7yz39KFs6LLhJff2KivKTl+7Jou51OomHDI79GBweHYxbu5axjAZmZ8qLXP/4BmzfLWrh33AE33rh/NE5ZkJUlq4ONHSsjjXr1ZCnEW26REQFIxNGoUbIA+wcfHPk5HRwcfha4NAzlAfXry1q66eniY2/QAP78Z3kTt29fsf6Li8vOn5Ii8wqrVsli7BddBE89JW8PX3aZrM+7ZYu4dZyl7+BQruGU/rGESpUk9PKrryTaZ9AgCfO86irxzd97b/wXsQ4VCQmSP+jdd2VlrOHDgyUiX39d3EoucsfBoVzDKf1jFW3aiDW+aRO8/76EY44eLb75c86RJGs7d5adPy1NEq2tXy8Tv5dfLvvd27gODuUaTukf66hUSd7qnTJF/P3PPCO58gcNgtRUmQh+662yh34mJclI4m+6apez9B0cyjWc0j+ekJIiSzguWiQ5/W+/Hb75Bm64QSz0nj3hlVckVv9wkZkpn07pOziUazilf7yibVux+n/6SeYABg0Sf//NN4vi7tYNXnhBlkc8FGzbJp/OvePgUK7hlP7xjoQEOP98We92/XpJyHbffTJRe/vt4ru/6CL5f9Wq0nmM0neWvoNDuYZT+uUJnicJ2UaOlDdzlyyR/D/Z2RL5c8YZknRtyBBZP9cOA83MFP9+nTq/XvkdHBx+driXsyoK0tMl3n/aNFkqMRKBunVlHuCqqyRW/4svZLLYwcHhmIVbGN3h8JGTI1k5P/hAMn6a0M+2bWWC2MHB4ZjFkSr9pKNZGIfjBDVrSurlfv3ExTN/vij/du1+7ZI5ODj8zHBKv6IjKUly/F944a9dEgcHh18AbiLXwcHBoQLBKX0HBweHCgSn9B0cHBwqEJzSd3BwcKhAcErfwcHBoQLBKX0HBweHCgSn9B0cHBwqEJzSd3BwcKhA+NXSMHielwVsKOPh9YDtR0muInD9Gud0XI7LcR25XDw09X0/pYzHgu/7x90GLDxachWB63gvv+NyXBWR6+fanHvHwcHBoQLBKX0HBweHCoTjVem/eBTlKgLXr3FOx+W4HNeRyx11/GoTuQ4ODg4OvzyOV0vfwcHBwaEMcErfwcHBoSKhLCE/wF+AZcASYBFwHjAbKwwJ6ACsU7llwB4gH9gUklsIxIAosFR/R1V+uZ7jW2A3EFFZX2WKdNukclGLK13ljGwhsFV/x0JcZitWrtWhY3eH5Ixs3gG4zHkKQlxZpXDl6/fSuPKOIlc0xFVoHWNve6zr2Kt84evLQN63sLmK4nDNtb6XxvV/wDPW70K9njDX9oNw+UAuJftAaeXafRCuGPABsNGqu2XAT3G4dh2Ay/TLWGjfYkr2t3hbIbA21GbhckaVa7F+D1+DD2Qr1+JSym9kY8i9+GdgTYjDvj7DNc/6P3x9Ma2vP1vtEZYtUr7FwPdxym36tQ/8oFxbQzJ7Q2UrjcvuTxtUbjclyxsu/2I956Y45d+L3CcFqs/utf6PIv1mMaLLCoBzVW6otmGRHr8U0aULLZkhwEo9bq/y1o2jj18D1uvxi4C2B9TfZVD4HYH5QLL+rgekIUo/A+ih+3+vlZkMVAd6AvdrxWUAPZRrjfLlAk2A3sAqRAn1Ua5x2lhNEeWXbnFtA/ooVxTYocduBrZYXFcDswhu3rYESmWMykaBG7TijdLLB3bqeVcR3MzVtFw2l7nxNhM8xMLlytJGXKWfUZUJcxVrQ27W75v1OqOIcglzmXI2BWYehGsQ0Jngwdme4MYz11gAfENwozYFhiEdb6vWl1GIiwluqJ5aP9nKlav1uJtAaSUDz+nxq7UNi4HfETxofgS66/6dVrnwjvVuAAAL40lEQVTygC8RJWi4niK4OU25PrHK9Z1y2eXK0+v/uhSu1QR9Il3lTPsarl3AbcBLWsa5emxU++1zyrsOuUd2AJ8RKP89wP8oVwHSDz/Wa65G8CBehyiTqP53AUG751nX2RPI1Hp6G7kHC7V+jQJaB9QHcnTfw5ZcLnCGnmMyomCLtP4663k2qdwu3eorb8Tiyre4IsB0i2uhci1E+lcu8hDPsbiKlecB/S8PmKPXOFC5IsAXyvUNouCLdX96KVz5Wo95WkedgDu0HjZomz2hv03592q9rtY2+xa4TjkSgapIn9qpx/+kx6cBb+k19wSaablmAzcq/1LA0/bZqsdfAcy29O23iH64Cuk/PUpR+n1/zjj9VGC77/t7AXzf3+77/mb976/Ag/q9HhDxfX+v7/t5vu9P14uz5VIRRWS4Mnzfn6qVuwB50oHcUL7v++YN3lyLa4TKpSINnAxUAe7USjJc/we8iygfgBOAysgDyIw2tgKPI508BlRSrreRmyKG3FQx3/eNkrW5dhBY3+lIQ4fLNYbAAt8JFPu+Pz8OVyaBBbJNy5KgXJXicK0BkrSO0g7AtRxRVHVVxgfesa7ZlKsy8CjSKUGU/h/1vyICpdUV+FBlYoii+1HLGUPatwpiwRucC1yL3CwRld2LRDR8pDJVgH9oPZq6z9VrHxXiuk7LYqyhbEQxmnJVVy67XHuRm3NEHC5TrhiiKOpaclUsLjOCPE3r+Amtr4jneRfqNe4GGiMPkF2IIilSuU3AYKT+k5U/FcjT/gXSPtWRPp6n+75CFI6xKBO1LN2QPuIh7dQVUUJmWdR1+r0zQbvmW3JovewFzkL6gJFBZaro9zX6f2fl9CyuTRZXMdDI4irUzx26H+Se8UJceyyuZGCkXu9FccqVDbRA6jViXWOYK4cABUBf4L8R5ZvqeZ6HeC1iVvkLkHqtQzAq2Yq0xbl6PclAoed5/YGTgO9VJ45RHrP4dKKWaah+f8IXrb0eqO15XipQC7ln0N81VT8MRPRXH44UZbD0ayBDiFXA34GLdf9sxKUzC7gEuBC5SffJATchjWjkeiDDl3yksQzXj8BEi+tdAis0ot8zkAZYCfxbucywthgYgDRYIdIpioG7CCxXYzUZLnuo1ZOSw8q/s//wbxPBENNwGQVmZIsJrEjDNSzEFUWsiJ0hriJKDtN3IS6s0riiei33EAx143HFtH77WFwzrPPHLNnd1vcoYnXEQvsiyA1suN7TOoyFtusJrPgoYq2bOjLXU6TtaHPlh+rd1/badQAuI7/uAOUyciuRfloal4/0oaWWTCekT8Wsus1H+ry5lp2IFWm7H8yIL8/iugR5WPihLey+2YEYH751TjNiMA+o7YhxY8pvHkp79PsavdZMZERhuyoiWq42+jsHUWDmXEaB56pcDnJ/ZSL3lSmXKU8R8qDztVzxuAoI3Hh5Ia5Cq1wRq1yGq5jg/jXuyjw9R2lcUcSlEwtxGdejGXmbbYF+bgdetep0r8XVh2Bkk6XfcxAduYzADbvZaqOztIwXqL67hMDo3ISkWQDRk58ho76dyCjgw1Is/ZWIfnga9cIcNUvf9/1c4GzgFr2YiZ7n3WSJjECs+ALEgt4nh1hftty9yFBnFVLp8bjGAKcjw+tbVC4BGb5nIcpkCPA3Av97gh63A1HyexCL6o/IkzSC3NwJus3R8y1CnsKjkEbI0//7I422Vs+xhMA/6ltcHnJjRhC3xzIC90pMZe5QLnPNRnGaB8gCq1zj9b/lSMMnEdzgYS6jUO4isNa+jMP1KvLgNu6V3YgSq6xcqwmstoXACq2beUBrPeYdvaZdep5FBDdSV8SKWYooXQ9RNn9QHh95eLTXY6YRPJzSkT4B0n+66PcliKKP6TWtJrDabK5JVrnmIBaUUWxdQuUC6RebVA7EeDBckwmUyl6kr8X0esbr/x8D9+m+xcD5BEp/PvCQtk0OYiVvR/pfVT3fRkRpGJ+4aaOYVe8gbVGMtFEEaacFVn0mEowuTb/8CFFEyXqNCXpMptZlPqIgr9JjEpA+9p2eYxXS74sR1+Nteh7jespElHoGYvEW6bWkI+1ZCRkhmbkUm2uonq9Q626hythcdyJtm6jbAq2XjcpVpFyPKFeulqtQ2yrMZdxzyXpOM+Loj7SlccPYaKuf6UBzPf844HLdXxcxSBdqfWxC7tGvfd9vq3IRRCmfr+WqzP7W+o0EbqS7gZd1vxmNXYWM7ozRFMYDyEjnHC3T/XFkApRlIjf0lOmLKLjZQAfd9xXiN54dkvteK2Y/OaTRDJex9LshnXOo4UIUcYbNZcntVM5iRDnkKFcuolTNsN4H/pPAcslHFOYmS8624gzXNkTxm2OGKlcL5dqjcruRDmfkdhIof8O1UT9zEQU3VP+/OQ7Xc5bcTgIfrc21VbcoMEHPVz8OVwdLzrjJ1lhcZrLKKDxjcc7X/4w1s4uSrhljkWYj1sk24k80+4hC3oEojUICN1xYbpfFtYtgNFBoletDPX6zVa7dlpzhyrK4Mq02trnmxOEqtOTsa3xJz20eNHbfilnlWq///yf7BwREEQttAkHf2IO4RSEYEaQjCtD87qS8EW3DRbr/auSBVogYHDmIqwg9V7a2dYrWQRFiGY5XHjPpWoxYn4brVeShtRDpa8bAiMf1EjKnZEYxxcpjuMYj/Wg+4pvPIzAebK4pwFiCyXdjFG0jcPmNV/kIwYjQjA7ODHE9quf8ELmPTLlWIr73QuRBl67lydXy+4iS3qG8YxDX0Q5EF2xH7klTFzFgtTX/GdPPZnq9e/ScRcD1KmdGT6mIos/R/amI/pqCzHcNBMYeRB93Js5o4Igsfc/zzvA8r7m1qy37Z8t8DJmgqhqS2xGSe03l4nHVRhr9Tt1srioW1xkq9xhQU/cVAqcQ+DK/Qir1FgI3zHSko1UjcAPdp3IP6XHr9L/mylUTsSqStEz3qNx25aqOdJJcRHknWeUy1oXNla3nH6xbAtIhwly+yo3U42qojM2Vi8yjFCCuNI9AEdlcPVSuUDkTkeGm4ToB6fR5iBIxlkUOMlKqq+fIR5RUHvLgRmXP0S1Ry/kDMtHW0+LK1GMaIMrY1P032kYgN8dpIa61iJvmGotrr5aroZYxH1EG7wF/UpkogSWUqNe4AelDNld2HK43LDlzjR2QydzGyKh3ITIaMZPItqutKbDC9/1xyFxWDsGDYTSiaHprPZ6uZSgigI/cC7ch90BM6yISkoshFqN5oLVArPZsS+YdpO1N/YIELkSQ/uIpR5bWnTEArkGUeYrWXQHwvMW1y+KahQRkxJAHUx5i3Ro//1VIP66H9A0fscK9EFdXRNmZSeLfErgYf7K4hunvZ5H22q2fmSGur5C5rh8I5moSEWXaX+UeIHBV2vX6H8jDIBFR3g21HjYiD8kxQD3f97OQUeipnuddifTxGIHuqo609+W6f7jOI8SAPb7vb0FGpKsB9He+7puq5ZhKCOr7R7n66DWWjjJY9mcjw1ETTvke0oCzUQte5VZoAxg580SPIk/rVspl/IDGp2ZCM43/bAWBry0cAmc+cxFlaVtRhaHfxlowvtp4YX6+VnKYK17YoO0nj/e/KVu4HPHOGw5FjccVtqhL4yogiJIpjStCMFQ0x4TLuW/S2jp/vPNtJnCjGK5waGQEcTtFrd/x6uxpSvr0C0qR22TtN5N3YZkdlPTpxytXlCBCx5QrHtcSAndcMWJZZ8fh+ix0zfFCNndadbpVZVaGeMx12fyZofYMz3OYfUsJon6iITnzuT7UZjZvoSU7ksAvHz5XocU1wbpW+342o6Ri5doe4jKfpi+uVi67Dsx9ttyq15GUfHgaJW6i0UrjilKyTxQj+iXP4rLDam29MxLp60bOuBFHIVFEEeBkxNVs38ubEWt9uf7ubkXpGLdmrsosRgyflZYefQTpa2uRh4vJovARkKbfZ2m7/wC8CdQ4kA53aRgcHBwcKhDcG7kODg4OFQhO6Ts4ODhUIDil7+Dg4FCB4JS+g4ODQwWCU/oODg4OFQhO6Ts4ODhUIDil7+Dg4FCB8P+8rO7k/wKTXwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#plotting the y_predict based on energy E_c.m. for example the column SN112EK42.5 represent energy 42.5 \n",
"Yy = np.array([Y[:,0],Y[:,1],Y[:,2],Y[:,3],Y[:,4],Y[:,5],Y[:,6],Y[:,7],Y[:,8],Y[:,9],Y[:,10],Y[:,11],Y[:,12],Y[:,13],Y[:,14],Y[:,15],Y[:,16],Y[:,17],Y[:,18],Y[:,19],Y[:,20],Y[:,21],Y[:,22],Y[:,23],Y[:,24],Y[:,25],Y[:,26],Y[:,27],Y[:,28],Y[:,29],Y[:,30],Y[:,31],Y[:,32]])\n",
"y_ex_std =df4.loc[:,error_features].values\n",
"y_ex_std = y_ex_std.reshape(33,1)\n",
"Y_ex = df3.loc[:, features].values\n",
"Y_ex = Y_ex.reshape(33,1)\n",
"plt.plot(features,Yy, 'r')\n",
"plt.errorbar(features,Y_ex, yerr = y_ex_std)\n",
"plt.fill(features, Yy - np.sqrt(np.diag(Ycov)), Yy + np.sqrt(np.diag(Ycov)), alpha=0.5, color='k')\n",
"plt.xlabel('$E_{c.m.}$ (MeV)')\n",
"plt.ylabel('$R_{n/p}$ & $DR_{n/p}$')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$Like = exp(-\\sum \\frac{(y(x) - y_{exp})^2}{2\\sigma^2})$ \n",
"\n",
"$like = e^{-(y(x)-y_{exp})(YcovM)(y(x)-y_{exp}).T}$\n",
"\n",
"$YcovM = Ycov - Y_{experror}I$"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(47,)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAQiElEQVR4nO3df4zkdX3H8efLOygIGg5ZyZWDHrX4gzR62JViaRsFtShGMLVVYulhac6mpcH6E01TtbYJpiqa2Nqcgl5bRClCMYCtF4QQo57uwYmHh6JI9eDKrSIV2kY9ePeP+R6ue7u3szszu/uB5yOZ7Pf7+X5n53WzN6985zvf73xTVUiS2vO4pQ4gSVoYC1ySGmWBS1KjLHBJapQFLkmNWrmYD3bEEUfU2rVrF/MhJal5W7du/X5VjU0fX9QCX7t2LRMTE4v5kJLUvCT/OdO4u1AkqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSoyxwSWqUBS5JjbLAJalRi3ompuZn7QXXLnWERXfXhacvdQSpGW6BS1Kj+i7wJCuS3JLkmm7+2CRbktyR5JNJDhxdTEnSdPPZAj8f2DFl/t3ARVV1HPBD4NxhBpMk7V9fBZ5kDXA68JFuPsApwBXdKpuAM0cRUJI0s363wN8PvBl4uJt/EnB/Ve3p5ncCR810xyQbkkwkmZicnBworCTpZ+Ys8CQvBXZX1dapwzOsWjPdv6o2VtV4VY2Pje3zfeSSpAXq5zDCk4GXJXkJcBDwRHpb5IclWdltha8B7hldTEnSdHNugVfVW6tqTVWtBV4FfK6qXg3cALyiW209cPXIUkqS9jHIceBvAV6f5Fv09olfPJxIkqR+zOtMzKq6Ebixm74TOHH4kSRJ/fBMTElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSo/q5qPFBSb6c5KtJbkvyzm78Y0m+k2Rbd1s3+riSpL36uSLPj4FTqurBJAcAn0/ymW7Zm6rqitHFkyTNZs4Cr6oCHuxmD+huNcpQkqS59bUPPMmKJNuA3cDmqtrSLfrbJLcmuSjJL8xy3w1JJpJMTE5ODim2JKmvAq+qh6pqHbAGODHJrwJvBZ4OPAc4nN5V6me678aqGq+q8bGxsSHFliTN6yiUqrqf3lXpT6uqXdXzY+CjeIV6SVpU/RyFMpbksG76YOAFwO1JVndjAc4Eto8yqCTp5/VzFMpqYFOSFfQK//KquibJ55KMAQG2AX8ywpySpGn6OQrlVuCEGcZPGUkiSVJfPBNTkhplgUtSoyxwSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktSofi6pdlCSLyf5apLbkryzGz82yZYkdyT5ZJIDRx9XkrRXP1vgPwZOqapnAeuA05KcBLwbuKiqjgN+CJw7upiSpOnmLPDuyvMPdrMHdLcCTgGu6MY30buwsSRpkfS1DzzJiiTbgN3AZuDbwP1VtadbZSdw1Cz33ZBkIsnE5OTkMDJLkuizwKvqoapaB6wBTgSeMdNqs9x3Y1WNV9X42NjYwpNKkn7OvI5Cqar7gRuBk4DDkuy9qv0a4J7hRpMk7U8/R6GMJTmsmz4YeAGwA7gBeEW32nrg6lGFlCTta+Xcq7Aa2JRkBb3Cv7yqrknydeATSf4GuAW4eIQ5JUnTzFngVXUrcMIM43fS2x8uSVoCnokpSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSoyxwSWpUP5dUOzrJDUl2JLktyfnd+DuS3J1kW3d7yejjSpL26ueSanuAN1TVzUmeAGxNsrlbdlFVvWd08SRJs+nnkmq7gF3d9ANJdgBHjTqYJGn/5rUPPMlaetfH3NINnZfk1iSXJFk15GySpP3ou8CTHAp8CnhdVf0I+BDwFGAdvS30985yvw1JJpJMTE5ODiGyJAn6LPAkB9Ar70ur6kqAqrq3qh6qqoeBDzPLFeqramNVjVfV+NjY2LByS9JjXj9HoQS4GNhRVe+bMr56ymovB7YPP54kaTb9HIVyMnA28LUk27qxtwFnJVkHFHAX8NqRJJQkzaifo1A+D2SGRdcNP44kqV+eiSlJjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmN6ueamEcnuSHJjiS3JTm/Gz88yeYkd3Q/V40+riRpr362wPcAb6iqZwAnAX+W5HjgAuD6qjoOuL6blyQtkjkLvKp2VdXN3fQDwA7gKOAMYFO32ibgzFGFlCTta177wJOsBU4AtgBHVtUu6JU88ORZ7rMhyUSSicnJycHSSpIe0XeBJzkU+BTwuqr6Ub/3q6qNVTVeVeNjY2MLyShJmkFfBZ7kAHrlfWlVXdkN35tkdbd8NbB7NBElSTPp5yiUABcDO6rqfVMWfRpY302vB64efjxJ0mxW9rHOycDZwNeSbOvG3gZcCFye5Fzgu8DvjSaiJGkmcxZ4VX0eyCyLTx1uHElSvzwTU5IaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUqH4uqXZJkt1Jtk8Ze0eSu5Ns624vGW1MSdJ0/WyBfww4bYbxi6pqXXe7brixJElzmbPAq+om4L5FyCJJmodB9oGfl+TWbhfLqtlWSrIhyUSSicnJyQEeTpI01UIL/EPAU4B1wC7gvbOtWFUbq2q8qsbHxsYW+HCSpOkWVOBVdW9VPVRVDwMfBk4cbixJ0lwWVOBJVk+ZfTmwfbZ1JUmjsXKuFZJcBjwPOCLJTuDtwPOSrAMKuAt47QgzSpJmMGeBV9VZMwxfPIIskqR58ExMSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1Kj5izw7qrzu5NsnzJ2eJLNSe7ofs56VXpJ0mj0swX+MeC0aWMXANdX1XHA9d28JGkRzVngVXUTcN+04TOATd30JuDMIeeSJM1hofvAj6yqXQDdzyfPtmKSDUkmkkxMTk4u8OEkSdON/EPMqtpYVeNVNT42Njbqh5Okx4yFFvi9SVYDdD93Dy+SJKkfCy3wTwPru+n1wNXDiSNJ6lc/hxFeBnwReFqSnUnOBS4EXpjkDuCF3bwkaRGtnGuFqjprlkWnDjmLJGkePBNTkhplgUtSoyxwSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktSoOS/osD9J7gIeAB4C9lTV+DBCSZLmNlCBd55fVd8fwu+RJM2Du1AkqVGDFngBn02yNcmGYQSSJPVn0F0oJ1fVPUmeDGxOcntV3TR1ha7YNwAcc8wxAz6cJGmvgbbAq+qe7udu4CrgxBnW2VhV41U1PjY2NsjDSZKmWHCBJzkkyRP2TgMvArYPK5gkaf8G2YVyJHBVkr2/5+NV9e9DSSVJmtOCC7yq7gSeNcQskqR58DBCSWqUBS5JjbLAJalRFrgkNcoCl6RGDePLrCRpXtZecO1SR1h0d114+tB/p1vgktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUqGbOxHwsnrn1WOTfWeqfW+CS1KiBCjzJaUm+keRbSS4YVihJ0twGuajxCuDvgRcDxwNnJTl+WMEkSfs3yBb4icC3qurOqvoJ8AngjOHEkiTNZZAPMY8Cvjdlfifw69NXSrIB2NDNPpjkGwM85kIcAXx/kR9zvsw4HGYcjhYyQhs5H8mYdw/0e35ppsFBCjwzjNU+A1UbgY0DPM5AkkxU1fhSPX4/zDgcZhyOFjJCGzlHnXGQXSg7gaOnzK8B7hksjiSpX4MU+FeA45Icm+RA4FXAp4cTS5I0lwXvQqmqPUnOA/4DWAFcUlW3DS3Z8CzZ7pt5MONwmHE4WsgIbeQcacZU7bPbWpLUAM/ElKRGWeCS1KhmCzzJJUl2J9k+y/JXJ7m1u30hybO68aOT3JBkR5Lbkpy/3DJOWb4iyS1JrlmOGZMcluSKJLd3z+dzl2HGv+j+ztuTXJbkoCXKeEaXb1uSiSS/OWXZ+iR3dLf1o8g3SMYk65J8sXseb03yylFlHCTnlOVPTHJ3kg8ux4xJjkny2e418/UkaxccpKqavAG/DTwb2D7L8t8AVnXTLwa2dNOrgWd3008Avgkcv5wyTln+euDjwDXL7Xns5jcBf9xNHwgctpwy0jvZ7DvAwd385cA5S5TxUH72mdMzgdu76cOBO7ufq7rpVcss41OB47rpXwR2jepvPUjOKcs/0L1uPrgcMwI3Ai+cst7jF5qj2S3wqroJuG8/y79QVT/sZr9E7zh1qmpXVd3cTT8A7KD3Ql82GQGSrAFOBz4yimyDZkzyRHr/iS/u1vtJVd2/nDJ2VgIHJ1kJPJ4RnavQR8YHq3vFAofws5PefgfYXFX3df+GzcBpyyljVX2zqu7opu8BdgNjo8g4SE6AJL8GHAl8dlT5BsnYfV/UyqraPGW9/11ojmYLfJ7OBT4zfbB763ICsGWR88xkesb3A28GHl6aODOamvGXgUngo91uno8kOWTpoj3ikYxVdTfwHuC79LYa/7uqRvrC3p8kL09yO3At8Efd8ExfSTGSDYp+zJJx6vIT6b3b+vZiZ5uWY5+cSR4HvBd401Jm22uW5/KpwP1JruxeN3+X3hcDLsijvsCTPJ/ei/ot08YPBT4FvK6qfrQU2aZk+bmMSV4K7K6qrUuZa6oZnseV9N5CfqiqTgD+B1jSrxSe4XlcRe8L1o6l99b/kCR/sFT5quqqqno6cCbwrm64r6+kWCyzZAQgyWrgn4HXVNWSbljMkvNPgeuq6nuz33PxzJJxJfBbwBuB59DbEDpnoY/xqC7wJM+ktwvijKr6wZTxA+iV96VVdeVS5euyzJTxZOBlSe6i9y2PpyT5lyWKOFvGncDOqtr77uUKeoW+JGbJ+ALgO1U1WVU/Ba6kt798SXVvv5+S5AiW6VdSTMu4d5fZtcBfVtWXljTcFNNyPhc4r3vdvAf4wyQXLmU+mPHvfUv1vsV1D/BvDPC6edQWeJJj6L1gz66qb04ZD739tjuq6n1Lla/LMmPGqnprVa2pqrX0vqLgc1W1JFuO+8n4X8D3kjytGzoV+PoSRJw1I71dJycleXz3dz+V3mceS5HxV7oMJHk2vd0QP6B3JvOLkqzq3jG8qBtbNhnT+6qMq4B/qqp/XYpsU82Ws6peXVXHdK+bN9LLuyTvCvfz9/4KsCrJ3s8QTmGA100z18ScLsllwPOAI5LsBN4OHABQVf8I/BXwJOAfuudxT/W+Fexk4Gzga0m2db/ubVV13TLKuGgGzPjnwKXdC/xO4DXLKWNVbUlyBXAzsAe4hRGd2txHxt+lt0X4U+D/gFd2H3Ldl+Rd9F7YAH9dVbN+OLYUGZP8Pr0PrJ+U5Jzu151TVdsYgQGey0UzQMaHkrwRuL4r+K3AhxecY5H/3ZKkIXnU7kKRpEc7C1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ16v8B4MaMLvLtaVQAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#creating the likelihood function\n",
"Y_exp = df3.loc[:, features].values\n",
"y_exp_std =df4.loc[:,error_features].values\n",
"covM = Ycov - y_exp_std*np.eye(33)\n",
"Z = Y-Y_exp\n",
"H = np.dot(Z,covM) #matmul?\n",
"J = np.dot(H,Z.T) #J = (Y-Y_exp)*covM*(Y-Y_exp).T, Should be inverse covM\n",
"like = np.exp(-sum(J))\n",
"#like = np.exp(-sum(((Y-Y_exp)**2)/(2*y_exp_std**2)))\n",
"plt.hist(like,5)\n",
"print(like.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$Posterior = \\frac{Like*S0*L*ms*mv}{\\sum(like*(S0+L+ms+mv)}\n",
"\\\\\n",
"\\\\\n",
"posterior \\propto Like*S0*L*ms*mv$"
]
},
{
"cell_type": "code",
"execution_count": 87,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.387185\n"
]
}
],
"source": [
"N = 10000 #number of iterations\n",
"\n",
"#size of the steps\n",
"met_step = 1\n",
"met_step2 = 0.05\n",
"met_step3 = 0.01\n",
"met_accept = 0\n",
"\n",
"#empty arrays to store the data from the mcmc\n",
"S0 = np.zeros(N)\n",
"L = np.zeros(N)\n",
"ms = np.zeros(N)\n",
"mv = np.zeros(N)\n",
"lik = np.zeros(N)\n",
"mu = np.zeros(N)\n",
"\n",
"#this is the likelihood function that I derived in the above line\n",
"data = like\n",
"\n",
"#starting points\n",
"S0[0] = 30\n",
"L[0] = 70\n",
"mu[0] = 0.8\n",
"lik[0] = 0.9\n",
"mv[0] = 1\n",
"\n",
"#Here is the random walk\n",
"for i in range(N-1):\n",
" if np.random.rand() > 0.5:\n",
" S0_c = S0[i] + uniform(0,met_step).rvs()\n",
" L_c = L[i] + uniform(0,met_step).rvs()\n",
" mv_c = mv[i] + uniform(0,met_step2).rvs()\n",
" else:\n",
" S0_c = S0[i] - uniform(0,met_step).rvs()\n",
" L_c = L[i] - uniform(0,met_step).rvs()\n",
" mv_c = mv[i] - uniform(0,met_step2).rvs()\n",
" \n",
" mu_c = norm(mu[i],met_step3).rvs() \n",
" ms_i = norm(0.7, 0.05).pdf(mu[i])\n",
" ms_c = norm(0.7, 0.05).pdf(mu_c)\n",
" lik_i = norm(mu[i],1.4).pdf(data).prod()\n",
" lik_c = norm(mu_c,1.4).pdf(data).prod()\n",
" \n",
" #Here we have set up the parameter range\n",
" constraints = [0.6<=ms_c<=1,\n",
" 0.6<mv_c<1.2,\n",
" 25.7<S0_c<36,\n",
" 32<=L_c<=120]\n",
" if all(constraints):\n",
" #the metroplois ratio\n",
" met_r = lik_c*S0_c*L_c*ms_c*mv_c/lik_i*S0[i]*L[i]*ms_i*mv[i]\n",
" if np.random.rand() < min(1,met_r):\n",
" S0[i+1] = S0_c\n",
" L[i+1] = L_c\n",
" mu[i+1] = mu_c\n",
" mv[i+1] = mv_c\n",
" ms[i] = ms_c\n",
" lik[i] = lik_c\n",
" met_accept = met_accept +1\n",
" else:\n",
" S0[i+1] = S0[i]\n",
" L[i+1] = L[i]\n",
" mu[i+1] = mu[i]\n",
" mv[i+1] = mv[i]\n",
" ms[i] = ms_i\n",
" lik[i] = lik_i\n",
" else:\n",
" S0[i+1] = S0[i]\n",
" L[i+1] = L[i]\n",
" mu[i+1] = mu[i]\n",
" mv[i+1] = mv[i]\n",
" ms[i] = ms_i\n",
" lik[i] = lik_i\n",
" \n",
" \n",
"M = met_accept/N\n",
"print(M)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"#Using the MCMC sampling figure 2 can be created using the plots below\n",
"#plt.hist(S0,30, range = (10,40))\n",
"#plt.hist(L,30, range=(30,120))\n",
"#plt.hist(ms,30,range =(0.6,1))\n",
"#plt.hist(mv,100, range = (0.6,1.2))\n",
"#plt.hist2d(ms,mv)\n",
"#plt.hist2d(ms,L)\n",
"#plt.hist2d(ms,S0)\n",
"#plt.hist2d(mv,ms)\n",
"#plt.hist2d(mv,L)\n",
"#plt.hist2d(mv,S0)\n",
"#plt.hist2d(L,ms)\n",
"#plt.hist2d(L,mv)\n",
"#plt.hist2d(L,S0)\n",
"#plt.hist2d(S0,mv)\n",
"#plt.hist2d(S0,L)\n",
"#plt.hist2d(S0,ms)"
]
},
{
"cell_type": "code",
"execution_count": 105,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
" C:\\Users\\danny\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:1: RuntimeWarning:divide by zero encountered in true_divide\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"inf\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD7CAYAAACG50QgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAVWElEQVR4nO3df4xl5X3f8ffHEHDb1GYxC8ULzoK8cUxUGdAIaC01Nlj8SuWlKiTrNvWaUm2dkshVW9VQRyLFdov7Rx1bSXBoIMZOwo+QWmxjErrmh6JKBrPENjYQvAt2zXY37DoLxC4KNeTbP+4z5LLcmbmze+fuzDzvlzS65zznOec+55nZ733u9zznbKoKSVIfXne4GyBJmh6DviR1xKAvSR0x6EtSRwz6ktQRg74kdWSsoJ/kmCR3JPnTJI8n+XtJjk2yLcmO9rqm1U2STyfZmeSRJGcOHWdzq78jyealOilJ0mjjjvQ/BfxRVf0E8A7gceAq4J6q2gDc09YBLgI2tJ8twPUASY4FrgHOBs4Crpn9oJAkTUcWujkryRuArwOn1lDlJE8A76qqPUlOBO6vqrcl+Y22fMtwvdmfqvqXrfxV9UY57rjjav369YdwepLUn4cffvh7VbV21LYjx9j/VGAf8FtJ3gE8DHwIOKGq9gC0wH98q78OeHpo/12tbK7yOa1fv57t27eP0URJ0qwk/3uubeOkd44EzgSur6ozgP/LX6dyRr7fiLKap/zVOydbkmxPsn3fvn1jNE+SNK5xgv4uYFdVPdjW72DwIfBMS+vQXvcO1T95aP+TgN3zlL9KVd1QVTNVNbN27chvJ5Kkg7Rg0K+qPwOeTvK2VnQe8BiwFZidgbMZuLMtbwXe32bxnAM839JAdwPnJ1nTLuCe38okSVMyTk4f4BeB30lyFPAUcDmDD4zbk1wBfBe4rNW9C7gY2Am80OpSVfuTfBR4qNW7tqr2T+QsJEljWXD2zuE0MzNTXsiVpMVJ8nBVzYza5h25ktQRg74kdcSgL0kdMehLUkfGnb0jaQrWX/XFV5a/c91PH8aWaLVypC9JHTHoS1JHTO9Iy5SpHi0FR/qS1BGDviR1xKAvSR0x6EtSRwz6ktQRg74kdcSgL0kdcZ6+dJgNz8cfp45z9nUoHOlLUkcM+pLUEYO+JHXEoC9JHTHoS1JHDPqS1BGDviR1xKAvSR0x6EtSR7wjV1rBvFNXi+VIX5I6MtZIP8l3gO8DLwMvVdVMkmOB24D1wHeAn6mqZ5ME+BRwMfAC8IGq+pN2nM3AL7XDfqyqbp7cqUh9GOdZPdJcFpPeeXdVfW9o/Srgnqq6LslVbf3DwEXAhvZzNnA9cHb7kLgGmAEKeDjJ1qp6dgLnIa0oBm4dLoeS3tkIzI7UbwYuGSr/XA08AByT5ETgAmBbVe1vgX4bcOEhvL8kaZHGDfoF/M8kDyfZ0spOqKo9AO31+Fa+Dnh6aN9drWyucknSlIyb3nlnVe1OcjywLcmfzlM3I8pqnvJX7zz4UNkC8Ja3vGXM5kmSxjHWSL+qdrfXvcAXgLOAZ1rahva6t1XfBZw8tPtJwO55yg98rxuqaqaqZtauXbu4s5EkzWvBkX6SvwW8rqq+35bPB64FtgKbgeva651tl63ALyS5lcGF3Oerak+Su4H/lGRNq3c+cPVEz0Zaxrx4q+VgnPTOCcAXBjMxORL43ar6oyQPAbcnuQL4LnBZq38Xg+maOxlM2bwcoKr2J/ko8FCrd21V7Z/YmUid80YtjWPBoF9VTwHvGFH+58B5I8oLuHKOY90E3LT4ZkqSJsE7ciWpIz57R1qFTPVoLo70JakjBn1J6ojpHWkJOU1Ty41BX1rlzO9rmOkdSeqIQV+SOmJ6R5ow8/hazhzpS1JHDPqS1BGDviR1xKAvSR0x6EtSRwz6ktQRp2xKE7BSpml6d64c6UtSRwz6ktQRg74kdcScvtQp8/t9cqQvSR0x6EtSRwz6ktQRg74kdcQLuZJew4u8q5dBXzpIK+UuXGmY6R1J6sjYQT/JEUm+muQP2vopSR5MsiPJbUmOauVHt/Wdbfv6oWNc3cqfSHLBpE9GkjS/xaR3PgQ8DryhrX8C+GRV3ZrkM8AVwPXt9dmqemuSTa3ezyY5DdgE/CTwZuBLSX68ql6e0LlIOkimqvox1kg/yUnATwO/2dYDnAvc0arcDFzSlje2ddr281r9jcCtVfViVX0b2AmcNYmTkCSNZ9yR/q8A/x742239TcBzVfVSW98FrGvL64CnAarqpSTPt/rrgAeGjjm8j7QiOCLWSrfgSD/JPwT2VtXDw8UjqtYC2+bbZ/j9tiTZnmT7vn37FmqeJGkRxknvvBN4b5LvALcySOv8CnBMktlvCicBu9vyLuBkgLb9jcD+4fIR+7yiqm6oqpmqmlm7du2iT0iSNLcFg35VXV1VJ1XVegYXYu+tqn8K3Adc2qptBu5sy1vbOm37vVVVrXxTm91zCrAB+MrEzkSStKBDuTnrw8CtST4GfBW4sZXfCHw+yU4GI/xNAFX1aJLbgceAl4ArnbkjSdOVwSB8eZqZmant27cf7mZIr+j9Qq6PZFgZkjxcVTOjtnlHriR1xKAvSR3xgWvSAnpP6Wh1MehLIxjotVqZ3pGkjhj0JakjBn1J6ohBX5I6YtCXpI4Y9CWpI07ZlDS24amsPpJhZXKkL0kdcaQv6ZD5DWDlcKQvSR0x6EtSR0zvSDooPp9oZXKkL0kdcaSvrjlanTwv6i5vjvQlqSOO9NUVR/bqnSN9SeqIQV+SOmLQl6SOGPQlqSMGfUnqiEFfkjpi0JekjiwY9JO8PslXknw9yaNJ/mMrPyXJg0l2JLktyVGt/Oi2vrNtXz90rKtb+RNJLliqk5IkjTbOSP9F4NyqegdwOnBhknOATwCfrKoNwLPAFa3+FcCzVfVW4JOtHklOAzYBPwlcCPx6kiMmeTKSpPkteEduVRXwg7b6I+2ngHOBf9LKbwZ+Gbge2NiWAe4AfjVJWvmtVfUi8O0kO4GzgC9P4kSkuXgXrvTXxnoMQxuRPwy8Ffg14Enguap6qVXZBaxry+uApwGq6qUkzwNvauUPDB12eB9Jq5APX1t+xrqQW1UvV9XpwEkMRudvH1WtvWaObXOVv0qSLUm2J9m+b9++cZonSRrTombvVNVzwP3AOcAxSWa/KZwE7G7Lu4CTAdr2NwL7h8tH7DP8HjdU1UxVzaxdu3YxzZMkLWCc2TtrkxzTlv8G8B7gceA+4NJWbTNwZ1ve2tZp2+9t1wW2Apva7J5TgA3AVyZ1IpKkhY2T0z8RuLnl9V8H3F5Vf5DkMeDWJB8Dvgrc2OrfCHy+Xajdz2DGDlX1aJLbgceAl4Arq+rlyZ6OJGk+48zeeQQ4Y0T5Uwzy+weW/yVw2RzH+jjw8cU3U5I0Cd6RK0kd8X/O0qrk3HxpNIO+Vg0DvbQw0zuS1BGDviR1xPSOpMPKRzVMl0Ff0lQY3JcH0zuS1BGDviR1xPSOVjSnaUqL40hfkjriSF/SsuHF3qVn0Jc0dablDh/TO5LUEYO+JHXE9I5WHFMD0sFzpC9JHTHoS1JHDPqS1BGDviR1xAu5WhG8eCtNhkFfy4p3ZEpLy6AvaVlyALA0zOlLUkcc6WvZMo8vTZ5BX9KyZ6pnckzvSFJHFgz6SU5Ocl+Sx5M8muRDrfzYJNuS7Giva1p5knw6yc4kjyQ5c+hYm1v9HUk2L91pSZJGGWek/xLwb6vq7cA5wJVJTgOuAu6pqg3APW0d4CJgQ/vZAlwPgw8J4BrgbOAs4JrZDwpJ0nQsGPSrak9V/Ulb/j7wOLAO2Ajc3KrdDFzSljcCn6uBB4BjkpwIXABsq6r9VfUssA24cKJnI0ma16Iu5CZZD5wBPAicUFV7YPDBkOT4Vm0d8PTQbrta2Vzl6pyzdLQYXtQ9NGNfyE3yo8DvA/+6qv5ivqojymqe8gPfZ0uS7Um279u3b9zmSZLGMFbQT/IjDAL+71TVf2/Fz7S0De11byvfBZw8tPtJwO55yl+lqm6oqpmqmlm7du1izkWStIBxZu8EuBF4vKr+69CmrcDsDJzNwJ1D5e9vs3jOAZ5vaaC7gfOTrGkXcM9vZZKkKRknp/9O4J8B30jytVb2H4DrgNuTXAF8F7isbbsLuBjYCbwAXA5QVfuTfBR4qNW7tqr2T+QsJEljWTDoV9X/YnQ+HuC8EfULuHKOY90E3LSYBmp18uKtdHh4R64kdcSgL0kd8YFrmhpTOtLh50hfkjriSF/SiuXduYvnSF+SOmLQl6SOGPQlqSPm9CWtCub3x+NIX5I6YtCXpI4Y9CWpI+b0taS8C1daXhzpS1JHHOlLWnWcyTM3g74mzpSOtHyZ3pGkjhj0Jakjpnc0EaZ0pJXBkb4kdcSgL0kdMb0jaVVz+uarOdKXpI4Y9CWpI6Z3dNCcsSOtPI70JakjBn1J6siCQT/JTUn2JvnmUNmxSbYl2dFe17TyJPl0kp1JHkly5tA+m1v9HUk2L83pSJLmM85I/7PAhQeUXQXcU1UbgHvaOsBFwIb2swW4HgYfEsA1wNnAWcA1sx8UkjQt66/64is/vVrwQm5V/XGS9QcUbwTe1ZZvBu4HPtzKP1dVBTyQ5JgkJ7a626pqP0CSbQw+SG455DPQVPX8j0VaDQ529s4JVbUHoKr2JDm+la8Dnh6qt6uVzVWuFcBAL60ek76QmxFlNU/5aw+QbEmyPcn2ffv2TbRxktS7gw36z7S0De11byvfBZw8VO8kYPc85a9RVTdU1UxVzaxdu/YgmydJGuVgg/5WYHYGzmbgzqHy97dZPOcAz7c00N3A+UnWtAu457cySToser2ou2BOP8ktDC7EHpdkF4NZONcBtye5AvgucFmrfhdwMbATeAG4HKCq9if5KPBQq3ft7EVdSdL0jDN7531zbDpvRN0CrpzjODcBNy2qdTpsehv9SL3wjlxJ6ogPXOuczxqX+uJIX5I6YtCXpI6Y3tErvHirXvWU5nSkL0kdMehLUkcM+pLUEYO+JHXEoC9JHXH2ToecpSP1y5G+JHXEoC9JHTG90wlTOpLAkb4kdcWR/irlyF7SKAZ9SRqy2p/DY3pHkjriSH+FW+2jEkmTZdBfRczjS1qI6R1J6ogjfUmaw4HfnldDCtWRviR1xKAvSR0xvbMCecFW0sFypC9JHTHoS1JHpp7eSXIh8CngCOA3q+q6abdhJTKlI2kSpjrST3IE8GvARcBpwPuSnDbNNkhSz6ad3jkL2FlVT1XV/wNuBTZOuQ2S1K1pp3fWAU8Pre8Czp5yG1YMUzrS8rIannU17aCfEWX1qgrJFmBLW/1BkicO4f2OA753CPsvFdu1OLZrcWzX4hxUu/KJJWjJqx1Kf/3YXBumHfR3AScPrZ8E7B6uUFU3ADdM4s2SbK+qmUkca5Js1+LYrsWxXYvTW7umndN/CNiQ5JQkRwGbgK1TboMkdWuqI/2qeinJLwB3M5iyeVNVPTrNNkhSz6Y+T7+q7gLumtLbTSRNtARs1+LYrsWxXYvTVbtSVQvXkiStCj6GQZI6sqKDfpLLkjya5K+SzHmVO8mFSZ5IsjPJVUPlpyR5MMmOJLe1i8uTaNexSba1425LsmZEnXcn+drQz18muaRt+2ySbw9tO31a7Wr1Xh56761D5Yezv05P8uX2+34kyc8ObZtof8319zK0/eh2/jtbf6wf2nZ1K38iyQWH0o6DaNe/SfJY6597kvzY0LaRv9MptesDSfYNvf+/GNq2uf3edyTZPOV2fXKoTd9K8tzQtqXsr5uS7E3yzTm2J8mnW7sfSXLm0LZD76+qWrE/wNuBtwH3AzNz1DkCeBI4FTgK+DpwWtt2O7CpLX8G+PkJteu/AFe15auATyxQ/1hgP/A32/pngUuXoL/GahfwgznKD1t/AT8ObGjLbwb2AMdMur/m+3sZqvOvgM+05U3AbW35tFb/aOCUdpwjptiudw/9Df38bLvm+51OqV0fAH51xL7HAk+11zVtec202nVA/V9kMLFkSfurHfsfAGcC35xj+8XAHzK4r+kc4MFJ9teKHulX1eNVtdDNWyMf/ZAkwLnAHa3ezcAlE2raxna8cY97KfCHVfXChN5/Lott1ysOd39V1beqakdb3g3sBdZO6P2HjfOokOH23gGc1/pnI3BrVb1YVd8GdrbjTaVdVXXf0N/QAwzug1lqh/JolQuAbVW1v6qeBbYBFx6mdr0PuGVC7z2vqvpjBoO8uWwEPlcDDwDHJDmRCfXXig76Yxr16Id1wJuA56rqpQPKJ+GEqtoD0F6PX6D+Jl77B/fx9tXuk0mOnnK7Xp9ke5IHZlNOLKP+SnIWg9Hbk0PFk+qvuf5eRtZp/fE8g/4ZZ9+lbNewKxiMFmeN+p1Os13/uP1+7kgye4PmsuivlgY7Bbh3qHip+mscc7V9Iv217P/nrCRfAv7OiE0fqao7xznEiLKap/yQ2zXuMdpxTgT+LoN7F2ZdDfwZg8B2A/Bh4NoptustVbU7yanAvUm+AfzFiHqHq78+D2yuqr9qxQfdX6PeYkTZgee5JH9TCxj72El+DpgBfmqo+DW/06p6ctT+S9Cu/wHcUlUvJvkgg29J546571K2a9Ym4I6qenmobKn6axxL+ve17IN+Vb3nEA8x16Mfvsfga9ORbbT2mkdCHGy7kjyT5MSq2tOC1N55DvUzwBeq6odDx97TFl9M8lvAv5tmu1r6hKp6Ksn9wBnA73OY+yvJG4AvAr/UvvbOHvug+2uEBR8VMlRnV5IjgTcy+Lo+zr5L2S6SvIfBB+lPVdWLs+Vz/E4nEcTGebTKnw+t/jdg9qk1u4B3HbDv/RNo01jtGrIJuHK4YAn7axxztX0i/dVDemfkox9qcGXkPgb5dIDNwDjfHMaxtR1vnOO+JpfYAt9sHv0SYORV/qVoV5I1s+mRJMcB7wQeO9z91X53X2CQ6/y9A7ZNsr/GeVTIcHsvBe5t/bMV2JTB7J5TgA3AVw6hLYtqV5IzgN8A3ltVe4fKR/5Op9iuE4dW3ws83pbvBs5v7VsDnM+rv/Euabta297G4KLol4fKlrK/xrEVeH+bxXMO8Hwb2Eymv5bqCvU0foB/xODT70XgGeDuVv5m4K6hehcD32LwSf2RofJTGfyj3An8HnD0hNr1JuAeYEd7PbaVzzD438Jm660H/g/wugP2vxf4BoPg9dvAj06rXcDfb+/99fZ6xXLoL+DngB8CXxv6OX0p+mvU3wuDdNF72/Lr2/nvbP1x6tC+H2n7PQFcNOG/94Xa9aX272C2f7Yu9DudUrv+M/Boe//7gJ8Y2veft37cCVw+zXa19V8Grjtgv6Xur1sYzD77IYP4dQXwQeCDbXsY/GdTT7b3nxna95D7yztyJakjPaR3JEmNQV+SOmLQl6SOGPQlqSMGfUnqiEFfkjpi0Jekjhj0Jakj/x/ujSNWeGU/JgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#using the ms and mv sampling from mcmc, figure 3 can be created\n",
"fI = 1/ms - 1/mv\n",
"plt.hist(fI, 100, range=(-1,1))\n",
"print(fI.mean())"
]
},
{
"cell_type": "code",
"execution_count": 97,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>S0</th>\n",
" <th>L</th>\n",
" <th>ms</th>\n",
" <th>mv</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>30.000000</td>\n",
" <td>70.000000</td>\n",
" <td>1.079819</td>\n",
" <td>1.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>30.000000</td>\n",
" <td>70.000000</td>\n",
" <td>0.684094</td>\n",
" <td>1.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>30.788374</td>\n",
" <td>70.141137</td>\n",
" <td>0.737421</td>\n",
" <td>1.023541</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>31.319453</td>\n",
" <td>70.521572</td>\n",
" <td>0.737421</td>\n",
" <td>1.026640</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>31.319453</td>\n",
" <td>70.521572</td>\n",
" <td>0.737421</td>\n",
" <td>1.026640</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" S0 L ms mv\n",
"0 30.000000 70.000000 1.079819 1.000000\n",
"1 30.000000 70.000000 0.684094 1.000000\n",
"2 30.788374 70.141137 0.737421 1.023541\n",
"3 31.319453 70.521572 0.737421 1.026640\n",
"4 31.319453 70.521572 0.737421 1.026640"
]
},
"execution_count": 97,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#reshaping the sampling from mcmc to look like the inputs\n",
"S0 = S0.reshape(N,1)\n",
"L = L.reshape(N,1)\n",
"ms = ms.reshape(N,1)\n",
"mv = mv.reshape(N,1)\n",
"SL = np.concatenate((S0,L),axis=1)\n",
"msv = np.concatenate((ms,mv),axis=1)\n",
"x_mcmc = np.concatenate((SL,msv),axis=1)\n",
"x_MCMC = pd.DataFrame(data = x_mcmc, columns = ['S0', 'L', 'ms', 'mv'])\n",
"x_MCMC.head()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#using the trained gaussian process to predict a new y based on the mcmc sampling. This new y will be the posterior\n",
"y_mcmc = gp.predict(x_mcmc, return_std=True)\n",
"#Reverse PCA on y_mcmc to look like y\n",
"Y_mcmc = np.dot(y_mcmc, pca.components_) + pca.mean_\n",
"Y_mcmcdf= pd.DataFrame(data = Y_mcmc, columns = features)\n",
"Y_mcmcdf.head()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#posterior mean\n",
"Y_mcmc_mean = Y_mcmcdf[features].mean()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#here figure 1 is created\n",
"#plotting the Y_mcmc based on energy E_c.m. for example the column SN112EK42.5 represent energy 42.5 \n",
"Yy_mcmc = np.array([Y_mcmc[:,0],Y_mcmc[:,1],Y_mcmc[:,2],Y_mcmc[:,3],Y_mcmc[:,4],Y_mcmc[:,5],Y_mcmc[:,6],Y_mcmc[:,7],Y_mcmc[:,8],Y_mcmc[:,9],Y_mcmc[:,10],Y_mcmc[:,11],Y_mcmc[:,12],Y_mcmc[:,13],Y_mcmc[:,14],Y_mcmc[:,15],Y_mcmc[:,16],Y_mcmc[:,17],Y_mcmc[:,18],Y_mcmc[:,19],Y_mcmc[:,20],Y_mcmc[:,21],Y_mcmc[:,22],Y_mcmc[:,23],Y_mcmc[:,24],Y_mcmc[:,25],Y_mcmc[:,26],Y_mcmc[:,27],Y_mcmc[:,28],Y_mcmc[:,29],Y_mcmc[:,30],Y_mcmc[:,31],Y_mcmc[:,32]])\n",
"plt.plot(features,Yy_mcmc, 'g')\n",
"plt.plot(features, Y_mcmc_mean, 'k')\n",
"plt.plot(features,Yy, 'r')\n",
"plt.errorbar(features,Y_ex, yerr = y_ex_std)\n",
"plt.xlabel('$E_{c.m.}$ (MeV)')\n",
"plt.ylabel('$R_{n/p}$ & $DR_{n/p}$')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Reference in New Issue View Git Blame Copy Permalink