{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Causal Structures - applied on the California School data set"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The imports"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pydataset\n",
"import pylab as pl\n",
"\n",
"from halerium import CausalStructure\n",
"from halerium import Evaluator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The California Test Score Data Set"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example we apply the `CausalStructure` class to the California Test Score Data Set. The preferred data format to use with `CausalStructure` is the pandas `DataFrame`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
distcod
\n",
"
county
\n",
"
district
\n",
"
grspan
\n",
"
enrltot
\n",
"
teachers
\n",
"
calwpct
\n",
"
mealpct
\n",
"
computer
\n",
"
testscr
\n",
"
compstu
\n",
"
expnstu
\n",
"
str
\n",
"
avginc
\n",
"
elpct
\n",
"
readscr
\n",
"
mathscr
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1
\n",
"
75119
\n",
"
Alameda
\n",
"
Sunol Glen Unified
\n",
"
KK-08
\n",
"
195
\n",
"
10.900000
\n",
"
0.510200
\n",
"
2.040800
\n",
"
67
\n",
"
690.799988
\n",
"
0.343590
\n",
"
6384.911133
\n",
"
17.889910
\n",
"
22.690001
\n",
"
0.000000
\n",
"
691.599976
\n",
"
690.000000
\n",
"
\n",
"
\n",
"
2
\n",
"
61499
\n",
"
Butte
\n",
"
Manzanita Elementary
\n",
"
KK-08
\n",
"
240
\n",
"
11.150000
\n",
"
15.416700
\n",
"
47.916698
\n",
"
101
\n",
"
661.200012
\n",
"
0.420833
\n",
"
5099.380859
\n",
"
21.524664
\n",
"
9.824000
\n",
"
4.583333
\n",
"
660.500000
\n",
"
661.900024
\n",
"
\n",
"
\n",
"
3
\n",
"
61549
\n",
"
Butte
\n",
"
Thermalito Union Elementary
\n",
"
KK-08
\n",
"
1550
\n",
"
82.900002
\n",
"
55.032299
\n",
"
76.322601
\n",
"
169
\n",
"
643.599976
\n",
"
0.109032
\n",
"
5501.954590
\n",
"
18.697226
\n",
"
8.978000
\n",
"
30.000002
\n",
"
636.299988
\n",
"
650.900024
\n",
"
\n",
"
\n",
"
4
\n",
"
61457
\n",
"
Butte
\n",
"
Golden Feather Union Elementary
\n",
"
KK-08
\n",
"
243
\n",
"
14.000000
\n",
"
36.475399
\n",
"
77.049202
\n",
"
85
\n",
"
647.700012
\n",
"
0.349794
\n",
"
7101.831055
\n",
"
17.357143
\n",
"
8.978000
\n",
"
0.000000
\n",
"
651.900024
\n",
"
643.500000
\n",
"
\n",
"
\n",
"
5
\n",
"
61523
\n",
"
Butte
\n",
"
Palermo Union Elementary
\n",
"
KK-08
\n",
"
1335
\n",
"
71.500000
\n",
"
33.108601
\n",
"
78.427002
\n",
"
171
\n",
"
640.849976
\n",
"
0.128090
\n",
"
5235.987793
\n",
"
18.671329
\n",
"
9.080333
\n",
"
13.857677
\n",
"
641.799988
\n",
"
639.900024
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
416
\n",
"
68957
\n",
"
San Mateo
\n",
"
Las Lomitas Elementary
\n",
"
KK-08
\n",
"
984
\n",
"
59.730000
\n",
"
0.101600
\n",
"
3.556900
\n",
"
195
\n",
"
704.300049
\n",
"
0.198171
\n",
"
7290.338867
\n",
"
16.474134
\n",
"
28.716999
\n",
"
5.995935
\n",
"
700.900024
\n",
"
707.700012
\n",
"
\n",
"
\n",
"
417
\n",
"
69518
\n",
"
Santa Clara
\n",
"
Los Altos Elementary
\n",
"
KK-08
\n",
"
3724
\n",
"
208.479996
\n",
"
1.074100
\n",
"
1.503800
\n",
"
721
\n",
"
706.750000
\n",
"
0.193609
\n",
"
5741.462891
\n",
"
17.862625
\n",
"
41.734108
\n",
"
4.726101
\n",
"
704.000000
\n",
"
709.500000
\n",
"
\n",
"
\n",
"
418
\n",
"
72611
\n",
"
Ventura
\n",
"
Somis Union Elementary
\n",
"
KK-08
\n",
"
441
\n",
"
20.150000
\n",
"
3.563500
\n",
"
37.193802
\n",
"
45
\n",
"
645.000000
\n",
"
0.102041
\n",
"
4402.831543
\n",
"
21.885857
\n",
"
23.733000
\n",
"
24.263039
\n",
"
648.299988
\n",
"
641.700012
\n",
"
\n",
"
\n",
"
419
\n",
"
72744
\n",
"
Yuba
\n",
"
Plumas Elementary
\n",
"
KK-08
\n",
"
101
\n",
"
5.000000
\n",
"
11.881200
\n",
"
59.405899
\n",
"
14
\n",
"
672.200012
\n",
"
0.138614
\n",
"
4776.336426
\n",
"
20.200001
\n",
"
9.952000
\n",
"
2.970297
\n",
"
667.900024
\n",
"
676.500000
\n",
"
\n",
"
\n",
"
420
\n",
"
72751
\n",
"
Yuba
\n",
"
Wheatland Elementary
\n",
"
KK-08
\n",
"
1778
\n",
"
93.400002
\n",
"
6.923500
\n",
"
47.571201
\n",
"
313
\n",
"
655.750000
\n",
"
0.176041
\n",
"
5993.392578
\n",
"
19.036402
\n",
"
12.502000
\n",
"
5.005624
\n",
"
660.500000
\n",
"
651.000000
\n",
"
\n",
" \n",
"
\n",
"
420 rows × 17 columns
\n",
"
"
],
"text/plain": [
" distcod county district grspan enrltot \\\n",
"1 75119 Alameda Sunol Glen Unified KK-08 195 \n",
"2 61499 Butte Manzanita Elementary KK-08 240 \n",
"3 61549 Butte Thermalito Union Elementary KK-08 1550 \n",
"4 61457 Butte Golden Feather Union Elementary KK-08 243 \n",
"5 61523 Butte Palermo Union Elementary KK-08 1335 \n",
".. ... ... ... ... ... \n",
"416 68957 San Mateo Las Lomitas Elementary KK-08 984 \n",
"417 69518 Santa Clara Los Altos Elementary KK-08 3724 \n",
"418 72611 Ventura Somis Union Elementary KK-08 441 \n",
"419 72744 Yuba Plumas Elementary KK-08 101 \n",
"420 72751 Yuba Wheatland Elementary KK-08 1778 \n",
"\n",
" teachers calwpct mealpct computer testscr compstu \\\n",
"1 10.900000 0.510200 2.040800 67 690.799988 0.343590 \n",
"2 11.150000 15.416700 47.916698 101 661.200012 0.420833 \n",
"3 82.900002 55.032299 76.322601 169 643.599976 0.109032 \n",
"4 14.000000 36.475399 77.049202 85 647.700012 0.349794 \n",
"5 71.500000 33.108601 78.427002 171 640.849976 0.128090 \n",
".. ... ... ... ... ... ... \n",
"416 59.730000 0.101600 3.556900 195 704.300049 0.198171 \n",
"417 208.479996 1.074100 1.503800 721 706.750000 0.193609 \n",
"418 20.150000 3.563500 37.193802 45 645.000000 0.102041 \n",
"419 5.000000 11.881200 59.405899 14 672.200012 0.138614 \n",
"420 93.400002 6.923500 47.571201 313 655.750000 0.176041 \n",
"\n",
" expnstu str avginc elpct readscr mathscr \n",
"1 6384.911133 17.889910 22.690001 0.000000 691.599976 690.000000 \n",
"2 5099.380859 21.524664 9.824000 4.583333 660.500000 661.900024 \n",
"3 5501.954590 18.697226 8.978000 30.000002 636.299988 650.900024 \n",
"4 7101.831055 17.357143 8.978000 0.000000 651.900024 643.500000 \n",
"5 5235.987793 18.671329 9.080333 13.857677 641.799988 639.900024 \n",
".. ... ... ... ... ... ... \n",
"416 7290.338867 16.474134 28.716999 5.995935 700.900024 707.700012 \n",
"417 5741.462891 17.862625 41.734108 4.726101 704.000000 709.500000 \n",
"418 4402.831543 21.885857 23.733000 24.263039 648.299988 641.700012 \n",
"419 4776.336426 20.200001 9.952000 2.970297 667.900024 676.500000 \n",
"420 5993.392578 19.036402 12.502000 5.005624 660.500000 651.000000 \n",
"\n",
"[420 rows x 17 columns]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data = pydataset.data(\"Caschool\")\n",
"data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The data set relates average test scores in California schools with various data about the schools, such as the amount of students and teachers etc."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Before we start, we split the data into a training and a test set."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"np.random.seed(123)\n",
"random_indices = np.random.choice([True, False], size=len(data), p=[0.75, 0.25])\n",
"data_train = data.iloc[random_indices]\n",
"data_test = data.iloc[~random_indices]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simple structure - only inputs and outputs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The columns all represent different data about schools in california. We want to predict the average reading score 'readscr', average math score 'mathscr' and average score 'testscr'.\n",
"We do not care for now about the details of the other columns and simply treat all other numerical columns as inputs."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"outputs = {'readscr', 'mathscr', 'testscr'}"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'avginc',\n",
" 'calwpct',\n",
" 'compstu',\n",
" 'computer',\n",
" 'elpct',\n",
" 'enrltot',\n",
" 'expnstu',\n",
" 'mealpct',\n",
" 'str',\n",
" 'teachers'}"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"inputs = set(data.columns[4:]) - outputs\n",
"inputs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we define a causal structure."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"causal_structure_1 = CausalStructure([[inputs, outputs]])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We provided the training data as scaling data. The CausalStructure will use these in order to apply the correct locations and scales when building the Graph. Alternatively, we could have standardized the data."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we simply train the causal_structure with the training data and test the predictive power on the test data."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"causal_structure_1.train(data_train, method=\"MGVI\")"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"pl.figure(figsize=(10,10))\n",
"for i, col in [(1, \"mathscr\"), (2, \"readscr\"), (3, \"testscr\")]:\n",
" ax = pl.subplot(2,2,i)\n",
" ax.scatter(data_test[col], test_data_predictions_1[col])\n",
" minval, maxval = data_test[col].min(), data_test[col].max()\n",
" ax.plot([minval, maxval], [minval, maxval], ls=\"--\", c=\"k\")\n",
" ax.set_xlabel(\"real value\")\n",
" ax.set_ylabel(\"prediction\")\n",
" ax.set_title(col)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that the predictions seem to follow the real test values rather ok."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Apart from making predictions we can also look at objectives. For example we can evaluate the performance using the `Evaluator`."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'teachers': None,\n",
" 'compstu': None,\n",
" 'mealpct': None,\n",
" 'calwpct': None,\n",
" 'str': None,\n",
" 'avginc': None,\n",
" 'readscr': 0.8020755684153762,\n",
" 'enrltot': None,\n",
" 'mathscr': 0.6574450121854509,\n",
" 'testscr': 0.7577271442716229,\n",
" 'computer': None,\n",
" 'expnstu': None,\n",
" 'elpct': None}"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"performance_1 = causal_structure_1.evaluate_objective(Evaluator, data=data_test,\n",
" inputs=inputs, outputs=outputs,\n",
" metric=\"r2\")\n",
"performance_1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the R2-score for the inputs columns is returned as `None`. The values for the outputs match with the visual impression of a decent fit."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To judge whether this makes sense we have to know what these columns mean."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Caschool\n",
"\n",
"PyDataset Documentation (adopted from R Documentation. The displayed examples are in R)\n",
"\n",
"## The California Test Score Data Set\n",
"\n",
"### Description\n",
"\n",
"a cross-section from 1998-1999\n",
"\n",
"_number of observations_ : 420\n",
"\n",
"_observation_ : schools\n",
"\n",
"_country_ : United States\n",
"\n",
"### Usage\n",
"\n",
" data(Caschool)\n",
"\n",
"### Format\n",
"\n",
"A dataframe containing :\n",
"\n",
"distcod\n",
"\n",
"disctric code\n",
"\n",
"county\n",
"\n",
"county\n",
"\n",
"district\n",
"\n",
"district\n",
"\n",
"grspan\n",
"\n",
"grade span of district\n",
"\n",
"enrltot\n",
"\n",
"total enrollment\n",
"\n",
"teachers\n",
"\n",
"number of teachers\n",
"\n",
"calwpct\n",
"\n",
"percent qualifying for CalWorks\n",
"\n",
"mealpct\n",
"\n",
"percent qualifying for reduced-price lunch\n",
"\n",
"computer\n",
"\n",
"number of computers\n",
"\n",
"testscr\n",
"\n",
"average test score (read.scr+math.scr)/2\n",
"\n",
"compstu\n",
"\n",
"computer per student\n",
"\n",
"expnstu\n",
"\n",
"expenditure per student\n",
"\n",
"str\n",
"\n",
"student teacher ratio\n",
"\n",
"avginc\n",
"\n",
"district average income\n",
"\n",
"elpct\n",
"\n",
"percent of English learners\n",
"\n",
"readscr\n",
"\n",
"average reading score\n",
"\n",
"mathscr\n",
"\n",
"average math score\n",
"\n",
"### Source\n",
"\n",
"California Department of Education http://www.cde.ca.gov.\n",
"\n",
"### References\n",
"\n",
"Stock, James H. and Mark W. Watson (2003) _Introduction to Econometrics_,\n",
"Addison-Wesley Educational Publishers,\n",
"http://wps.aw.com/aw_stockwatsn_economtrcs_1, chapter 4–7.\n",
"\n",
"### See Also\n",
"\n",
"`Index.Source`, `Index.Economics`, `Index.Econometrics`, `Index.Observations`\n",
"\n",
"\n"
]
}
],
"source": [
"pydataset.data(\"Caschool\", show_doc=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that in the fully connected model (all inputs influence all outputs) the biggest influences seem to be the total enrollment 'enrltot', the amount of teachers 'teachers' adn the percentage of students that get subsidized meals 'mealpct'. The influence of the other tests on 'testscr' is zero.\n",
"\n",
"We might now ask the question whether it makes sense that these quantities are the direct main influences or if they are only confounders. If we take a look at the correlations of the inputs, we see that the inputs have significant correlations that can lead to confounding and the explain-away-effect."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
teachers
compstu
calwpct
mealpct
str
avginc
enrltot
computer
expnstu
elpct
\n",
"
\n",
"
teachers
\n",
"
1.000000
\n",
"
-0.205573
\n",
"
0.105457
\n",
"
0.142855
\n",
"
0.272616
\n",
"
0.023591
\n",
"
0.997424
\n",
"
0.936511
\n",
"
-0.082112
\n",
"
0.357818
\n",
"
\n",
"
\n",
"
compstu
\n",
"
-0.205573
\n",
"
1.000000
\n",
"
-0.171575
\n",
"
-0.217232
\n",
"
-0.301937
\n",
"
0.183565
\n",
"
-0.213354
\n",
"
-0.041737
\n",
"
0.274414
\n",
"
-0.272306
\n",
"
\n",
"
\n",
"
calwpct
\n",
"
0.105457
\n",
"
-0.171575
\n",
"
1.000000
\n",
"
0.735591
\n",
"
0.017576
\n",
"
-0.511519
\n",
"
0.099840
\n",
"
0.072180
\n",
"
0.081346
\n",
"
0.350795
\n",
"
\n",
"
\n",
"
mealpct
\n",
"
0.142855
\n",
"
-0.217232
\n",
"
0.735591
\n",
"
1.000000
\n",
"
0.155836
\n",
"
-0.708045
\n",
"
0.147787
\n",
"
0.076777
\n",
"
-0.031435
\n",
"
0.679559
\n",
"
\n",
"
\n",
"
str
\n",
"
0.272616
\n",
"
-0.301937
\n",
"
0.017576
\n",
"
0.155836
\n",
"
1.000000
\n",
"
-0.222932
\n",
"
0.305878
\n",
"
0.243671
\n",
"
-0.590435
\n",
"
0.226483
\n",
"
\n",
"
\n",
"
avginc
\n",
"
0.023591
\n",
"
0.183565
\n",
"
-0.511519
\n",
"
-0.708045
\n",
"
-0.222932
\n",
"
1.000000
\n",
"
0.010920
\n",
"
0.076274
\n",
"
0.298755
\n",
"
-0.362238
\n",
"
\n",
"
\n",
"
enrltot
\n",
"
0.997424
\n",
"
-0.213354
\n",
"
0.099840
\n",
"
0.147787
\n",
"
0.305878
\n",
"
0.010920
\n",
"
1.000000
\n",
"
0.929746
\n",
"
-0.100289
\n",
"
0.365512
\n",
"
\n",
"
\n",
"
computer
\n",
"
0.936511
\n",
"
-0.041737
\n",
"
0.072180
\n",
"
0.076777
\n",
"
0.243671
\n",
"
0.076274
\n",
"
0.929746
\n",
"
1.000000
\n",
"
-0.061414
\n",
"
0.292201
\n",
"
\n",
"
\n",
"
expnstu
\n",
"
-0.082112
\n",
"
0.274414
\n",
"
0.081346
\n",
"
-0.031435
\n",
"
-0.590435
\n",
"
0.298755
\n",
"
-0.100289
\n",
"
-0.061414
\n",
"
1.000000
\n",
"
-0.049224
\n",
"
\n",
"
\n",
"
elpct
\n",
"
0.357818
\n",
"
-0.272306
\n",
"
0.350795
\n",
"
0.679559
\n",
"
0.226483
\n",
"
-0.362238
\n",
"
0.365512
\n",
"
0.292201
\n",
"
-0.049224
\n",
"
1.000000
\n",
"
\n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_train[inputs].corr().style.background_gradient(cmap='coolwarm', vmin=-1, vmax=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see for example that the total enrollment 'enrltot' and the amount of teachers 'teachers' are very strongly correlated (also with the amount of computers)."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
teachers
compstu
calwpct
mealpct
str
avginc
enrltot
computer
expnstu
elpct
\n",
"
\n",
"
testscr
\n",
"
-0.157794
\n",
"
0.278550
\n",
"
-0.628098
\n",
"
-0.877528
\n",
"
-0.227190
\n",
"
0.734062
\n",
"
-0.167217
\n",
"
-0.084573
\n",
"
0.150482
\n",
"
-0.662471
\n",
"
\n",
"
\n",
"
readscr
\n",
"
-0.189992
\n",
"
0.285396
\n",
"
-0.615917
\n",
"
-0.887392
\n",
"
-0.252007
\n",
"
0.720649
\n",
"
-0.199871
\n",
"
-0.117787
\n",
"
0.178198
\n",
"
-0.704534
\n",
"
\n",
"
\n",
"
mathscr
\n",
"
-0.116787
\n",
"
0.260359
\n",
"
-0.617021
\n",
"
-0.832906
\n",
"
-0.191503
\n",
"
0.720223
\n",
"
-0.125351
\n",
"
-0.045295
\n",
"
0.114616
\n",
"
-0.591257
\n",
"
\n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_train.corr().loc[inputs, outputs].T.style.background_gradient(cmap='coolwarm', vmin=-1, vmax=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Advanced structure"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the information we have about the data set we can now try to make a more realistic model of the involved columns.\n",
"\n",
"We start with a number of assumptions which result in dependencies."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"dependencies = []"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First of all we see in the data documentation that 'testscr' is just the average test score."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"dependencies += [\n",
" [[\"mathscr\", \"readscr\"], \"testscr\"]\n",
"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Assumption 1: We believe the student-teacher-ratio is what actually matters for the test performance. The student teacher ratio of course depends on the number of students and the number of teachers."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"dependencies += [\n",
" [\"str\", [\"mathscr\", \"readscr\"]],\n",
" [[\"teachers\", \"enrltot\"], \"str\"]\n",
"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Assumption 2: We believe the amount of teachers depends on the number of students and the expenses per student."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"dependencies += [\n",
" [[\"expnstu\", \"enrltot\"], \"teachers\"]\n",
"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Assumption 3: The amount of computers per student might indicate the learning possibilities at the school and therefore influence the test results. The amount of computer per student of course is determined by the number of students and the number of computers."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"dependencies += [\n",
" [\"compstu\", [\"mathscr\", \"readscr\"]],\n",
" [[\"computer\", \"enrltot\"], \"compstu\"]\n",
"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Assumption 4: The amount of computers probably depends on the amount of students and the funding per student."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"dependencies += [\n",
" [[\"expnstu\", \"enrltot\"], \"computer\"]\n",
"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Assumption 5: We assume poorer students achieve worse test results on average. In the data we have the percent qualifying for reduced-price lunch 'mealpct', the percent qualifying for CalWorks 'calwpct' and the average income in the school district 'avginc'. Since the latter is a bit indirect we assume 'mealpct' and 'calwpct' are direct influences on the test results while both of them are influenced by 'avginc'."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"dependencies += [\n",
" [[\"mealpct\", \"calwpct\"], [\"mathscr\", \"readscr\"]],\n",
" [\"avginc\", [\"mealpct\", \"calwpct\"]]\n",
"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Assumption 6: We assume that being a non-native speaker influences the test results directly as well as the economic situation of the student. The percent of English learners is 'elpct'."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"dependencies += [\n",
" [\"elpct\", [\"mathscr\", \"readscr\"]],\n",
" [\"elpct\", [\"mealpct\", \"calwpct\"]]\n",
"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we put in these assumptions into a causal structure."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"causal_structure_2 = CausalStructure(dependencies)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and train it"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"causal_structure_2.train(data_train, method=\"MGVI\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us have a look at the performance on the test data."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"test_data_inputs = data_test[inputs]\n",
"test_data_predictions_2 = causal_structure_2.predict(test_data_inputs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First visually..."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmwAAAJcCAYAAABE7/iIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAACUy0lEQVR4nOzdeZyO9f7H8ddnxmCso9Ay2SpRUpS2o8hSiCIp6rQvquO0qBTtnRKl06nOrz2V0yKRRE6ULJ2jRUSL4lQiRokYimG27++P+76ne2bubWbudeb9fDw8zFz3dV/358J8fa7v9jHnHCIiIiKSvNISHYCIiIiIhKaETURERCTJKWETERERSXJK2ERERESSnBI2ERERkSSnhE1EREQkySlhk5RkZmvNrHei4xARqSwzc2Z2cKLjkNSghE2Snpm9aGb3JToOERGRRFHCJhKAmdVKdAwiklySuV0wD/2fXo3pL1dixjtsOcrMvjCznWY20cz2MbN3zOw3M5tnZk285041s5/NbLuZfWBmHbzHhwN/Bm42s9/NbJbfR3TyXnu7mU0xs7re9zQ1s7fNLNfMtprZf3wNmZm1MLPpZrbZzH41s//zHr/YzBab2T/MbCtwdxz/qEQkSXnbsVvM7Atgp5mdaGYfetuXz83sZL9zLzGzb7zt2xozu7LMtUaZ2U9mttHMLi3z2mlm9rX3vTlmdpPfawPNbIWZ7TCz782sr/f4QjMba2aLgV3AgTH8o5AEU8ImsXYWcApwCHA68A5wK9AUz7+/a73nvQO0BZoDnwGvADjnnvF+/aBzroFz7nS/a58D9AXaAEcAF3uP3whsAJoB+3g/z5lZOvA2sA5oDWQDr/ld7zhgjTeGsVG4dxGpHs4F+uNJiN4C7gP2Am4C3jCzZt7zfgEGAI2AS4B/mNlRAN4k6yY87WFboOwc3InAlc65hsDhwHzv+44F/gWMArKAbsBav/ddAAwHGuJp26SaUsImsfZP59wm51wO8B/gE+fccufcHuBNoDOAc+5559xv3uN3A0eaWeMw137MObfRObcVmAV08h4vAPYDWjnnCpxz/3GeornHAvsDo5xzO51zu51z//W73kbn3D+dc4XOubzo3L6IVAOPOefWA+cD/3bO/ds5V+ycew9YCpwG4Jyb7Zz73nksAt4FTvJe4xzgBefcV865nZTvxS8ADjOzRs65bc65z7zHLwOed8695/3MHOfcKr/3veicW+lttwpicveSFJSwSaxt8vs6L8D3Dcws3czGe7v6d/DH02PTMNf+2e/rXUAD79cTgO+Ad73DEqO9x1sA65xzhUGutz7M54lIzeRrG1oBZ3uHQ3PNLBc4Ec8DImbWz8w+9k7FyMWTyPnasf0p3caU7Q07y3v+OjNbZGYneI+3AL6PIDap5pJ2AqXUKOcBA/EMEawFGgPbAPO+7ipyMefcb3iGRW/0zoVbYGaf4mnYWppZrSBJW4U+R0RqDF/bsB54yTl3RdkTzKwO8AZwIfCWc67AzGbwRzv2E57ky6dlqQ9w7lNgoJllAH8FXveevx44KILYpJpTD5skg4bAHuBXoB5wf5nXN1GBybRmNsDMDjYzA3YARd5fS/A0muPNrL6Z1TWzrtG4ARGpEV4GTjezPt6RgbpmdrKZHQDUBuoAm4FCM+sHnOr33teBi83sMDOrB9zle8HMapvZn82ssXdY09dugWdu2yVm1svM0sws28zax+FeJckoYZNk8C88wwM5wNfAx2Ven4hnbkeu94k1nLbAPOB34CPgCefcQudcEZ6FDwcDP+JZmDA0KncgItWedx7bQDwLmTbj6f0aBaR5e/avxZOYbcMzcjDT773vAI/gWUzwnfd3fxcAa73TQq7CM18O59wSvAsYgO3AIjxDs1LDmGcutoiIiIgkK/WwiYiIiCQ5JWwiIiIiSU4Jm4hIGWbWzruzvO/XDjO73sz2MrP3zOxb7+9N/N4zxsy+M7PVZtYnkfGLSPWjOWwiIiF4K2Tk4KmEMQLY6pwb793fr4lz7hYzOwyYzB+bM88DDvEudBERqbKU3oetadOmrnXr1okOQ0TiaNmyZVucc83Cnxk1vYDvnXPrzGwgcLL3+CRgIXALnpWDr3krdfxgZt/hSd4+CnVhtWEiNUtV2q+UTthat27N0qVLEx2GiMSRmcW7XuIwPL1nAPs4534CcM79ZGbNvcezKb0dzQbvsXLMbDie2o+0bNlSbZhIDVKV9ktz2EREgjCz2sAZwNRwpwY4FnC+iXPuGedcF+dcl2bN4tlRKCKpTAmbiEhw/YDPnHO+GribzMxXN3I/4Bfv8Q2ULjt0ALAxblGKSLWnhE1EJLhz+WM4FDw711/k/foi4C2/48PMrI6ZtcFTbWNJ3KIUkWovpeewiYjEirfe4ynAlX6HxwOvm9lleMqbnQ3gnFtpZq/jKa1WCIzQClERiSYlbCIiATjndgF7lzn2K55Vo4HOHwuMjUNoIlIDaUhUREREJMkpYRMRERFJckrYRERERJKcEjYRSZiioiJUHk9EUlVRUfzWFilhE5GE2LNnD4MHD+aee+5JdCgiIhX2r3/9ixNOOIFt27bF5fOUsIlI3OXl5XHmmWcyc+ZMmjdvHv4NIiJJZOLEiVx88cU0atSI2rVrx+Uzta2HiMTd1VdfzZw5c3j22We5/PLLEx2OiEjE3nvvPS6//HL69u3L9OnTyczMjMvnKmETkbi7/fbb6du3L8OGDUt0KCIiFdKjRw/+/ve/M2LECOrUqRO3z9WQqIjExY4dO/jHP/6Bc46DDz5YyZqIpJRnnnmGjRs3UqtWLW644Ya4JmughE1E4mDbtm2ccsop3HzzzXz++eeJDkdEpELuvfderrzySh577LGExRCzhM3M2pnZCr9fO8zsejPby8zeM7Nvvb838XvPGDP7zsxWm1mfWMUmIvHz66+/0qtXL1asWMEbb7xBp06dEh2SiEhEnHPccccd3HnnnVxwwQWMHZu46nMxS9icc6udc52cc52Ao4FdwJvAaOB951xb4H3v95jZYcAwoAPQF3jCzNJjFZ+IxN4vv/xCjx49+Prrr3nrrbc444wzEh2SiEhEnHOMHj2a++67j8suu4wXXniB9PTEpSXxGhLtBXzvnFsHDAQmeY9PAgZ5vx4IvOac2+Oc+wH4Djg2TvGJSAx8+eWXrF+/ntmzZ9O3b99EhyMiErGdO3fyzjvvcPXVV/PMM88kNFmD+K0SHQZM9n69j3PuJwDn3E9m5tuEKRv42O89G7zHSjGz4cBwgJYtW8YsYBGpvD179lCnTh169erFDz/8QFZWVqJDEhGJSHFxMUVFRTRo0ID//Oc/NGrUCDNLdFix72Ezs9rAGcDUcKcGOFauZo1z7hnnXBfnXJdmzZpFI0QRiaJ169bRsWNHXnvtNQAlayKSMoqKirjiiisYOnQoRUVFNG7cOCmSNYjPkGg/4DPn3Cbv95vMbD8A7++/eI9vAFr4ve8AYGMc4hORKFmzZg3dunVj8+bNHHjggYkOR0QkYoWFhVxyySU8//zzdOzYkbS05NpIIx7RnMsfw6EAM4GLvF9fBLzld3yYmdUxszZAW2BJHOITkSj43//+R7du3fj99995//33OfZYTUEVkdRQUFDA+eefz0svvcS9997LPffckzQ9az4xncNmZvWAU4Ar/Q6PB143s8uAH4GzAZxzK83sdeBroBAY4ZwrimV8IhIdW7ZsoXv37hQVFbFgwQKOOOKIRIckIhKxK6+8kilTpvDggw8yatSoRIcTUEwTNufcLmDvMsd+xbNqNND5Y4HEbXIiIpXStGlTRo4cyYABAzjssMMSHY6ISIVcddVVHH300YwYMSLRoQSlWqIiUmnLly8nLS2NI488kptvvjnR4YiIRCwvL4+ZM2cydOhQjj322KSfxpFcM+pEJGUsWbKEnj17ctlll+FcuQXdIiJJa+fOnQwYMIBzzz2Xr776KtHhREQ9bCJSYR9++CF9+/aladOmvPHGG0k3OVdEaqYZy3OYMHc1G3Pz2D8rk1F92jGoc+ktXX/77Tf69+/P4sWLmTRpEocffnjMPiualLCJSIUsWrSI/v37s//++zN//nwOOOCARIckIsKM5TmMmf4leQWe9Yo5uXmMmf4lQEkitX37dvr168eSJUt49dVXGTp0aMw+K9o0JCoiFfLYY4/RsmVLFi1apGRNRJLGhLmrSxIon7yCIibMXV3y/fz581m2bBmvv/56pZO1SD8r2tTDJiIRcc5hZrz88sv8/vvvqNKIiCSTjbl5QY/72q8zzzyT7777jhYtWgQ8NxqfFSvqYRORsGbNmsWJJ55Ibm4umZmZStZEJOnsn5UZ8HjT9DyOO+445s+fD1DlZC3UZwU7Hg1K2EQkpDfeeIPBgwdTUFCg1aAikrRG9WlHZkZ6qWMZu7fz06tj+Oqrr6LafgX6rMyMdEb1aRe1zyhLQ6IiEtRrr73G+eefz3HHHce///1vGjdunOiQREQC8k32963c3Nt+Z+Obt/Pb1s3MmTOHbt26xeyztEpURBJm2rRp/PnPf+bEE0/k7bffpmHDhokOSUQkpEGdsxnUOZtNmzZxwgknsDP3V959911OOOGESl0v1NYdvs+KFw2JikhAXbp04YILLuDf//63kjURSSnNmjWjf//+zJs3r0rJ2pjpX5KTm4fjj607ZizPiW6wEVLCJiKlvP/++xQXF9O6dWtefPFF6tevn+iQREQisnr1atatW0daWhr//Oc/OeaYYyp9rURs3RGKEjYRKfHII4/Qu3dvnn766USHIiIS0ozlOXQdP582o2fTdfx8Hps2n+7du3PuuedGZYFBIrbuCEUJm4gA8OCDDzJy5EjOOussLrvsskSHIyISVNnhyh/+t5IbLjqT/CKYOHFiVMrlJWLrjlCUsIkI9957L7fccgvDhg3jtddeo3bt2okOSUQkKP/hyj0/f8emybdBem1aX/Qghx56aFQ+IxFbd4SihE2khluzZg1jx47lwgsv5OWXX6ZWLS0eF5Hk5j8smbvweax2JvucN56ttfaO2qKAQZ2zGTe4I9lZmRiQnZXJuMEd47oy1J9aZpEa7sADD+STTz6hY8eOpKXpGU5Ekt/+WZnkeJO2pgNH4wp2U6tRc4CoFmGP99Ydoah1FqmBnHOMHDmS559/HoAjjzxSyZqIpIx+e//K1pkP4AoLSM9sVJKsQWJXcsaSWmiRGqa4uJgRI0bwyCOPsHLlykSHIyJSIe+99x73X3cRWXs2UZy/K+A5iVrJGUtK2ERqkKKiIoYPH86TTz7JLbfcwkMPPZTokEREIvbvf/+b008/nYMPPpilH/2HlvvvG/C8RK3kjCUlbCI1hHOOSy+9lIkTJ3LnnXcybty4qCx9FxGJh7fffptBgwbRoUMHFixYQPPmzZNuJWcsKWETqSHMjEMOOYR7772Xe+65R8maiKSUb3+vTf02ndh80s2c8eznzFiek3QrOWNJq0RFqrn8/HzWrFlD+/btue222xIdTsowsyzgOeBwwAGXAnnAU0BdoBD4i3Nuiff8McBlQBFwrXNubgLCFql2vvzyS74v3ItnVhbT+My7gD/qekJyreSMJfWwiVRje/bsYciQIfzpT3/i119/TXQ4qeZRYI5zrj1wJPAN8CBwj3OuE3Cn93vM7DBgGNAB6As8YWbpgS4qIpGbNGkSnTp1YtS4/0uqup6JoIRNpJrKy8tj4MCBzJo1i/vvv5+999470SGlDDNrBHQDJgI45/Kdc7l4etoaeU9rDGz0fj0QeM05t8c59wPwHXBsXIMWqWaee+45LrnkEnr16kV+iy4Bz6mOq0GDUcImUg3t3LmTAQMG8O677zJx4kSuuuqqRIeUag4ENgMvmNlyM3vOzOoD1wMTzGw98BAwxnt+NrDe7/0bvMfKMbPhZrbUzJZu3rw5Zjcgksoef/xxrrjiCvr168fMmTM5oFmTgOdVx9WgwShhE6mGJkyYwMKFC5k0aRKXXnpposNJRbWAo4AnnXOdgZ3AaOBqYKRzrgUwEm8PHBBoBYcLdGHn3DPOuS7OuS7NmjWLfuQiKW716tVce+21DBw4kOnTp1O3bt2Aq0ENz1y2ruPnhy1HNWN5Dl3Hz6fN6NkRnZ+MtOhApBoaM2YM3bt3p0ePHokOJVVtADY45z7xfj8NT8J2InCd99hUPIsSfOe38Hv/AfwxXCoiFdCuXTveffddunXrRkZGBvBHmakJc1eTk5uH8ccTUdkFCDOW5zBh7mo25uaxf1YmPdo3441lOSVz4MqenyrUwyZSTWzbto2LL76YX3/9lTp16ihZqwLn3M/AejPzbebUC/gaTxLW3XusJ/Ct9+uZwDAzq2NmbYC2wJI4hiyS0pxzjB07lnfeeQeAXr16lSRrPoM6Z7N4dE+yszLLdV/7FiDMWJ7DmOlfkpObh8OTnL3y8Y/VYsGCethEqoEtW7Zwyimn8PXXX3PBBRfQq1evCl+j7FPpqD7tUurpMwauAV4xs9rAGuAS4C3gUTOrBewGhgM451aa2et4krpCYIRzrijwZUXEn3OO22+/nfvvv5/hw4fTr1+/kOcHW2iQk5vH9VNWlL9+Ba+TrJSwiaS4TZs20bt3b7777jveeuutSidrY6Z/mfJDBtHknFsBlF2a9l/g6CDnjwXGxjgskaRVmYc+5xyjRo3i73//O1dccQVPPvlk2Gvtn5VJThSSrVRbsKAhUZEqSuRk1p9++omTTz6Z77//nrfffpu+fftW6joT5q6uFkMGIpIYgYYix0z/MmR76Jzjuuuu4+9//zsjRozgqaeeIi0tLey1Ai1ACKfsqqBULF+lhE2kCirTSEVTcXExdevWZc6cOZXqWfMJNjSQakMGIpIYlX3oy8/P54YbbuCf//wnaWlpEV2rbDmqcDIz0vnz8S1TvnyVhkRFqiBUwxLLxuCnn36iefPmZGdns2zZspKGrrKCDTGk2pCBiCRGRR76ioqK+OWXX9hvv/144oknMLNStY0juZZ/Oaqu4+cHHSLNrkbzcdXDJlIFieiZ+v777zn++OO57jrP7hJVTdYg8BBDKg4ZiIhHvKdqBHu4K3u8sLCQiy66iOOPP57c3FzS0tJKJWsVuZZPsPbrkaGdWDy6Z7VI1kAJm0iVVLRhqarVq1fTrVs3du7cyWWXXRa165YdYojHkEF12MhSJBklYqpGJA99BQUFnHfeebzyyitceeWVZGVlVfpa/hLRfiWChkRFqmBUn3alVldC+J6pym6fsXLlSnr16oVzjgULFtCxY8eo3IOP/xBDrGlVqkjsJGKqhv/GtoHatj179jBs2DBmzJjBQw89xI033ljpawV7T3VvO5SwiVRBRRuWyiYq+fn5DBgwgLS0NN5//30OPfTQKN9JfCVq7p9ITZCoRUShkqa7776bGTNm8Nhjj3HNNddU6Vrxlix7VMY0YTOzLDylWw7Hs3fdpUAe8BRQF88Gk39xzi3xnj8GuAwoAq51zs2NZXwi0VCRhqWyiUrt2rWZNGkS++23H23btq1SvMlAq1JFYicZFxGNHj2ao48+miFDhiQshspIptGAWM9hexSY45xrDxwJfAM8CNzjnOsE3On9HjM7DBgGdAD6Ak+YWcU2WhFJchVNVD755BOeffZZALp161YtkjWI/9w/kZokmouIqjLXdOfOnYwePZq8vDwaN26ccskaJNcelTHrYTOzRkA34GIA51w+kG9mDmjkPa0xfxRIHgi85pzbA/xgZt8BxwIfxSpGkXiryJPv4sWL6devH/vssw/nn38+mZnVJ5mpzNw/EYlMZeaA+fgP/zXOzGBnfiEFRZ7iThXpXfrtt9847bTT+PDDD6nb6kje275PwocUKyOZRgNiOSR6ILAZeMHMjgSWAdcB1wNzzewhPD18f/Kenw187Pf+Dd5jpZjZcLz1+1q2bBmr2EViItJEZeHChQwYMIDs7Gzmz59fkqwly1yKqqrKfygiEl5l5oCVHf7LzSsod04kUzhyc3Pp168fn376KTfc/zhTNjYhr8CT4KTaAqNkGl6OZcJWCzgKuMY594mZPQqMxtOrNtI594aZnQNMBHpTvnIEBKjZ6px7BngGoEuXLsFquookjbJJ1llHZ7Ng1eagicq8efM444wzaNOmDe+//z777rtvyXWSZS5FNCTTpGIRCTz8F0io3qWtW7dy6qmn8sUXXzBt2jQmrGpUkqz5pNICo2QaDYhlwrYB2OCc+8T7/TQ8CduJeHraAKbiWZTgO7+F3/sP4I/hUpGUFCjJemNZTsg9gr766ivatm3LvHnzaNasWclxrawUkViKdJgvVO/SL7/8wqZNm3jzzTfp378/I0fPrtJnJVoyjQbELGFzzv1sZuvNrJ1zbjXQC/gaz1Bpd2Ah0BP41vuWmcCrZvYwsD/QFlgSq/hE4qEiSdaOHTto1KgR119/PVdddRV169Yt9XoyzaUQkeon2PCfv7K9S74RhA2bfiW7+V7c3Lc93377bUn7lUxDipWVLKMBsV4leg3wipl9AXQC7geuAP5uZp97vx8O4JxbCbyOJ6mbA4xwzoXvmxVJYpEmWdOmTaNNmzYsX74coFyyBlpZKSKxFWh1aUaa0aReRsAKAr4RhHXrN7DxXzfw9eznGTP9S+Z882vJ6tKc3Lxy8520wKhyYroPm3NuBdClzOH/AkcHOX8sMDaWMYnEUyRPl6+++ioXXnghxx9/PAcddFDQayXTXAoRqX4qOvw3Ye5qfvv1Zza9ditFO3Op2/pI8gqKuGfWSnYXFJe0VQ7PJHVH9SrGHm+qdCASQ+GSrEmTJnHppZdy0kkn8fbbb9OgQYOg10qmuRQiUj0FGv4Ltjp93bq1/Dz5VorzfmOfc+6lTnZ7ALbtKr+61JesLR7dMx63US0pYROJoVBJ1rx587jkkkvo1asXb731FvXq1YvoekrQRCRegq1OL8jfw5Ypt+P27GSfYWOps1/4Tb0137ZqlLCJxFiwJKt79+7cf//9XH/99QHnrImIREMk+zcGOyfYwqlHFqzlmtF3Me3bQtzerUtey8xIp06ttIB7uGm+bdUoYROJs+eff57+/fuzzz77MHr06ESHE5HqsmGvSE3g//OaVS+D33cXUlAcvFpBqD0ey/aK5W9eR+H2n9l48HFMGH8lXQO0DYDm28aAEjaRGCqb6By48T1e/ud4br75Zh544IGI3pPo5Ki6bdgrUp2V/XkNNJ+s7NZCobYf8l84lf/LGja9djuWUZeDOnmKFIWappFM7Vh1EOttPURqLF/DmZObR7FzrJz9PC//czzd+g1i7NjAi6H93+P4IzmqSMHlaEum4sciElplqhWE2n7It9XHnp++ZdPkW7FatWl1/v3c0r9jyOsP6pzNqD7t2D8rk425eUyYuzqh7Vh1oB42kSoKN/fDOUfuf15mx0dTqH94LwpO/Au1agX+0UvGagbasFckcSra416ZagWhth8a1Dmb1V8s47YJt5NWtwFHDP87d5x7ctj2SD3z0aeETWqsaAw9RjL3w+XnsWv1Yhoc2Ye9+ozgpx35Qa+XjMlRddipXCQVVSbpqUy1gnDbD21f/TGtD9iX+fPn07Jly4hiT8aHz1SnIVGpkaI19BiqUdqvcV1ccRFpdeqx7wUPsVefEZilhUx0krGaQaDdzzWBWCT2KjMdIWC1gnQjKzNwtQKfOrX+SAea1Mtg3OCO9D+8OQBjx45lyZIltGzZsqSCQZvRs+k6fn7QNjMZHz5TnRI2qZGiNS8rWOOTs20njT6bxLa3H8IVF5FetwFmaWETnWRMjgZ1zmbc4I5kZ2WGbPBFJLoqk/QE+nmdMORIVtx1Kj+M78/i0T1L/ez6Hl79t+HYXVDMio8/4NBDD+W7777DzNhrr70q9KCbjA+fqU5DolIjRevpL9DwgysuYveCJ5m7dA6DLxnBT03q89P23RENu4arZpCoFaTasFck/io7HaEiP6+BHl63rvqYv42/nyM6HEbjxo1DnusrRVW2XVIpvehTwiZxlSxbVkRrXlbZRskVF5H7zqPs+Go+d911F3fddRdmZUsfhxassdUkXpHqr+weahlpVrKHGkQ/6Sn7kLrrfx+x+a0HqN28NfPnz2evvfYKeq7Ptl0FJduH+NqlcYM7Mm5wx6Ro76sLJWwSN8mUcETr6a9sj1jewqfZ8dV8xo4dy6233hrVmDWJV6R6C7aHWmZGGrsLimOS9Pg/vOatXcHmt8ZTe9+D6XT5g6WStbLnhuJrl8oOv0rVKGGTuElkwhGoZy9aT3/+PWIfftiEJUtO4frrr4/yHWgSr0h1F2wPtd0FxfxjaKeYtJP+D6919m9Hw8792bfnhYwZdHTIc8NRuxR9StgkbhKVcATr2Rs3uCOLR/es8vV3797NO++8w5lnnsmf/vQn/vSnP1X5moFoew2R6i1YW+iAG1//HIj+aMSgztl8uug93t/ejE1k0nHIdUEfXgPNsd25p1B1Q+NEq0QlbhK1aiiWO/Xv2rWLgQMHctZZZ7Fy5coqXy+UZFxBKiLRE6otLHIuJlVPnn76ae4feQk99ywOuIq0rEGds1k8umfJuXef0UHtUpyoh03iJlGrhirasxfpwoidO3dy+umns3DhQg4cfBMDXlrL/lmbYjaxNtwKUhFJHpVZYDWqTztGTlmBC/J6NKaQ+MdlX8/hh1n/R//+/bnzzjsrdT21S/GjhE3iJlE/2BUZSox0YcSOHTvo378/H374IfuecROFB3cPeX60aHsNkeRX2QVWgzpns3TdVl75+MegSVtVppD4x7X9k+nkLnyeBu1O4KI7/0ndunUrvYpf7VJ8KGGTuIrnD7av8cnJzcOgVAMYrGcv0oUR7777Lh9//DEHD72NPS2PC3t+NO5DT68iqaGiC6zK/oz/+fiWTP5kPUWufNpWlSkkvriKdv/Ob0tnUK/9Sew14EYemf8DGRm1k2YVvwSmOWxSLfnvyA2eZM23G1qonfrDDZ86bwM6ZMgQVq9eTX6ZZC3cdSoqWiW0RCR+glZAyc0rV84p0M/4G8tyOPe4FlGfG5azbRfOOdLrNmDfCx6i6ek3Yem12JibF9O5vhIdStikWgrU+Dg8yVqoSbXBnl7TzGh57as0btOR+ydOB+DAAw+M+UIKNaIiqSfUz3/Zh65gP+MLVm2Oakk45xyFn7zKtvnP4ZyjVqPmWFp6SbzaNij5KWGTaqmyjU+glZgA+b9t5efJY/g951ueev+bksY21is31YiKpJ5g7YiP/0NXqJ/xsisyq5Ks3XjjjeQsmkx6cQH+E0R87ZVqfyY/zWGTaqkqNfjgj4URaWbs2bGZTa/dTtFvm2k25C7SWh5ZMhfF//yc3DzSzUo1xlWd+6G910RST9l2IRBfohbrn/Hi4mIGnHsp77w+iYZHn06jXsNJSzOc8/Ta+c+JDbWKf8byHO6eubJkz7Um9TK46/QOmt8WR+phk2qpoj1fM5bn0HX8fNqMns2EuasZ1acdP4zvT8Hv29j06hiKfv+V5mffQ2arI4HST8WDOmeXfJ5vknC05ppp7zWR1OTrHcsO03MV659xX7LW6NjBNOk1HDNPspYG7MovZOSUFXQdPx8g6BDsjOU5jJr6eakNcrftKmDUtM81nzaO1MMm1VJFthAJtQT/gP2as631kTQ4vDd1stuXvKfs02+sym5pj6PEMbMs4DngcDxjSJc65z4ys2uAvwKFwGzn3M3e88cAlwFFwLXOubkJCVySSrj9Jyv7Mx7p6vFv0lrR+E/DaHzinzGzkuPFELBge6DqLxPmri5VgN6noMiplnEcKWGTaivSLUQCJVs7flnPfa9v5/ZzTmRM3nVhN/uN5Vwz7XGUMI8Cc5xzQ8ysNlDPzHoAA4EjnHN7zKw5gJkdBgwDOgD7A/PM7BDnXPiii1Lt1c1IK2lDsjIzuPuM0kOJFf0ZD7fPW2FhIZ988gldu3bFtTqWrFbHhr1mqAfMUO2Y5tPGj4ZEpcYr2+AU/LqeTa+O5quX/sbATvtHtFJLE3arFzNrBHQDJgI45/Kdc7nA1cB459we7/FfvG8ZCLzmnNvjnPsB+A4I/7+kVGu+xMrXkwWwp7C4Qu/3TdXw3w4kVI9+QUEB5557Lt27d+d///tfheINlnyFasfUxsWPetikwqrbRq7+k37zN69l05TbAWg/5AbMLKKn30SV3ZKYORDYDLxgZkcCy4DrgEOAk8xsLLAbuMk59ymQDXzs9/4N3mPlmNlwYDhAy5YtY3YDkniRTpUI1KYCQXvRgu7ztmUHZ599Nm+99RYPP/wwhxxyCFmZawMWZw8kWPI1qk87Rk39vNywaEa6qY2LI/WwSYVUx41cfZN+839Zw6bJt2KWRqsLHuTui/pGfI1BnbOjumeSJFwt4CjgSedcZ2AnMNp7vAlwPDAKeN08E4MswDUCVhdyzj3jnOvinOvSrFmzmAQvySHcVIkZy3Po/Ld3uX7KinJt6j2zVgZN9gIlVq4wn+1vj+ett97i//7v/xg5ciQAd5/RgYy08v8808scC/WAOahzNhPOPpKszIySY03qZTBhyJFq4+Io4h42M8sGWvm/xzn3QSyCkuQVq8n1ieSL+6IhY7BatTnkkgepvVc2I6esKFkxqnp6qa0S7dcGYINz7hPv99PwJGwbgOnOU/JiiZkVA029x1v4vf8AYGP07kBSUagtO8rOQ/OXV1AU8Dh4kr1/DO1U7r35q//D9v8t4emnn2b48OElx4Mtagh0LFydU7VviRVRwmZmDwBDga/xrIACz9OjErYaJpk2co3m0Oygztl89/G7TP/kOx75ZHvJEILq6aW+yrRfzrmfzWy9mbVzzq0Gennf/z3QE1hoZocAtYEtwEzgVTN7GM+ig7bAkhjdkqSIUFMlAj38RmL/rMyASdhN99xA9h3DOPbY8lMngyVbodq06jb1pTqItIdtENDON9FWaq5k2cg13CqpSP33v//lySef5IUXXqBZs2b866vyT7yp3oMolW6/rgFe8a4QXQNcgmdo9Hkz+wrIBy7y9ratNLPX8SR1hcAIrRCtPiqbvITasmPklBVh35+RbhQUla9K4Lt2z4MacvHFF3Pfffdx2GEH4OnYrbpota8SXZEmbGuADEAJWw2XLJProzE0u2DBAgYMGECLFi3Izc2lefPmEfcg6ukzpVSq/XLOrQC6BHjp/CDnjwXGVjQ4SW5VTV6C9W4Fe/j1V792LerXqRWwncnNzaVv374sW7aMCy64gMMOO6yitxZUVdtXtY+xEWnCtgtYYWbv49foOeeujUlUkrSisZFrNH6Yw5V7Cefdd99l4MCBHHTQQbz//vs0b94ciKwHUU+fKUftl1RarObtBnr4LWt7XgEr7jq13PGtW7dyyimn8OWXXzJt2jQGDhxY6TgCCboKNTePruPnh2y71T7GTqQJ20zvL5EqTT6t6g+zr55dMKGGZn2J4nfL/sPmGWNpdWBbFixYgP9KvUCNqAE92v9xTnVceFHNqf2SSovVvN1I6o0Gas+2bNlC7969WbVqFTNmzOC0006rUhzBPjdQTMYfD8vB2m61j7ETUcLmnJvkncdxiPfQaudcZBu7iPipyg9zqFVV4GlMQtUK9b03rV5j6mQfhp02hlOfXE7uroJST4tL123llY9/LNmTwQFvLMuhS6u9GNQ5O6kWXkh4ar+kKmI5b9f38BuobQs21aRevXoccMABTJgwgVNOOSXgdas6ihHswbXsPjWB2m61j7ET6SrRk4FJwFo8f28tzOwibeshFRXJD3OwxibQvkT+HMF76SbMXc2On9eSsXcL6uzXluZD76PYrFwtPYAFqzaHbJiSZeGFREbtl1RFJPN2q5ogRTLVZOPGjdSvX5/GjRsza9asUnVB/UVjSDJQPJFOQwl2rgO6jp+v+WxVEOmQ6N+BU73L2/EuZ58MHB2rwKR6CpfsBGtslq7bWqq8SyDZIRKm/334Dlvefpi9+4+kQYceARs7X1IWLqlMloUXEjG1X1Jp4ZKpaM3ZCjXV5Mcff6Rnz54cfPDBzJkzJ2iy5oszGkOSZePpOn5+RA+qoebmaT5b1URa6SDD19gBOOf+h2fVlUiF+KoK+PNPdoI1NpM/WR/yuqESphdffJEtb/+dOi0Op17b40Nex9cgB+I7rqoGKUftl1TJoM7ZLB7dkx/G92fx6J4Rz9mKhqff/ohDOh3Hmg0/s/HA08NWlYnVkGS4ttvHv30MJJp/NjVNpD1sS81sIvCS9/s/46mtF5KZZQHPAYfj6RG91Dn3kZldA/wVz35Fs51zN3vPHwNchmdzy2udc3MrcC+SAsI9rQZrVIpcwCo/gKdEyl2ndwiYMD3zzDNceeWVHHl8N/K6j2RPmP+nffGE60HTrt8ppVLtl0gkYpUgzView5gX3uXbF2/GFexhn2Fj2dGoddgeqlhN2ajIDgG+9rHN6NkB67NpPlvlRJqwXQ2MAK7FMwfkA+CJCN73KDDHOTfEO+m3npn1AAYCRzjn9phZcwAzOwwYBnTAs1P4PDM7RJtPVj+hkp1gjU26WcCkLSvTk6xNmLuakVNWlGpEVq1axVVXXUX//v2ZNm0ac775lbtnrgxaCDkj3di5p5CRU1aQVS+DOrXS2J5XUOmtRyRpVLb9EgkrFgnSjOU5jH7jC3544wFcYT77nDuW2s0PBMIPb8ZyykZFH1Q13ze6zIXouajShc0aAZ8DBzq/D/HuBv6Mc25emfPHADjnxnm/nwvc7Zz7KNhndOnSxS1dujQW4UuCBFstddbR2byxLCfi474hyjlz5tCzZ09q164NBJ+HkWaepLCg2AW8TtkYtSlk4pjZMudcoA1tU47asNTk3wZk1cvg992FpdqOjDSjQd1a5VagR8rXThVs24grLKB2s1alXjfgh/H9I4pv/6xMerRvxoJVm+PeZgVrz2vyFJKqtF8he9jM7HXn3Dlm9iXlV/TinDsixNsPBDYDL5jZkXiGIK7Ds7T+JDMbC+wGbnLOfQpkAx/7vX+D91jZmIYDwwFatmwZKnxJQaG63bu02qvc8UDzR37+7+vc9tOXDHpuFH379i31WrCu+GIHxWUeXgI9yWpTyNRRxfZLJKCybcC2XQVkpBtZmRlszyugcWYGO/MLA65Aj6SNWLFiBV/NeJKs7heR0WT/gOf4iscHe3D07wlLZJsVjY3W5Q/hhkSv8/4+oJLXPgq4xjn3iZk9Coz2Hm8CHA8cA7xuZgfieWgoK1Aj+wzwDHieTisRl8RItHqeQhUqLnvcvx6fc47tiyezffGrFG49FRhV7hqRlIPxV/ZcbQqZUqrSfokEFKgNKChy1K9TixV3nUrX8fPLTbuItI1YunQpp556KruLMyg+ZiDp9ZuUOyczI50e7ZtFnIQlus3SfN/oCblK1Dn3k/fLvzjn1vn/Av4S5tobgA3OuU+830/Dk8BtAKY7jyVAMdDUe7yF3/sPADZW7HYkUXxPcTm5eTj+aEDCrWiqKt9cCOccuR/8i+2LX6X+4b3pOPSmgOcHW+kUbJF8epnl89oUMnVUsf0SCShcG1CResRdx8+nzejZdB0/n/EvzqRXr15kZWXxz1feokFW03LXaFIvg3GDO7Jg1eaIV6aqzao+It3WI9B2yv1CvcE59zOw3sx8Mx17AV8DM4CeULIfUm1gC57SMcPMrI6ZtQHaAksijE8SLNZL24MZ1acddWulsW3BRHZ8PJUGnfpywMCR3NwvcCHkQZ2zOevo7JJELN2Ms47ODriSCcqvTg235YckpQq3XyI+ZROrrHqBV5r72oBI2oiyD7jff/kpt105jPpZe7No0SKG9z++3NZBjwztxPI7T61wtRW1WdVHuDlsV+N5Ej3IzL7we6kh8GEE178GeMW7QnQNcAmwE3jezL4C8oGLvIsSVnoXJHyNZ7uPEVohmjoS9RQ3qHM2xcXFXD1nNxx9OoefdR03920ftAt+xvIc3liWU5KIFTnHG8tyaFIvI+DGvGX3EtKmuakjCu2X1HCB5n9lpBkZ6UZBUekFSr42IJI2ouwDrsvfTXrWvmT/+QFatPAMNFVmNX2gJExtVvURbg7bq8A7wDg88898fnPObQ13cefcCiDQaojzg5w/Fhgb7rpSObFc3ZiI5dvFxcVs2bKFwUe3YOAnb5OWlhZyB3AI3hNYp1YamRnpYRs1TaJNKVVqv0QCzlcrdmRlZlC/Tq2gE/597w1aZsrbVhb9vo30Bk3IPKgLddt0Zktx6ekaPoFWfQZaHR8oCVObVX1EtK2HmR0PrHTO/eb9viFwmN/8tITQkvjIxXp5daDrl13aHs2l5UVFRVx22WV88MEHLF++nMaNG5eKJVjjFGwjRwP+MbSTGrUUUNFl8cnafoHasGQXqr0Ita1GOF3Hz+fbpYvY8tYDNB14M/UOPg7w9OgvHt2z1LmhtjpKxFYdUjUx29bDz5N4Fgz47AxwTJJYrFcKlX2KC7S0/eWPfyw5v7JLy2csz+HBf3/NF6+OZefXixh21Y00atSo1OuhVk+F6gnUaqZqS+2XVEoktY8r85B3UsYaPnzzfmo3b0OdbM9822A9ZMHa7gWrNpdL7qR6i3TRgflvfuucKybyZE+SQDzmmPnX26tfp1apOR6BBFqUUHaCr/8q0xnLcxg9dTmfv3QPO79eRFb3i1jRtDdvrfhjMXG4xQ+R1sOTakXtl1RKj/bNyq0g97UXM5bnMGrq56VWxo+a+nnYlfFTpkzhoVuuot3hR9LpyoeoldkwZD3iYNsQ5eTmBW0rYyFU2yzxEWmjtcbMrsXzVAqeibxrYhOSxEK855hFmgj6nxeud2zC3NX8tGASu1YvpknPy2l0zKByvYThElPN56iR1H5JhfkWKPk/dhpw1tGenvhO97xbqroBeOa33T1zZdD25Msvv+S8887jxBNP5O2336Zhw4ZhYzACbEjqjcXXpsd6M1xtGJ4cIu1huwr4E5CDZ7+04/BWG5DUEEnPUjSfoCJNBP2XyIfrHduYm0ejY8+i6emjaHTMoJJz/JO0SJaw+/cELh7dUw1O9af2SyosUHvkgAWrNgMErUkc7DhAx44dmThxIv/+97/DJmu+GIKNU5Q9HsttlBK1bZOUFlHC5pz7xTk3zDnX3Dm3j3PuPOfcL7EOTqJnUOfscvv6+HfBR3vj20AJYiC/7y4s+YxgvWMbNm/jzjvvZN/66aRnNqT+Yd1Lve6fjAX6XMMztCE1k9ovqYxoTiN57rnnWLFiBQAXX3wx9evXr1IMoc6PxdClNt9NDuH2YbvZOfegmf2TwGWiro1ZZBJ1oSbVV3RRQrjJtoGGHrfu3ENeQXGp6xQUu5LPCDRsW5yfx/a37uO+H77gzv97mSn5e4Vcyj6oczZL123llY9/LPkH64A3luXQpdVe6k2rQdR+SVWEm0YSbO/GJmU21n300Ue5/vrrueiii3jxxRejEkO6WblNvQEaZ2bEZOgyEds2SXnheti+8f6+FE/x9rK/pJqoyBNUpL1xZYced5dJ1sp+RtneseI9u9gy9S5+X/slL7/8Mnf/5byQvYQ+C1ZtjutwgSQttV9SaeGmkdx1egcy0ksvSchIN+46vUPJ9xMmTOD666/nzDPP5JlnnolaDOce1yJwiT0jJkOXWqyVHEL2sDnnZnl/nxSfcCRRKvIEVdktQsJ9hn+v3PqfN7Nt+j3k//wtr732GmeffXbJOeGeFNV9L6D2SyrPN4KQV1BU0puVXcENcu+77z7uuOMOhg4dyksvvURGRuCSVqEE+gzffpaBYhs5ZUXA61S17dNireQQbkh0FoEXqADgnDsj6hFJQlSkfEllE6Ie7ZuVGqqEP1Y6dR0/v6QBGNQ5m6+//ppeU7bx0tSpDBo0qEL3ou57AbVfUjllV0QWOVfSFpZNUII9QBYWFrJ48WLOP/98XnjhBWrVqvwuMv6fES62CXNXR7Xti2V1HKm4cP+KHvL+PhjYF3jZ+/25wNoYxSQJUJEnqGAJUePMDLqOn1+yca4ZpaoclF0iD3/8b+obVs3btZNhf2rLYYcdxvfff0+9evUqfC/RqJ2nhqpaUPslFRZsBOGeWcG36/BxzrFr1y7q16/Pm2++SUZGBunp4RdfVTU23+hGNOuGaiuP5BNuSHQRgJnd65zr5vfSLDP7IKaRSdxFutN/oEYhI83YmV9YsqTdf2l7Tm5euZ61QH7ftoXLBl/Fmmuv4LbbbqtUsua7D6h8970aqupB7ZdURrCRgm27CpixPCdoG+Cc44YbbuCDDz5g0aJFNGjQIG6xxWKfyVhXx5GKi7SftpmZHeicWwNgZm0A7ZNQQwVqFHb5laEKJFyyVvjbFja9dhtFv23hT3/6U1RirGyjooaq2lH7JRELNoIABG0DiouLueaaa3jiiSe49tprI962I1qxld1nMhrtlOYCJ59IN84dCSw0s4VmthBYAFwfq6Ak9YRK1sIp3PELm14dQ9HvW+lw6QP06NEjipFVXLI1VCoJU2VqvyRioYYPA7UBxcXFXHnllTzxxBOMGjWKRx55BLOyBa2iF1u8VmtGsgm5xFdEPWzOuTlm1hZo7z20yjm3J3ZhSTILNGQYrHyKv0DnuMJ8Nk2+jaK8HbT88/3ce+WQGERcMZEuWojHPDcNz1ad2i+piEGds7l75sqAFQsCJSu33347zz33HLfffjt/+9vfYpas+WKD+KzWjOZ8OImOiBI2M6sH3AC0cs5dYWZtzaydc+7t2IYnyShYyZZQSVtmRjpnHZ3NglWbyy1P33Xieezf+mDuufT0cg1PIib/R9JQxSuR0vBs1an9koq6+4wOEScrV155Jfvssw/XXXddXGKL1pBnJJ8D2sojmUQ6h+0FPBtNnuD9fgMwFVCDVwMFGxp0eDazDbRKtOwP+qpVq1i7di33je4L9Ax4vUT1LkXSUMUrkUq24dkUpfZLKiRcG5Cfn8/EiRO58soradWqVdyStXiLV3IokYk0YTvIOTfUzM4FcM7lWSz7fSWpBRsyzM7KZPHowMmXv6+++opevXqRmZnJ6tWrqVOnTsDzEtm7FK6hilcipT3lokLtlwAV67EP1gbs2bOHs88+m1mzZtG2bVt69+4d67BFgMgTtnwzy8Q74mVmBwGaA5IA8RoiDPU5VZnbsGLFCnr37k2dOnWYO3du0GQNkrt3KV6JlOaRRIXaL4m4x97X9uXk5pWrJNCn/V4MHjyYOXPm8PjjjytZk7iKdJXoXcAcoIWZvQK8D9wcs6gkoEhreMb6cwZ1zo6opmdZS5cupWfPntSrV49FixbRrl3opCOZVyn1aN+Msl00sUikKvtnLaWo/ZKQPfY+/m0fUFJgPSc3j1umfMoJPfowd+5cnn32Wf7yl7/EL3gRIuhhM7M0oAme3cKPxzO3/Drn3JYYxyZlxGuIMJLPiXRug39PXf6H/6J2vYZ88MEiWrduHfa9yda75P/kXXaBhQFnHR2b+R6aR1J5ar/EJ5Ie+0Btn8+OjWv4bvmnvPjii1x44YUhP0uVUiQWwiZszrliM/urc+51YHYcYpIg4jVEGK3P8T2t7tqTj6WlU+uE86nX5UxWbMugbL4WqoELdDzeDWLZ4ZRAJbYWrNocs8+XylH7JT6RTGMI1Ma54iIsLZ06+7cje/hzXHjhn0N+jrbikViJdEj0PTO7ycxamNlevl8xjUzKidcQYbQ+Z8Lc1Wz97jM2ThxBQe7PmKVRULthqSEICD0EO6hzNotH9+SH8f1LFjR0/tu7XD9lRYWHhquyAW2oJ2+fZJhbJwGp/ZKINp0t28YV7f6dn1+5md+/eBeAFtn7hf2cSIZeRSoj0kUHl+LpRCg7aH9gdMORUOI1RBjpPmTheri+/+y/bH5zLLWy9iMt44/FBWUTm2AN3I2vfw788VRa9sm17Pllh4b9Y2ycmcHO/EIKiv6Yk1KRp95IkrFkmFsnAan9koi26/Fv+4rydvDLlDvI37KOtMzGEbe1ybxYSlJbpAnbYXgauxPxNHz/AZ6KVVASWLw2Mgz3OZF0+b/99tv8Mv0+au3dgn2G3kt6vcYl1y+b2ARryIqcK3XdcL1c/tcpG2OgXcsrMv8vVH1B0MrNJKf2q5qp7JSIcPNBfa+NfeNjPp98K4Vbc2g++A4OPupEerRvxoS5qxk5ZUXIz9RWPBIrkSZsk4AdwGPe78/1HjsnFkFJcPHc5TrY54RblPD+++8zePBgDmp3GNbvNvJr1Ss5L1BiEyoZ8r9uuCdU/wYxkiFMiPypN1Cvo2/hQbYmFSe7SrVfZpYFPAccjuev+lLn3Efe124CJgDNfAsYzGwMcBlQBFzrnJsb9TuRmM8RO7VdE26dfgf2288cevH97Gp6KDv3FDLl0/UR9dAn22IpqT4iTdjaOeeO9Pt+gZl9HouAJPmF6/Lv0qULl19+OePGjWPBmt/DPgkHauACXTdUYle2QYw0EYv0qVdlWlJaZduvR4E5zrkhZlYbqAdgZi2AU4AffSea2WHAMKADsD8wz8wOcc6Ff2qQCqnqavlwvXP16tXjmFMHs3NrI3Y2PRSoWA+92gqJlUgTtuVmdrxz7mMAMzsOWBy7sARivzS8stcPljhlbvqSXbt60LhxY5544gkABnVuHPaavtdvfP3zkn2Pyn4eBE/ssjIzuPuMDqU+J9wQJlT8qVfba6SsCrdfZtYI6AZcDOCcywfyvS//A88+bm/5vWUg8Jq3qPwPZvYdcCzwURTvQ6jaHLFQvXOd9ypky5YtHH300Xy3b0+sbvjrBftMtRUSC5GuEj0O+NDM1prZWjyNUHcz+9LMvohZdDVYrDfJrcr1A6222rNyHqsm3crYsWMrFc+gztn8/ZwjQ67iCrSJ7CNDO7HirlMD9tqVvVZGmtGkXoY2oK15KtN+HQhsBl4ws+Vm9pyZ1TezM4Ac51zZHrpsYL3f9xu8x8oxs+FmttTMlm7erK1gKqoqq9iD9c7dO3kB3bp14+yzz6agoCDqPfQi0RBpD1vfmEYh5dwza2VMN8mtyrBC2S7/tNXz+PntR+jTpw+33357pWOKZCgh0idXDUuIn8q0X7WAo4BrnHOfmNmjwN14et1ODXB+oNqk5buLAefcM8AzAF26dAl4jgRXlTligRKxgq05fD75VhrXdrz33ntkZGTEpIdepKoiSticc+tiHYj8YcbyHLbtKj9nAqK3NLyqS899idNjjz3GdQ88woABA5g6dSp169atUlzRHErQsIRApduvDcAG59wn3u+n4UnY2gCfe2vHHwB8ZmbHes9v4ff+A4CNlY1ZgqvKw1jZRKxgy3o2vXYrhmPBgkUcccQRQOCkMCPNaFC3Frm7CvQAKAkRaQ+bxFGoDRaj1QUfjaXn27ZtY+zYsQwePJjJkydTu3btqMSmsi6SaM65n81svZm1c86tBnoBnznnevnO8Q6vdnHObTGzmcCrZvYwnkUHbYEliYi9Jqjsw1jZRGzHp2+CwcMvTi9J1nzXB/XQS3JRwpaEQvVyRasLPtg2FT3aN4vo/c45mjRpwkcffUSLFi3IyMiISlyBJgVfP2UFd89cWW5hgUiMXQO84l0huga4JNiJzrmVZvY68DVQCIzQCtHk42s/Hpyzip+27+bws0dySecshvc/PuC5am8kmShhS0LBer+yMjOiOly4dN1WXvn4x5KJNg54Y1kOXVrtFfRznHPcdddd5OfnM27cOA48MLqbxQfbPy03r0D1+CSunHMrgC4hXm9d5vuxQOVW3Ujc7F+QQ8acv7F0+nT23nvvRIcjErFIV4lKHAWreXf3GR2i+jkLVm0uNys6VM075xxjxozh3nvvZfPmzbgAW3BUVajeRdXjE5Gq+PDDD+nduzcbNmxg586diQ5HpELUw5aE4jV/oiILD5xzDDx/OLNefY4GnU9jQbMzOfDWd6K+y3+41VmqxycilbFo0SL69+9PdranGssBBxyQ6JBEKkQJW5KKx/yJiiw8OOPPV/D25Ik0PPoMmvS6gmLvLgbRLgsTruqB9j0SkVACLVra67fv6devH61bt+b9999nv/32S3SYIhWmIdEaLNjQa6CFDf+zbBodN4Qmva7Au6VBiWgOVfo2x21Sr/wiBu17JCKhBNsQ/Lvd9enbty8LFy5UsiYpK6YJm5llmdk0M1tlZt+Y2Ql+r91kZs7MmvodG2Nm35nZajPrE8vYJHDlAP/d/4uKiliyxLMzQX6LY2ly8sXlkjWfaA5VDuqczfI7T+WRoZ2CxiYiUlbZRUt7Nq5m1558Jn62nenTp9O8efMERlc5M5bn0HX8fNqMnk3X8fOjVu1GUk+sh0RVPDkKYrkvWbCh18LCQi688EKmTp3KypUrw84ti8VQpZbVi0hF+LdRO1cvZsvMB2n8p2FY13MTGFXlhap9qrax5olZD5tf8eSJ4Cme7JzL9b7sK57sv8ywpHiyc+4HwFc8uUaLdU3RQPLz8xk2bBiTJ0/mvvvu45BDDgk4fOqjoUoRSQbp3hGAnV8vYstbD1Bnv0No1OWMkuOpJlQJQal5YtnD5l88+UhgGXAdnh3Dc5xzn5cZXssGPvb7PmDxZDMbDgwHaNmyZWwiTyLhfmCj3fO2Z88ezj77bGbNmsXDDz/MyJEjgdIrV3Ny80g3o8g5srMy6dG+GRPmrmbklBXaEVxEEqbIOX7/8n1+fedR6hxwGM3PupO0OvUoCrAFUSpUVKlqCUGpXmKZsMWkeHJNK5wc7AfT19MW7a7yl19+mVmzZvF///d/jBgxotRrgYYo1WUvIsmieUY+P857irotO9Js8B2k1fbUNs4uM2UjULs1csoKlq7byn2DOsY97mCiUUJQqo9YJmwqnhwFwX5g081i0vN26aWX0r59e7p27RrR+aF6AJWwiUg83Tr4WEb+9ABFjbNJy6gDBJ6yEajdcsArH/8YsNJLonrjAm1zpCkoNVfM5rA5534G1puZ71+Wr3hyc+dca29Zlw3AUd5zZwLDzKyOmbVBxZOB4FtvBOrihz96uCoy5+33339n6NCh/O9//8PMIk7WQF32IpJ4jzzyCI8//jiDOmfzjxGDadEsK+Tq8mDtk4Ny88MSMY/YJ9xKfqlZYr1KVMWTqyhY1QPfXLKyQvW8Bfoh37FjB/369eOTTz7hnHPO4ZBDDqlQfMF6ABtnZtB1/Pyknh8iIqnvwQcf5JZbbmHIkCH85S9/iWh1eahV72WTuUSPImi1vPjENGFT8eToCPYDG6irPFiFgEBPlNu2baNv37589tlnTJkyhbPOOitsLGWHBnq0b8Yby3JKfW5GmrEzv5DcvAJA89pEJDbuvfde7rzzToYNG8ZLL70UdJ/Iskb1acfIKSvKT5Km/PwwjSJIslClgxQVrKu87ORan7KN0NatW+nduzefLV/BwefexU2f1g27KWOgoYE3luVw1tHZpeJoULcWBUWlm0ItRReRaLrjjju48847ufDCC3n55ZepVSvy/odBnbP58/Ety610CzQ/LNgEf038l3hTLdEUVpGet7KNUO3atSmq3YB9z7qdvP07A+F7woINDUz+ZD3FzpUMfY6csiJgvHoiFZFo2Wuvvbj88st5+umnSUsL3PcQarHAfYM60qXVXmEXE2jivyQLJWzVTLA5b77jP//8M/Xr16dhw4bUO+NOcrfvLvX+UHMzgiVcvgUQvoQvq14G23YVlDtPT6QiUhXOOb7//nsOPvhgRo4ciXMu6DBoJFsORTI/LFybKhIvStiqmVBPlBs2bKBnz560bduW2bNn81OZZM0nWGIWrjwVeBK+OrXSys2n0xOpiFRFcXExI0aM4OWXX+aLL76gTZs2IeesRXOxgCb+SzLQHLZqJNTy83Xr1tG9e3d+/vlnbrvtNqDiczNClafytz2vQEvRRSRqioqKuOKKK3jqqacYMWIErVu3DvseLRaQ6kY9bNVIsCfKeycvYMuU29ixYwfz5s3j2GM9JVorOjej7NBAmrc8VVn7Z2XqiVREoqKwsJBLL72Ul156iTvvvJO77747otWgqhIg1Y162KqRQE+Ozjm+evk+fv/9d95///2SZA0qtynjoM7ZLB7dkx/G9+fv5xwZcFNfDX2KSLQ89dRTvPTSS9x3333cc889Fdq6Q+2TVCfmguyYnwq6dOnili5dmugwoq6yZVC6jp8f8Ilyr4ItvHhRZ4444oiYxto4MwMzyN1VoIm5EjNmtsw5F3R/x1RSXduwaCooKGD27NkMGjSowu9NhQLvUrNUpf1SwpZkyq5sAs9TYSRzwPzfm795LbtW/Zd9e1zI+LOOiHkjVZW4RSpCCVv1t2fPHm6++WZuvfVW9tlnn0SHIxI1VWm/NCSaZEKtbArHN8TZeNcGNk2+lbyv5jG6x/5xSZiqEreIiE9eXh4DBw7kscceY/78+YkORyRpaNFBkqnqyqbswo38+NJoDmiWxfz58znooIOiGV5QWpElUvNEe8hx586dnHHGGSxYsICJEydy7rnnRjFakdSmhC0J+Dd6oVZehvPhhx/Sr18/mjZtyvz582nVqlUswg1IK7JEapZINqatiN9++40BAwbw3//+l0mTJtHw8J50HT9f889EvJSwxVGgp1EoXUoqULIW6cqmLVu20KJFC+bMmcMBBxwQ3eDDUPkWkZolmhvTAuzevZvc3FxeffVV6hxyYlSTQZHqQAlbnAR7Gq2bkVau0QNINytVnzNUI7VlyxaaNm3KGWecwWmnnVahIsjRovItIjVLRadBBBs+3b59O/Xq1aNZs2YsW7aMWrVq0XX8/KgmgyLVgRK2OAn2NBooWQModo4fxvcPe905c+Zw9tlnM2XKlIQlaz7aLFek5qjINIhgD6w7tm3lkZsuokOHDrz00ksl7ZfmxIqUp4QtTira0ASb++X/lFr3p+V8N/leOh7egeOOOy7sNbUnkYhES0WmQQR6YP0991eu/vNfKM79iXHjxpV6TXNiRcrTth5xEqyhycrMiHg3bv9aob+vXsyql+8mrWkrCvrcxn9+DFzIPdB7y9YZFRGpqIpUSin7wFr4+1Y2TR5D3q8bGfPIJP62PIM2o2fTdfx8ZizPUZUCkQCUsMVJsAbo7jM6RNzo+Z5S8zevZctbD1Bnv0PYZ+h9/Obqhk2+tE+aiESbf6m6xaN7Bu2x939gdc6x+Y17KfxtCwf+eSyTNzQq9yAJVLhsnkh1pyHROAk3KT+Shsj3lJrRtBV7nfoX6h/ajbQ69YDwE3I1J0REEsV/+NTMaNJrOHVqpdGwzeFs21VQ6lxfWxYqAdT0DqmJlLDFUVUn5dda8x92NmxB7Watadipb7nXQyVfmhMiIokyqHM2P61fy0MvTKO4XW8OOvwoRvVpx8gpKwKeH6oti/b+byKpQkOiKeKJJ57gu6kPsPPTN4OeEyr50pwQEUmU1atXc9/V57D9Py/x6U3HlfSeBWuzQrVlmt4hNZUSthTwyCOPMGLECE4//XSee/YZsjIzyp0TLvmqyARhEZFoWblyJd27d6ewsJAFCxbQtGnTktcq8yCp6R1SU9WIIdFUnu/wwAMPMHr0aM466yxeffVVateuzdnHHVipe9I+aSIST59//jm9e/cmIyOD999/n0MPPbTU65XZcFvTO6SmqvYJWyrPdygsLGTevHkMGzas1KaSoORLRJLf0qVLyczM5P3336dt27YBz6loW6YyeFJTVfuELdr17uLBOcfu3bvJzMxk5syZ1K5dm/T09PBvFBFJAnl5eWRmZnLZZZcxdOhQGjRoELVrqwye1FTVPmFLtfkOzjlGjx7NokWLeP/996lfv36iQxIRidjixYsZMmQI06ZNo2vXrlFN1nw0wiA1UbVfdFCZVUiJ4pxj5MiRPPjggxx11FFkZiZfjCIiwSxcuJA+ffrQqFEjWrdunehwRKqVap+wpcp2FsXFxYwYMYJHH32U66+/nscff5y0tGr/1yOStMwsy8ymmdkqM/vGzE4wswne778wszfNLMvv/DFm9p2ZrTazPgkMPSHmzZvHaaedRqtWrVi0aBHZ2eoBE4mmap8RpMp2FnfccQdPPvkkN998Mw8//DBmluiQRGq6R4E5zrn2wJHAN8B7wOHOuSOA/wFjAMzsMGAY0AHoCzxhZjVm4uny5csZMGAAbdu2ZeHChey7776JDkmk2qn2c9ggNeY7XHrppTRp0oQbb7xRyZpIgplZI6AbcDGAcy4fyAfe9TvtY2CI9+uBwGvOuT3AD2b2HXAs8FFVY0mFbYmOOOIIbrnlFq699lr23nvvRIcjUi1V+x62ZFZQUMDEiRMpLi7moIMO4qabblKyJpIcDgQ2Ay+Y2XIze87Myq4AuhR4x/t1NrDe77UN3mPlmNlwM1tqZks3b94cMgjftkRli6PPWJ5TmXuKulmzZrFhwwbS09O55557lKyJxJAStgTJz89n2LBhXH755cyfPz/R4YhIabWAo4AnnXOdgZ3AaN+LZnYbUAi84jsU4Bou0IWdc88457o457o0a9YsZBDJXIbp1VdfZdCgQdx2222JDkWkRlDClgB79uxhyJAhTJ8+nUceeYTevXsnOiQRKW0DsME594n3+2l4EjjM7CJgAPBn55zzO7+F3/sPADZWNYhk3ZZo0qRJnH/++XTr1o3HH388obGI1BRK2OIsLy+PgQMHMmvWLJ588kmuu+66RIckImU4534G1puZbzl5L+BrM+sL3AKc4Zzb5feWmcAwM6tjZm2AtsCSqsaRjNsSPfvss1xyySX07t2b2bNnx2SfNREpTwlbnK1YsYJFixYxceJErrrqqkSHIyLBXQO8YmZfAJ2A+4H/AxoC75nZCjN7CsA5txJ4HfgamAOMcM4VBbxqBSTbtkT5+fk88cQT9OvXj5kzZ1KvXr2ExCFSE9WIVaLJoKioiPT0dE444QTWrFnDfvvtl+iQRCQE59wKoEuZwweHOH8sMDaaMSRTGaaioiJq167NvHnzaNCgAXXq1Il7DCI1mRK2ONi+fTv9+/fnyiuv5IILLlCyJiIRS4ZticaNG8dHH33EtGnTtBJUJEE0JBpj27Zto3fv3ixZskRzPUQkpTjnuOeee7j11ltp2LChqq+IJJB62GJoy5YtnHLKKXz99ddMnz6dAQMGJDokEZGIOOe44447GDt2LBdffDHPPfcc6ek1pniDSNKJ6eNSTa7Ft2vXLnr27MmqVauYOXOmkjURSSl/+9vfGDt2LFdccQUTJ05UsiaSYLHuYfPV4htiZrWBenhq8Y1xzhWa2QN4avHdUqYW3/7APDM7JBorrRKhXr16DB06lOOPP55evXolOhwRkQoZMGAAu3btYty4cRoKFUkCMUvYkqkWXzxt2LCBX3/9lSOPPFI7gItISikuLmbOnDmcdtppHH300Rx99NGJDklEvGL52BSTWnwVqcMXb2vXrqVbt24MHjyYgoKCRIcjIhKxoqIirrjiCvr378+iRYsSHY6IlBHLhC0mtfgqUocvnr7//nu6d+/Otm3beO2118jIyEh0SCIiESksLOTiiy/m+eef56677qJbt26JDklEyojlHLZAtfhGQ6lafL1iXYsvHlavXk3Pnj3Zs2cP8+fPp3PnzokOSUQkIgUFBZx//vm8/vrrjB07lltvvTXRIYlIADHrYUuWWnzxMH78eAoLC1mwYIGSNRFJKYsWLWLq1Kk89NBDStZEklisV4n6avHVBtYAlwCfAnXw1OID+Ng5d5VzbqWZ+WrxFRKlWnzx8OSTT5KTk8NBBx2U6FBERCqkd+/efPnll3To0CHRoYhICDFdq+2cW+Gdb3aEc26Qc26bc+5g51wL51wn76+r/M4f65w7yDnXzjn3TqhrJ9pnn33GKaecwrZt26hbt66SNRFJGXl5eQwaNIj33nsPQMmaSArQ5jqVsGTJEnr16sW3337L9u3bEx2OiEjEdu7cSf/+/Zk5cyY5OTmJDkdEIqSErYIWL15M79692WuvvVi0aBGtW7dOdEgiIhH57bff6NevH4sWLeJf//oXF198caJDEpEIqZZoBfz3v/+lb9++ZGdnM3/+fLKzy20TJyKSlH7//XdOPfVUPv30UyZPnsw555yT6JBEpALUw1YBrVq14uSTT2bRokVK1kQkpdSrV49OnToxdepUJWsiKUg9bBFYvnw5RxxxBC1atODtt99OdDgiIhHbvHkzu3btolWrVjz55JOJDkdEKkk9bGG89dZbHHfccYwfPz7RoYiIVMjPP/9Mjx49OO200ygqSoldkkQkCPWwhTB16lTOO+88jjrqKEaMGJHocEREIpaTk0OvXr1Yv349s2bNIj09PdEhiUgVqIctiFdffZVhw4Zx3HHH8d5775GVlZXokEREIvLjjz/SvXt3cnJymDNnDj179kx0SCJSRephC2Dz5s1ceeWVdOvWjVmzZtGgQYNEhyQiErEbb7yRLVu28N5773H88ccnOhwRiQIlbAE0a9aMefPm0bFjR+rVq5focEREKuSZZ55h3bp1dOrUKdGhiEiUaEjUz+OPP86zzz4LwHHHHadkTURSxqpVq7jooovYvXs3TZo0UbImUs0oYfP6xz/+wV//+lfeeecdnHOJDkdEJGJfffUVJ598MnPmzGHDhg2JDkdEYkAJGzB+/HhuuOEGhgwZwpQpUzCzRIckIhKRzz//nB49epCWlsaiRYs4+OCDEx2SiMRAjU/Y/va3vzFmzBjOO+88Jk+eTEZGRqJDEhGJyLJly+jRoweZmZl88MEHtG/fPtEhiUiM1PiErXbt2lx00UX861//olYtrcEQkdSRkZHBgQceqJ41kRqgRmYozjnWrVtH69atGT16NM45DYOKSMpYu3YtrVq14ogjjuDTTz9V+yVSA9S4HjbnHNdddx2dOnVi7dq1AGrsRCRlLFiwgA4dOvD4448Dar9EaooalbAVFxdz9dVX889//pPLLruMVq1aJTokEZGIvfvuu5x22mm0adOGs88+O9HhiEgc1ZiEraioiMsvv5ynn36a0aNH89BDD+nJVERSxuzZszn99NNp164dCxYsYJ999kl0SCISRzUmYXvyySd54YUXuOuuu7j//vuVrIlIyvj5558ZMmQIHTt2ZP78+TRr1izRIYlInNWYRQfDhw+nefPmnHPOOYkORUSkQvbdd1+mTp3KiSeeSFZWVqLDEZEEqDEJW+3atZWsiUjKGjBgQKJDEJEEqjFDoiIiIiKpSgmbiIiISJJTwiYiIiKS5JSwiYiIiCQ5JWwiIiIiSU4Jm4iIiEiSU8ImIhKAmWWZ2TQzW2Vm35jZCWa2l5m9Z2bfen9v4nf+GDP7zsxWm1mfRMYuItWPEjYRkcAeBeY459oDRwLfAKOB951zbYH3vd9jZocBw4AOQF/gCTNLT0jUIlItKWETESnDzBoB3YCJAM65fOdcLjAQmOQ9bRIwyPv1QOA159we59wPwHfAsfGMWUSqt5SudLBs2bItZrYu0XEE0RTYkuggqiDV4wfdQ7KI9j20iuK1gjkQ2Ay8YGZHAsuA64B9nHM/ATjnfjKz5t7zs4GP/d6/wXusHDMbDgz3frvHzL6KQfzxVh3+nfpUl3upLvcB1edemlKF9iulEzbnXNJWQDazpc65LomOo7JSPX7QPSSLFL2HWsBRwDXOuU/M7FG8w59BWIBjLtCJzrlngGcgZf9syqku9wHV516qy31A9bkX7320ruz7NSQqIlLeBmCDc+4T7/fT8CRwm8xsPwDv77/4nd/C7/0HABvjFKuI1ABK2EREynDO/QysN7N23kO9gK+BmcBF3mMXAW95v54JDDOzOmbWBmgLLIljyCJSzaX0kGiSeybRAVRRqscPuodkkar3cA3wipnVBtYAl+B5yH3dzC4DfgTOBnDOrTSz1/EkdYXACOdcUQSfkap/NmVVl/uA6nMv1eU+oPrcS5Xuw5wLOM1CRERERJKEhkRFREREkpwSNhEREZEkp4StkoKUrZng/f4LM3vTzLL8zk+6sjWB7sHvtZvMzJlZU79jSXUPweI3s2u8Ma40swf9zk+q+CHov6NOZvaxma0ws6Vmdqzf+Ul1D2bWzhun79cOM7teJZw8qlN5q+rQ5kHqt3s+1aH980n1dhDi1BY65/SrEr/w7HJ+uffr2kAWcCpQy3vsAeAB79eHAZ8DdYA2wPdAejLeg/frFsBcYB3QNFnvIcjfQQ9gHlDHe7x5ssYf4h7eBfp5j50GLEzme/C7l3TgZzwbQz4IjPYeH53sPwtx/vtNyT+b6tDmBbsP79cp0e6F+ftIqfYvzL2kZDvojTEmbaF62CrBgpStcc6965wr9J72MZ69mCAJy9YEuwfvy/8Abqb0xp9JdQ8h4r8aGO+c2+M97tsnK6nih5D34IBG3tMa88d+Xkl3D2X0Ar53zq1DJZyqVXmr6tDmQeq3ez7Vof3zqYbtIMSoLVTCVjn+ZWuWm9lzZla/zDmXAu94v84G1vu9FrRsTRwFvAczOwPIcc59Xub8ZLuHYH8HhwAnmdknZrbIzI7xnp9s8UPwe7gemGBm64GHgDHe85PxHvwNAyZ7vy5VwgnwL+GUzPcQTcH+flPxz6Y6tHmQ+u2eT3Vo/3yqWzsIMWoLlbBVjq9szZPOuc7ATvzK1pjZbXj2YnrFdyjANRK9n0qge7gbuA24M8D5yXYPwf4OagFNgOOBUXj2zDKSL34Ifg9XAyOdcy2AkXifPEnOewDAPHuVnQFMDXdqgGNJcQ8xELKdCCCZ/2yqQ5sHqd/u+VSH9s+n2rSDENu2UAlb5QQrW4OZXQQMAP7svAPVJGfZmmD30Ab43MzW4onzMzPbl+S7h2DxbwCmO48lQDGegrvJFj8Ev4eLgOneY1P5o5s8Ge/Bpx/wmXNuk/d7lXCqXuWtqkObB6nf7vlUh/bPpzq1gxDDtlAJWyW4IGVrzKwvcAtwhnNul99bkq5sTZB7+Mw519w519p5CtRuAI7ynptU9xDs7wCYAfQEMLND8Exg3UKSxQ8h72Ej0N17rCfwrffrpLsHP+fyxxAAqIRTqL/flPuzqQ5tHqR+u+dTHdo/n2rWDkIs28JQKxL0K+QqkE7AUuALPD8kTfBMGlwPrPD+esrv/NvwrAJZjXflS6J/BbqHMq+vxbtaKhnvIcjfQW3gZeAr4DOgZ7LGH+IeTgSW4VlB9AlwdJLfQz3gV6Cx37G9gffxNLLvA3sl8z3E+e83Jf9sqkObF+w+yrye1O1emL+PlGr/wtxLSrWD3rhi2haqNJWIiIhIktOQqIiIiEiSU8ImIiIikuSUsImIiIgkOSVsIiIiIklOCZuIiIhIklPCJknJzNaaWdNkuY6ISEWoDZNoU8ImMWUe+ncmIilJbZgkC/0jlKgzs9Zm9o2ZPYFn88YWZjbKzD41sy/M7B6/c2eY2TIzW2lmw8Nc92oze9Dv+4vN7J+RXMcb01d+399kZnd7vz7IzOZ43/8fM2tf5T8EEUlZasMkGSlhk1hpB/zLeYr5tsNTduNYPDtaH21m3bznXeqcOxroAlxrZnuHuOY0YLDf90OBKZW4TlnPANd4338T8EQF3isi1ZPaMEkqtRIdgFRb65xzH3u/PtX7a7n3+wZ4Gr8P8DRMZ3qPt/Ae/zXQBZ1zm81sjZkdj6fMRztgsffliK/jz8waAH8CppqZ73CdiO5QRKoztWGSVJSwSazs9PvagHHOuaf9TzCzk4HewAnOuV1mthCoG+a6U4BzgFXAm845F+F1Cindo+x7PQ3Idc51iuSmRKTGUBsmSUVDohIPc4FLvU+CmFm2mTUHGgPbvA1Ue+D4CK41HRgEnMsfQwmRXGcT0NzM9jazOsAAAOfcDuAHMzvbG5uZ2ZGVvVERqZbUhknCKWGTmHPOvQu8CnxkZl/imcfREJgD1DKzL4B7gY+DX6XkWtuAr4FWzrkl3sNhr+OcKwD+BnwCvI3n6dbnz8BlZvY5sBIYWJn7FJHqSW2YJANzziU6BhEREREJQT1sIiIiIklOCZuIiIhIklPCJiIiIpLklLCJiIiIJDklbCIiIiJJTgmbxI2ZrTWz3lW8xsVm9t9oxSQiIpIKlLBJjefdaFI/CyIikrT0n5TEhZm9BLQEZpnZ72Z2s5kdb2YfmlmumX3uLc/iO/9ib82938zsBzP7s5kdCjwFnOC9Rq733NPM7GvvuTlmdpPfdQaa2Qoz22Fm35tZX+/xhWY21swWA7uAA+P2hyEiIlJB2jhX4sbM1gKXO+fmmVk28AVwAZ5dvnsBrwHt8SRQPwHHOOdWm9l+wF7OuZVmdrH3Gif6Xfcn4Bzn3H/MrAnQxjn3mZkdC7wHDAHeB/YDGjrnVnlr9R0I9ANW4/lZKIj9n4KIiEjFqYdNEuV84N/OuX8754qdc+8BS4HTvK8XA4ebWaZz7ifn3MoQ1yoADjOzRs65bc65z7zHLwOed8695/2MHOecfzmXF51zK51zhUrWREQkmSlhk0RpBZztHQ7N9Q5vngjs55zbCQwFrgJ+MrPZ3oLIwZyFJ9FbZ2aLzOwE7/EWwPch3re+ynchIiISB0rYJJ78x9/XAy8557L8ftV3zo0HcM7Ndc6dgmcYcxXwbIBr4D33U+fcQKA5MAN43e8zDoowHhERkaSlhE3iaRN/TO5/GTjdzPqYWbqZ1TWzk83sADPbx8zOMLP6wB7gd6DI7xoHmFltADOr7V2Q0Ng7rLnD79yJwCVm1svM0swsO0xPnYiISFJSwibxNA643Tv8ORQYCNwKbMbTGzYKz7/JNOBGYCOwFegO/MV7jfnASuBnM9viPXYBsNbMduAZRj0fwDm3BLgE+AewHViEZyhWREQkpWiVqIiIiEiSUw+biIiISJJTwiYiIiKS5JSwiYiIiCQ5JWwiIiIiSa5WogOoiqZNm7rWrVsnOgwRiaNly5Ztcc41S3QcIiLxlNIJW+vWrVm6dGmiwxCRODKzdYmOQUQk3jQkKiIiIpLklLCJiIiIJLmYJWxm1s7MVvj92mFm15vZXmb2npl96/29id97xpjZd2a22sz6xCo2ERERkVQSs4TNObfaOdfJOdcJOBrYBbwJjAbed861Bd73fo+ZHQYMAzoAfYEnzCw9VvGJiIiIpIp4DYn2Ar53zq3DUz9ykvf4JGCQ9+uBwGvOuT3OuR+A74Bj4xSfiIiISNKKV8I2DJjs/Xof59xPAN7fm3uPZ+MpAO6zwXusFDMbbmZLzWzp5s2bYxiyiIiISHKIecJmZrWBM4Cp4U4NcKxcZXrn3DPOuS7OuS7NmmkrJhEREan+4tHD1g/4zDm3yfv9JjPbD8D7+y/e4xuAFn7vOwDYGIf4RERERJJaPBK2c/ljOBRgJnCR9+uLgLf8jg8zszpm1gZoCyyJQ3wikiD//e9/Wbt2baLDEBFJejFN2MysHnAKMN3v8HjgFDP71vvaeADn3ErgdeBrYA4wwjlXFMv4RCRx5s2bx6mnnsqIESMSHYqISNKLaWkq59wuYO8yx37Fs2o00PljgbGxjElEEu+dd97hzDPP5JBDDuGFF15IdDgiIklPlQ5EJO4WL17MYYcdxoIFC2jevHn4N4iI1HDmXLmFmCmjS5cuTsXfRVLHrl27qFevHs45du3aRf369St8DTNb5pzrEoPwRESSlnrYRCQuJk+eTNu2bfn2228xs0olayIiNZUSNhGJuUmTJnH++edz8MEHs++++yY6HBGRlKOETURi6rnnnuOSSy6hR48e/Pvf/6Zhw4aJDklEJOUoYRORmHnrrbe44oor6NOnD7NmzdIwqIhIJSlhE5GYOeWUU7jrrruYMWMGmZmZiQ5HRCRlKWETkaj717/+xY4dO6hXrx533303derUSXRIIiIpTQmbiETVvffey0UXXcRjjz2W6FBERKqNmFY6EJGawznHnXfeyX333ceFF17ImDFjEh2SiEi1oYRNRKrMOcfo0aN58MEHufzyy3n66adJS1MHvohItKhFFZEq+/XXX3n11Ve5+uqrlayJiMSAethEpNKKi4sBaNq0KUuXLqV58+aYWYKjEhGpfvQYLCKVUlRUxPDhwxkxYgTOOfbZZx8layIiMaKETUQqrLCwkEsuuYSJEyfSrFmzRIcjIlLtaUhURCqkoKCACy+8kNdee417772X22+/PdEhiYhUe0rYRKRCLr74Yl577TUefPBBRo0alehwRERqBCVsIlIh5513HscddxzXXnttokMREakxlLCJSFh5eXl88MEH9OnTh/79+yc6HBGRGkeLDkQkpJ07dzJgwAAGDBjADz/8kOhwRERqJPWwiUhQv/32G/3792fx4sW8+OKLtGnTJtEhiYjUSErYRCSg7du3069fP5YsWcKrr77K0KFDEx2SiEiNpYRNRAJ6/fXXWbp0Ka+//jqDBw9OdDgiIjWaEjYRKcU5h5lx+eWXc9JJJ9G+fftEhyQiUuNp0YGIlNi0aRM9evTgiy++wMyUrImIJAn1sIkIAD/99BM9e/Zk3bp1bNmyJdHhiIiIHyVsIsKGDRvo2bMnP/30E3PmzKFbt26Vus6M5TlMmLuajbl57J+Vyag+7RjUOTvK0YqI1DxK2ERSTLSTopycHLp168avv/7Ku+++ywknnFDpuMZM/5K8giLPdXPzGDP9SwAlbSIiVaQ5bCIpxJcU5eTm4fgjKZqxPKfS19x777057rjjmDdvXqWTNYAJc1eXJGs+eQVFTJi7utLXFBERDyVsIikkmknRt99+y9atW6lbty6TJ0/mmGOOqVJsG3PzKnRcREQip4RNJIVEKyn6+uuvOemkk7jggguiERYA+2dlVui4iIhETgmbSAqJRlL0xRdfcPLJJ5OWlsZDDz0UrdAY1acdmRnppY5lZqQzqk+7qH2GiEhNpYRNJIVUNSn67LPP6NGjB3Xq1GHRokUceuihUYttUOdsxg3uSHZWJgZkZ2UybnDHoAsOZizPoev4+bQZPZuu4+dXaR6eiEh1Z865RMdQaV26dHFLly5NdBgicVXZVaLOObp06cKvv/7KggULElrIveyKUvAknqESPB8zW+ac6xLrGEVEkokSNpEaZO3ataSlpdGyZcuExtF1/HxyAsy7y87KZPHoniHfq4RNRGoiDYmKpLBIhhUXLlzIX//6V4qLi2ndunXCkzXQilIRkYrSxrkiKSqSjWrfe+89Bg4cSJs2bdi+fTtNmjRJWLz+9s/KDNjDphWlIiKBqYdNJEWF25Pt3//+N6effjpt27Zl4cKFSZOsgVaUiohUlHrYRFJUqGHFmTNnMmTIEDp27Mi7777L3nvvHefoQvP1AKruqIhIZJSwiaSoUMOK9es7unbtyptvvklWVlb8g4vAoM7ZStBERCIU0yFRM8sys2lmtsrMvjGzE8ysk5l9bGYrzGypmR3rd/4YM/vOzFabWZ9YxiaS6nq0b1buWEHuz/Ro34xevXoxf/78pE3WRESkYmLdw/YoMMc5N8TMagP1gNeBe5xz75jZacCDwMlmdhgwDOgA7A/MM7NDnHNFwS4uUpMtWLW51Pe/f/k+v77zKG/k/Y0urfbScKOISDUSsx42M2sEdAMmAjjn8p1zuYADGnlPawxs9H49EHjNObfHOfcD8B1wLCISkP8ctt9WzOHXfz9C3ZZHsLNpe8ZM/5Kc3Dwcf6weVSUBEZHUFcsh0QOBzcALZrbczJ4zs/rA9cAEM1sPPASM8Z6fDaz3e/8G77FSzGy4dyh16ebNm8u+LFJj+LbA+O2zt9k69/+oe+BRNB9yJxm1M0OuHhURkdQTy4StFnAU8KRzrjOwExgNXA2MdM61AEbi7YEDLMA1ypVhcM4945zr4pzr0qxZ+Tk8IjXFqD7tsC1r2PreU2S2PZ7mZ95OvcxMioJUL9GmtCIiqSuWc9g2ABucc594v5+GJ2E7EbjOe2wq8Jzf+S383n8AfwyXikgZgzpnw4jBjN79K7v360z23g0Z1acdE+au1qa0IiLVTMwSNufcz2a23szaOedWA72Ar/EMlXYHFgI9gW+9b5kJvGpmD+NZdNAWWBKr+ERSlXOOBx98kJ49ezLomGMY9PLd5c4JVFhdm9KKiKSuWK8SvQZ4xbtCdA1wCfAW8KiZ1QJ2A8MBnHMrzex1PEldITBCK0RFSnPOcfvtt3P//fdzzTXXcMwxx5Q7R5vSiohUP+aCzHdJBV26dHFLly5NdBgiceGc4+abb+ahhx7iiiuu4KmnniItLTHV5WYsz0lYQmhmy5xzXeLyYSIiSUKVDkRSgHOO66+/nscee4wRI0bw2GOPlSRr/slT48wMzCB3V0HMEqlIis6LiEh0qfi7SAooLCzkhx9+4IYbbuCf//xnqWTNf8+13LwCtu0qiOn+a+GKzouISPSph00kzioynFhUVMSOHTto0qQJb7zxBrVq1cLsjx1wAiVP/nyJVDR7vkIVnRcRkdhQD5tIHM1YnsOoaZ+XqkIwatrnAXvBCgsLueiii+jevTt5eXlkZGSUStYgsiQp2olUsO1BtG2IiEjsKGETiaN7Zq2koKj0Qp+CIsc9s1aWPlZQwHnnnccrr7zCsGHDyMysfJIU7URqVJ92ZGaklzqmbUNERGJLCZtIHG3bVRD2eH5+Pueccw5Tp07loYce4tZbbw16vUDJkz9fIjVjeQ5dx8+nzejZdB0/v0rz2gZ1zmbc4I5kZ2ViQHZWJuMGd9SCAxGRGNIcNpEkc+ONNzJjxgwee+wxrrnmmpDnlt1zLdAqUSDqqzoHdc5WgiYiEkdK2ETiKCszg9y88r1sWZkZJV+PHj2a4447jvPPPz+ia4ZLnrqOnx90VaeSLhGR1KAhUak2ojnsFyt3n9GBjLTSCwcy0ozRp7Rh/PjxFBUVkZ2dHXGyFgmt6hQRSX3qYZNqIVU2cw1UNuqvJ2bzxOjLWLx4MSeddBJdu3aN6mfun5WpYvAiIilOCZtUC6E2c02mhA1KD2Hm5ubSr18/Pv30UyZPnhz1ZA08CxNUDF5EJLUpYZNqIZbDfrGqm7l161b69OnD559/ztSpUznzzDOrfM1AVAxeRCT1KWGTaiFWw36xHGr97rvvWLNmDdOnT2fAgAFVulY4WtUpIpLalLBJtRCrYb9Ihlor2gO3e/du6taty7HHHssPP/xAo0aNqhRjNMWqN1FERKpGq0SlWojVZq7hhlrLFl8PV3B948aNdO7cmWeeeQYg6ZK1ityLiIjEj3rYpNqIxbBfuKHWiix2uPa593jy5osp2rmNuz/I5cfmX3LfoI5RjbcqUmnhhohITaMeNpEQwtXNjHSxw1+ffocnbrqAol3b2eece6l9QAde/vhHbp/xZWwCrwTt1yYikryUsImEEG6oNdiiBv/jv/32G0/dfBFuz072GXofdbLbl7w2+ZP1MY2/IiK5FxERSQwNiYqEEWqoNZLFDg0bNqTxcUOok30otfc5sNT7i5yLTdCVoP3aRESSlxI2ESq/OjLYHmcAnUc+z8+bt3Lg4UfR+Oj+FAfIzdLNyh9MEO3XJiKSvJSwSY1X1b3WyvbAzView8gn3+LHl8eQltmI2vs/TnpaesD3nntciyjcQfRovzYRkeSkOWxS44VaHVkZtzz1Jj++dAtWqzbNz7oDS0unGMjMSCvpUUs34/zjWybVKlEREUle6mGTGi+aqyMfmDSTb1+8hbTMhuwzbCwZWfuWvLa7oJgfxvevdJwiIlJzqYdNarxgqyDTzGgzejZdx8+PePPYh//vSdLrNWbf88aVStZCfY6IiEg4Stikxgu01xp4VnBGuuN/cXExAJknX8W+50+gVqPmAT9HRESkMpSwSVKZsTyHruPnV7hnqyrK7rUWaOVmqDlt7777LscccwybNm0ie++GpNfPKndOVmaGJvOLiEilKWGTpJHIWpaDOmezeHRPfhjfn+Ige6MFmtM2e/ZsTj/9dIqKikhLSwtaGeHuMzrEJG4REakZlLBJ0oj2as3KCjWnzT95nDFjBmeeeSYdO3Zk/vz5LN6QX3IPvl66QEXoE9GLKCIiqU2rRCVpJEsty0A7/oNnTptvf7aMjSs4++yz6dKlC++88w4Lf9hZ6j1FzpVUCSibrFVlzzcREamZ1MMmSSNZaln65rSFmsvWpUsXzj//fObOnUtWVlbEvYPJ0osoIiKpRQmbJI1g87+qsrqyssOPgzpnB5zLlrd2BRt+/Y199tmHF154gUaNGgGR9w4mSy+iiIikFg2JStKoSi3LQLVAgSoNP+6flUmOXyL124p32Dr3cZp0v5gZy7uUukbZc/2Ph7pmsPNERET8mQuyIi4VdOnSxS1dujTRYUiClZ0XBp6euboZaWzbVVDu/OysTBaP7hnRdUdOWYEDdiybxbZ5T5N50DE0GzSGA5o2LnWNYDEEWnAQyXkSnJktc851SXQcIiLxpCFRSVm+4c7rp6wIOC8sULIGnp62SIZHB3XO9iRrS6Z7krW2x9PszFuxWrXLDWGW3cst0OrQipwnIiLiT0OikpIC9VRVRKTDo83SdvHj4snUa38STQfciKV7fmQCDWEO6pwdUeIV6XkiIiI+StgkZfjPU0szoyjMcH5WZgZ7CouDJnW+1Zmhkqfbzv4TN2z9B0WN9sPS/lgQsSu/kBnLc5R4iYhIXChhk5RQtkctXLLmX11gwtzVASf6Q+DVmc45br31Vpo1a8YNN9wAVw7g7pkryc37Y4h1264C7Z8mIiJxozlskhIC7V8WjP+8MF/JqewI93hzznHTTTcxfvx4/ve//+GcY1DnbOrXKf9so/3TREQkXpSwSUqIZJ+yzIx0HhnaicWje5br9Ypkj7fi4mKuvfZaHn74Ya655hqefPJJzLt5bkV66ERERKIt4iFRM8sGWvm/xzn3QSyCEikr2P5l6WYUOxd2z7ZI9ngbMWIETz31FDfeeCMTJkwoSdZmLM/BgECDsNo/TURE4iGifdjM7AFgKPA14BuXcs65M8K8Lwt4Djgcz/93lzrnPjKza4C/AoXAbOfczd7zxwCXeT/jWufc3FDX1z5sqS/QhreBkq5Y7F9W9rM77viEve137rvvvpJkDaDr+PlBe9jAMwQbqGZoZTYAlvC0D5uI1ESR9rANAto55/ZU8PqPAnOcc0PMrDZQz8x6AAOBI5xze8ysOYCZHQYMAzoA+wPzzOwQ51zl9m2QpFeRQuhVqYIQ6rN37cmn4Nf15NCarRmdGTe4Y6lkDcIPe5aNWwXeRUQk2iJN2NYAGUDECZuZNQK6ARcDOOfygXwzuxoY70v+nHO/eN8yEHjNe/wHM/sOOBb4KNLPlOQSrpcpVCH0QIlNNPcvmzB3Nbt272HLrIfIW7OU/a94iryGTQN+drDh2GBxV/S+REREwol00cEuYIWZPW1mj/l+hXnPgcBm4AUzW25mz5lZfeAQ4CQz+8TMFpnZMd7zs4H1fu/f4D1WipkNN7OlZrZ08+bNEYYv8ebrZcrJzcPxRy+Tf3WBRBZCz9myg81vjWfX6v+SddL51GrYNOhnR1p83vdeFXgXEZFoizRhmwncC3wILPP7FUot4CjgSedcZ2AnMNp7vAlwPDAKeN08Y1AW4BrlJtg5555xznVxznVp1qxZhOFLvIXqZfIJNmG/KhP5feWq2oyeHbD81IzlOZxw7xx+eXMsed9+zF6nXEWjYwaF/OxBnbNpUi8j7Gf73huL+xIRkZotooTNOTcJmMwfidqr3mOhbAA2OOc+8X4/DU8CtwGY7jyWAMVAU+/xFn7vPwDYGOmNSHKJpJcp0FYbBvRoX7lEPFyvnu/11QumkbdmGXv1+SsNjxpQ8v6y23z4u+v0DuVi9ef/3ki2EBEREamIiBI2MzsZ+BZ4HHgC+J+ZdQv1Hufcz8B6M/P9L9ULzyrTGUBP73UPAWoDW/D04g0zszpm1gZoCyyp2O1Isoikl2lQ52zOOjq7VNeqA95YlhO2MHsg98xaGbJXz9fr17DLGexz7lgadupbcl6TehkhV5yWLdqelZlBk3oZAQu4q8C7iIhEW6SLDv4OnOqcWw0lidZk4Ogw77sGeMW7QnQNcAmeodHnzewrIB+4yHn2FllpZq/jSeoKgRFaIZq6RvVpF3AbjrK9TAtWbS437l2ZCfozluewbVdBwNc25uaxY8cOPn/5PrK6X0Sthk2p2/KIUufUq10r7OdVZNGDCryLiEg0RZqwZfiSNQDn3P/MLOykHufcCiDQfknnBzl/LDA2wpgkiUW6DUe0JuiHKhHVgD2ceuqp7PpmKfXan1SywKAqnyciIhJPkSZsS81sIvCS9/s/E37RgdRwgXqZym71kVUvI2DPWEUn6AdLuIrydvD91Dsp3LKOWx56htc27aOKBSIiknIiTdiuBkYA1+KZF/4BnrlsIhGbsTyHUVM/p6DYkzLl5OaRBmSkGwVFf6RRlZmgH2ivtKJd29n02m0UbM1h9swZnHbaaaTP+JJXPv6xVNKmBQEiIpLsIl0lusc597BzbrBz7kzn3D8qUfVAari7Z64sSdZ8ioFaaVblCfqj+rQrvy+MGZZRh8Muuo/TTjsNgPsGdeQfQztpQYCIiKSUkD1sZva6c+4cM/uSwHuiHRHgbSIB5eYFXhSQV1DM4tE9q3TtQZ2zWbpuK698/COFv28jrW4D0jMb0eaSfzD2rCPKnasETUREUkm4IdHrvL8PCHmWSBRVtnD6fYM60qp2HtdeMJz0fQ/hyAvuiHnRdRV5FxGReAiZsDnnfvJ++Rfn3C3+r5nZA8At5d8lEliTIAsM/KsIVKVw+po1axj7l3OoU7STuc+P57jjjiu5ZiySKhV5FxGReIm0NNUpAY71i2YgUv3ddXoHMtJLzzTLSDfuOr1DyfeRlLQK5Ntvv6V79+78mrudVheMZ9ibW+g6fj63z/gybE3TyqpsrCIiIhUVbg7b1cBfgIPM7Au/lxriqSsqErFI9marzL5sxcXFnHnmmez4fRdNzxnL9vqeCmc5uXnlVoRC5TbmrUhM0d7TTcOuIiISbg7bq8A7wDg8hdt9fnPObY1ZVFItBEs0QiUbgbbnAGicmUHX8fMDJi1paWm88MILDH/1S7bV2afU+wLtuQbRSaqCxeq/p1tVky0Nu4qICIQZEnXObXfOrQUeBbY659Y559YBBWZ2XDwClNQUrhB7MIEKp2ekGTvzC8td6+HJc3n44YcBOOaYY8gtk6yFEo2NcsMVea/sn4E/DbuKiAhEPoftSeB3v+93eo+JBBQu0ZixPIeu4+fTZvRsuo6fX5LEBCqc3qBurVIb6wLk/vgNN186hEceeYTt27cDwZOwsvuzVWaj3EDxhivyHo1kK17DriIiktwirXRg3gLtADjnis0s0vdKiqrKcF6oRCPcMF/ZYdM2o2eXusaenG/Y9PpdpGc25IMPPqBx48ZA8ILzZx2dzYJVm2M2LBnsWtFItiIZdhURkeov0qRrjZldyx+9an8B1sQmJEkGVZ07FSrRCNXzFOja/tfavf4rfpl2D+n1m3Dk8L/TunXrkvMiLThfURWNN1DcZY9HKlgSqlJaIiI1S6RDolcBfwJygA3AccDwWAUliVeV4bwZy3PYuaew3HFfolHRnif/uWKFuT9Tq2EzWl34ILcPPancuYM6Z7N4dE9+GN+fxaN7RmVifmV7ysLNcYtEuGFXERGpGSLqYXPO/QIMi3EskkQqm6SU7ZnzaVIvg7tO78CgztlMmLu6wj1PtQp2AXVo0LE32Uf34p4zO8ctaalsT1m0evxUSktERMLtw3azc+5BM/sngWuJXhuzyCShKpukBOqZA6hXu1ZJ0lGRYb4Zy3P46/hn2fjWQzQfcid1D+hAvovv9Mke7ZuV288t0p4yJVsiIhIN4f7n+8b7+9JYByLJpbJzpyLpmQvX8+S/2CHvfx+y6a0HqL3PgWQ0bQVEb+PbSMxYnsMby3JKJWsGnHW0EjEREYmfcLVEZ3l/nxSfcCRZVHY4L9KeuWA9T/5Dqju/+YAtsx6izn6H0Pyce0irU7/kvHhtaxGox9ABC1Ztjsvni4iIQPgh0VkE3ywe59wZUY9IkkZlhvOquqrRlyDtyVnlSdayD6X5kLtIq1Ov1Hn7Z2XGpWST9kETEZFkEG5I9CHv74OBfYGXvd+fC6yNUUySwqo60d6XCNXe/xCanHwJDTr1I6123VLnZGak06N9s7iUbKrMXD7V/hQRkWgzv/1wg59k9oFzrlu4Y/HWpUsXt3SpptdVJ23PGUPe3odQq3HzUsfTzSh2riQBCrbSNDsrk8WjewLRSZwCrXrNzEgPurVGRc+XijOzZc65LomOQ0QkniJdbtfMzA50zq0BMLM2QLPYhSWpqipJ0hW33Mt3U8fToFM/9u4zouR4oIRn5JQVAa/h66GLVtH0ivYYVnaTXRERkVAiTdhGAgvNzFfdoDVwZUwikoSqSsJVlSTp4pF3MOmR+8g85AT26v3Hnsz++7f5CzdUGUniFOm9VmQun+a8iYhILERU6cA5NwdoC1zn/dXOOTc3loFJ/PkSrpzcPBx/JFy+wuzhVLY6wtixY5n0yH3Ua38Szc64BUvPKHnNf/82f+GqCIRLnKp6r8EEm9um2p8iIlIVESVsZlYPGAX81Tn3OdDSzAbENDKJu6qUo4LK9S7t3r2bN998k/odetD09Juw9NKdvsHeG65kU1a9jIDv8x2vzL3OWJ5D1/HzaTN6Nl3Hzw+Y3EWjHJWIiEhZkdYSfQHIB07wfr8BuC8mEUnCVHU4ryK9S8458vPzqVu3LvPnz+eI88ZgaekRvddnUOdsRvVpx/5ZmWzMzWPC3NUlSVSwtTS+4xW910h75FT7U0REYiHShO0g59yDQAGAcy4Pz4bvUo1UdTgv0t4l5xw33HADp59+Ovn5+TRq1Iib+x1W4Z6pUEnU9ryCgO/xHa/ovVakRy4WBehFRKRmizRhyzezTLyb6JrZQcCemEUlCRFJwhVqWDCS3qXi4mL++te/8sgjj9C+fXsyMjIifm9ZoZKocAlZRYcutZhAREQSKdJVoncBc4AWZvYK0BW4OFZBSWJEUuNz1NTPKSj2jCvm5OYxaurnpd4bakVlcXExV155Jc899xyjRo3igQcewOyPjtqKVlYIlUT9Y2inkBUXKrpdR2U20BUREYmWsAmbmaUBTfBUOzgez1Dodc65LTGOTRIgVNJ098yVJcmaT0Gx4+6ZKyNKtG688Uaee+45br/9dv72t7+VStYqI1QSFUlCVpEEsaolt0RERKoibMLmnCs2s786514HZschJklSuUHmhQU7XtbFF1/M/vvvz6hRo6IST7gkqjK1UIOpasktERGRqoh0SPQ9M7sJmALs9B10zm2NSVRSbeTn5/PGG28wbNgwjjzySI488sioXTveSVQ0E0AREZGKiDRhuxTPgoO/lDl+YHTDkViIVjHyJvUy2LarfG9akyB7nu3Zs4ezzz6bWbNm0aZNG44//vgKf2Y4SqJERKQmiHSV6GHA48DnwArgn0CHGMUkURTNHf3vOr0DGeml551lpBt3nV7+n0JeXh6DBg1i1qxZPP744zFJ1kRERGqKSBO2ScChwGN4krVDvcckyVW1eoG/QZ2zmTDkyFJbb0wYcmS5Hq6dO3dy+umnM3fuXJ599ln+8peyHbMiIiJSEZEOibZzzvlPPlpgZp/HIiCJrmjvHxZsCNJ/2LXeLyv536IPePHFF7nwwgsr9TkiIiLyh0gTtuVmdrxz7mMAMzsOWBy7sCRa4rF/mG/YdVd+IWbGzuYdaHXVszTq2KvUOVphKSIiUjmRDokeB3xoZmvNbC3wEdDdzL40sy9iFp1UWTyKkU+Yu5rff9vOpldvIe/7TwEorN+0ZNg1knl0kRRWFxERqaki7WHrG9MoJGbisfXF+o0/88vrd5K/ZR2uuLjkuG/YNdQ8ukGds0sSOt85voTOP34REZGaLKKEzTm3LtaBSOxUZuuLSIcwf/nlF36dejv5WzbQfPAdZB54dMlrvmHXcPPowiV0IiIiNV2kQ6KVYmZZZjbNzFaZ2TdmdoLfazeZmTOzpn7HxpjZd2a22sz6xDK26ioaQ4uRbgWSm5vLySefTGHuT7QYek+pZM1/2DVcIXYVVhcREQktpgkb8CgwxznXHjgS+AbAzFoApwA/+k40s8OAYXj2d+sLPGFm6eWuKEFFa8+1SLcCady4Mf369ePdOXN47KYLS7b7yMrMoG5GGiOnrKDr+Pn0aN8s5Dy6cAmdiIhITRezhM3MGgHdgIkAzrl851yu9+V/ADfjqZ7gMxB4zTm3xzn3A/AdcGys4quOKrLnWqieuHA9XuvWreO7777DzPj73/9O9+7dGdQ5m8Wje/KPoZ3YU1jMtl0FJUnjG8tyOOvo7FL7t40b3LFkuDMeCyNERERSWaSLDirjQGAz8IKZHQksA64DegE5zrnPzUrtmp8NfOz3/QbvsVLMbDgwHKBly5axiTxFRTq0GG6Sf6itQNasWUOPHj3Iyspi+fLlpKWVzvmDJY0LVm1m8eieAeNTYXUREZHQYpmw1QKOAq5xzn1iZo8Cd+PpdTs1wPkW4Jgrd8C5Z4BnALp06VLu9Zos0j3Xwk3yH9WnXamEDjw9Xue1z6Bbt27k5eXx5ptvlkvWoPLz0VQTVEREJLhYzmHbAGxwzn3i/X4angSuDfC5dz+3A4DPzGxf7/kt/N5/ALAxhvFVO5EOLYZLqgZ1zmbc4I6lhjD/0qku9119Dvn5+SxYsICjjjoq4DU0H01ERCT6YtbD5pz72czWm1k759xqPEOhnznnSra/9yZtXZxzW8xsJvCqmT0M7A+0BZbEKr7qKNKhxUh64sr2eJ155pkALFy4kMMOOyxoDMF65zQfTUREpPJiOSQKcA3wipnVBtYAlwQ70Tm30sxeB74GCoERzrmiYOdLYJEMLVYmqXrxxRfZvHkzBx98cNjPB81HExERiSZzLnWngXXp0sUtXbo00WGkpEg2xl2yZAnjxo3jlVdeoV69egmKVKQ0M1vmnOuS6DhEROIp1j1skqTC9cR9+OGH9O3bl6ZNm7J161YlbCIiIgkU641zJcXMWJ7D4Vc+wokn96KgdiNufWIKBxxwQKLDEhERqdGUsEmJGctzuPbhl/n6+dHUatSMvYeN4++Lf61UeSsRERGJHg2JSokJc1dTVKcxdbLb0/T0UaTXzwpahD3S4vAiIiJSdUrYaqBAyVbb2rnkbNtFRtMW7DNsbKnzK1opIZLPU3InIiISOQ2J1jCBCsT/9f6nOLJTJ9JWvRfwPZFWSrjeW+zdfwg1WgXpRUREajIlbDVM2WRr59eLyHnjfuplt+PeG6+qUqUEKJ+QVaQgfVmhCtSLiIjUJErYahj/ZOv3r95ny9t/p84Bh9HkzLv4c7dDy5WkGje4Y8BKCaH4J2SVrS2qnjkREZE/aA5bDeMrS1W4YzO/zvkndVt2pNngOzigeROg8pUSyvIlZJEWpC8rXIF6ERGRmkQ9bDWMr0B8rUbN2Ofsv9HsrDupX79+hWp9+heHD8aXkEVakL6syvbMiYiIVEdK2GqYdR9M44zG68nOyiSz1RG0aJYVcNgznEGds1k8uiePDO0UMiHzT+5CDbOWFawHLlzPnIiISHWkIdFqrOx2Ggf/9D7/eux+hg0bxuLJV0TlMyIp9h7JMGtZlSlQLyIiUl2p+Hs1VXavtNzFk9n+31c4sc9AFrw9jVq1kj9X1/5tEoiKv4tITZT8/2tLpfgm7Tvn2P6fl9n+0RTqH96Tom4jUiJZg8r1zImIiFRHqfE/t1SY/+T84oLdNDjiVPbq+1d+2pGfwKhERESkMpSwVVP7Na7LjzkbqdVgL5r0vBxwmKVp0r6IiEgK0irRaqi4uJjGy//Fz5Ouo2jnNswMszRN2hcREUlRStiqmaKiIoYPH86cqf/ijCHn0mL//Sq0nYaIiIgkHw2JJkhlV0CGel9hYSGXXnopL730EnfeeSd33303ZhbrWxEREZEYU8KWAGW33PDVyQRCJm3h3jdhwgReeukl7rvvPm677bYY34WIiIjEi4ZEEyBUncyqvO+aa67hlVdeUbImIiJSzShhS4DK1skMVETdFebz1cyn+e2332jQoAHnnXdeVGIUERGR5KGELQEqWyczvcx8tOKCPfwy/T52fPQ68+bNi1p8IiIiklyUsCXAqD7tQhZMD6bIr4xYcf5uNr9xD7t/WM7e/a7lhk9q03X8fGYszwn7+TOW59B1/HzajJ4d8XtEREQkcbToIAEiKZgeSHZWJjm5eRTv2cUv0+5hT8437N1/JA0O74kjssULlV3wICIiIomjhC1BKlMnc1SfdoyZ/iU7cndQuOMXmg64kfqHdS91jm8RQrBrh1q4oIRNREQkOSlhSyG92zbGnXk4D71bm1qXP4Vl1Al4XqjFC5Vd8CAiIiKJozlsKWLLli2cdNJJLH71ERaP7sm6vw8muxKLFyq74EFEREQSRwlbCvjll1/o0aMHq1atonfv3iXHK7N4obILHkRERCRxNCSa5H766Sd69erF2rVrefvtt+nVq1fJa5VZvFDZBQ8iIiKSOOb8topINV26dHFLly5NdBgxU1hYSKdOnVi3bh2zZ8+mW7duiQ5JJOHMbJlzrkui4xARiSf1sCWxWrVqcf/999O0aVP+9Kc/JTocERERSZAakbDNWJ6TUkOA33//PStWrOCss87ijDPOSHQ4IiIikmDVPmFLtY1iV69eTc+ePSkqKqJPnz40aNAg0SGJiIhIglX7VaKhNopNNl9//TXdu3ensLCQefPmKVkTERERoAYkbKmyUewXX3zBySefTFpaGgsXLuTwww9PdEgiIiKSJKp9wpYqG8XOmjWLOnXqsGjRIg499NBEhyMiIiJJpNonbMm+UWxhYSEAt956KytWrKBt27YJjkhERESSTbVP2AZ1zmbc4I5kZ2ViQHZWJuMGd0yKBQeLFy+mffv2fP3115gZe++9d6JDEhERkSRU7VeJgidpS4YEzd/ChQsZMGAA2dnZNG7cONHhiIiISBKLaQ+bmWWZ2TQzW2Vm35jZCWY2wfv9F2b2ppll+Z0/xsy+M7PVZtYnlrEl0rx58zjttNNo1aoVixYtIjs7uZJJERERSS6xHhJ9FJjjnGsPHAl8A7wHHO6cOwL4HzAGwMwOA4YBHYC+wBNmlh7wqinso48+YsCAAbRt25aFCxey7777JjokERERSXIxS9jMrBHQDZgI4JzLd87lOufedc4Vek/7GDjA+/VA4DXn3B7n3A/Ad8CxsYovUTp16sSVV17J/PnzadasWaLDERERkRQQyx62A4HNwAtmttzMnjOz+mXOuRR4x/t1NrDe77UN3mOlmNlwM1tqZks3b94ci7hjYu7cueTm5pKZmcmjjz6qBQYiIiISsVgmbLWAo4AnnXOdgZ3AaN+LZnYbUAi84jsU4Bqu3AHnnnHOdXHOdUmVHqpXX32V/v37c8cddyQ6FBEREUlBsUzYNgAbnHOfeL+fhieBw8wuAgYAf3bOOb/zW/i9/wBgYwzji4tJkyZx/vnnc9JJJzFu3LhEhyMiIiIpKGYJm3PuZ2C9mfl2qO0FfG1mfYFbgDOcc7v83jITGGZmdcysDdAWWBKr+OLh2Wef5ZJLLqFXr17Mnj1btUFFRESkUmK9D9s1wCtmVhtYA1wCfArUAd4zM4CPnXNXOedWmtnrwNd4hkpHOOeKglw36e3atYv777+fvn37Mn36dOrWrZvokERERCRF2R8jkqmnS5cubunSpYkOoxznHGbG+vXrad68OXXq1El0SCLVhpktc851SXQcIiLxVO1LU8Xb+PHjufrqq3HO0aJFCyVrIiIiUmVK2KLEOcff/r+9+4+1uq7jOP58XViyfloY5grUmsH6J4py1HaLC/2gZJLN1BYbKlq5ZpPtVmIsyqaVNFtrK2tREzWL6ySaWyj9ETYmSkoQFC0tFSzo0krvYJXApz/OFz3Rvfzy3nu+557nY/vufs/nfL6fvd/3wM5r3+8593v99SxZsoSBgQEOHmzbq7mSJKlmDGzDoJTC0qVLWbZsGZdeeikrV65k/PiOuE2rJEkaBQa2YbB06VJuvPFGrrzySlasWMG4cWPujlqSJKmFDGzDoLu7m8WLF3PLLbfQ1eWvVJIkDS/TxUk6dOgQGzZsAGDu3LncfPPNhjVJkjQiTBgn4eDBg1xxxRV0d3ezefPmVpcjSZLGOD8Zf4IOHDjAZZddxu23386yZcuYPn16q0uSJEljnIHtBDz77LMsWLCAVatWccMNN3Dddde1uiRJktQBDGwnYM2aNaxatYrly5fT29vb6nIkSVKHMLCdgAsvvJAHHniAmTNntroUSZLUQfzSwTHs37+fiy+++LkvFxjWJEnSaDOwHcW+ffuYN28efX19bN++vdXlSJKkDuUl0SEMDAxw3nnnsWHDBlauXMmCBQtaXZIkSepQBrZBPP3008ydO5dNmzZx5513ctFFF7W6JEmS1MEMbIOYMGECp59+On19fVxwwQWtLkeSJHU4A1uT/v5+urq6mDhxIqtXryZJq0uSJEnySweH7d69m56eHubPn08pxbAmSZJqwzNswFNPPcWcOXPYuXMn99xzj2FNkiTVSscHtieffJLZs2ezZ88e1q5dS3d3d6tLkiRJ+h8dH9gWLVrE3r17WbdunX8UV5Ik1VLHB7YVK1bQ39/PjBkzWl2KJEnSoDrySwc7duygt7eXQ4cOMWXKFMOaJEmqtY4LbNu2bWPWrFncdttt7Nq1q9XlSJIkHVNHBbYtW7bQ09NDV1cX69evZ8qUKa0uSZIk6Zg6JrA9/PDD9PT0MGHCBNavX8+0adNaXZIkSdJx6ZjANjAwwBlnnMH999/POeec0+pyJEmSjlvHfEt01qxZbN26lXHjxrW6FEmSpBPSMWfYAMOaJElqSx0V2CRJktqRgU2SJKnmDGySJEk1Z2CTJEmqOQObJElSzRnYJEmSas7AJkmSVHMGNkmSpJozsEmSJNVcSimtruGkJekHnmh1HSfoNGBvq4sYBmOlDxg7vXRKH2eWUl49WsVIUh20dWBrR0l+XUp5W6vreKHGSh8wdnqxD0kau7wkKkmSVHMGNkmSpJozsI2+77W6gGEyVvqAsdOLfUjSGOVn2CRJkmrOM2ySJEk1Z2CTJEmqOQPbMEtyapK7kuxI8vsk70iyvHq8NcnqJKc2zV+S5NEkf0jy/haW/n8G66Xpud4kJclpTWO17GWoPpJcXdW6PclNTfPbpo8k05NsTPKbJL9Ocm7T/Lr2MbWq9/D2TJJrkrwqybokf6x+vrLpmFr2IkmjppTiNowbcCtwRbX/IuBU4H3A+Grsa8DXqv03AVuAU4CzgceAca3u4Wi9VPuTgXtp/NHi0+reyxCvSQ/wC+CUanxSm/ZxH/CBauyDwC/r3scRPY0DdgNnAjcB11bj17bL/xM3Nze30dg8wzaMkrwceBewAqCU8p9Syj9LKfeVUg5U0zYCr6v25wM/LqX8u5TyZ+BR4Nwj122FoXqpnv4G8Fmg+RsrtezlKH1cBXy1lPLvavxv1SHt1kcBXl5NewXwl2q/ln0MYg7wWCnlCRo131qN3wp8qNpvl14kacQY2IbX64F+4IdJNif5fpKXHDHncuDn1f5rgZ1Nz+2qxupg0F6SnA88VUrZcsT8uvYy1GvyRqA7yYNJ1id5ezW/3fq4BlieZCfwdWBJNb+ufRzpEuDOav/0UspfAaqfk6rxdulFkkaMgW14jQfeCnynlPIWYB+NSzsAJPk8cAC44/DQIGvU5e+sDNbLF4HPA18YZH5dexnqNRkPvBKYCXwGWJUktF8fVwGLSymTgcVUZ+Cobx/PSfIi4Hyg71hTBxmrVS+SNNIMbMNrF7CrlPJg9fguGm+yJFkIzAM+VkopTfMnNx3/Op6/pNVqQ/VyNrAlyeM06n0kyWuoby9D9bELuLs0PAQconHT8XbrYyFwdzXWx/OXCuvaR7MPAI+UUvZUj/ckOQOg+nn4MnU79CJJI8rANoxKKbuBnUmmVkNzgN8lmQt8Dji/lLK/6ZCfAZckOSXJ2cA5wEOjWvQQhujlkVLKpFLKWaWUs2i8kb61mlvLXoZ6TYCfArMBkryRxof499J+ffwFeHc1Nhv4Y7Vfyz6O8FGevxwKjZoXVvsLgTVN43XvRZJG1PhWFzAGXQ3cUV3u+RNwGbCJxjfc1jWuurGxlPLJUsr2JKtovPEeAD5VSjnYoroHM1gvg6p5L4P1sQ/4QZJtwH+AhdWZz3brYw3wzSTjgX8BH4favx4keTHwXuATTcNfpXFpehHwJPARqH8vkjQavDWVJElSzXlJVJIkqeYMbJIkSTVnYJMkSao5A5skSVLNGdgkSZJqzsCmWkryeJLT6rKOJEmtZGDTiEqD/84kSXoBfCPVsEtyVpLfJ/k28AgwOclnkmxKsjXJl5rm/jTJw0m2J/n4Mda9KslNTY8vTfKt41mnqmlb0+PeJF+s9t+QZG11/K+STHvBvwRJkoaRgU0jZSqwsrpR+VQatxM6F5gOzEjyrmre5aWUGcDbgE8nmXiUNe8CPtz0+GLgJyexzpG+B1xdHd8LfPsEjpUkacR5ayqNlCdKKRur/fdV2+bq8UtpBLj7aYSrC6rxydX43wdbsJTSn+RPSWbSuGfmVGBD9fRxr9MsyUuBdwJ91W3DoHEbMUmSasPAppGyr2k/wFdKKd9tnpBkFvAe4B2llP1JfglMOMa6PwEuAnYAq0sp5TjXOcD/nlE+/HwX8M9SyvTjaUqSpFbwkqhGw73A5dXZLJK8Nskk4BXAP6qQNQ2YeRxr3Q18CPgoz18OPZ519gCTkkxMcgowD6CU8gzw5yQfqWpLkjefbKOSJI0EA5tGXCnlPuBHwANJfkvjs2gvA9YC45NsBb4MbBx6lefW+gfwO+DMUspD1fAx1ymlPAtcDzwI3EPjDN1hHwMWJdkCbAfmn0yfkiSNlJRSWl2DJEmSjsIzbJIkSTVnYJMkSao5A5skSVLNGdgkSZJqzsAmSZJUcwY2SZKkmjOwSZIk1dx/AZ/I82BlVe9UAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import pylab as pl\n",
"pl.figure(figsize=(10,10))\n",
"for i, col in [(1, \"mathscr\"), (2, \"readscr\"), (3, \"testscr\")]:\n",
" ax = pl.subplot(2,2,i)\n",
" ax.scatter(data_test[col], test_data_predictions_2[col])\n",
" minval, maxval = data_test[col].min(), data_test[col].max()\n",
" ax.plot([minval, maxval], [minval, maxval], ls=\"--\", c=\"k\")\n",
" ax.set_xlabel(\"real value\")\n",
" ax.set_ylabel(\"prediction\")\n",
" ax.set_title(col)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and by the R2 score."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"performance_2 = causal_structure_2.evaluate_objective(Evaluator, data=data_test,\n",
" inputs=inputs, outputs=outputs,\n",
" metric=\"r2\")"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'performance of simple model:'"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"{'readscr': 0.8020755684153762,\n",
" 'mathscr': 0.6574450121854509,\n",
" 'testscr': 0.7577271442716229}"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"'performance of structured model:'"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"{'readscr': 0.7736833098673347,\n",
" 'testscr': 0.733377201980135,\n",
" 'mathscr': 0.6371819769563171}"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(\"performance of simple model:\")\n",
"display({\n",
" key: value for key, value in performance_1.items() if key in outputs\n",
"})\n",
"display(\"performance of structured model:\")\n",
"display({\n",
" key: value for key, value in performance_2.items() if key in outputs\n",
"})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that the R2 scores are actually slightly worse than for the simple causal structure, which only had inputs and outputs. Maybe some of our assumptions were incomplete or even wrong.\n",
"\n",
"So what have we gained if the prediction performance is not better? We can see that once we start comparing predicions on partial data."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Comparing the performance of the simple and the advanced structure."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we established the performance when predicting the outputs based on all inputs is the same in both models. When all input data are available and the training set and the test set are split randomly it does not matter that much whether the learned behavior is based on real causality or confounders.\n",
"\n",
"However, if we make predictions from incomplete data, this changes."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Example 1: predicting based on subsidized meals and CalWorks\n",
"Let's try a prediction based on only the amount of subsidized meals and the percentage in the CalWorks program."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'prediction based on simple model:'"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"{'readscr': 0.6342623417400537,\n",
" 'mathscr': 0.5509301033970446,\n",
" 'testscr': 0.590015453901932}"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"'prediction based on structured model:'"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"{'readscr': 0.7078249997685941,\n",
" 'testscr': 0.6881614958989846,\n",
" 'mathscr': 0.6110603334844447}"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(\"prediction based on simple model:\")\n",
"display({\n",
" key: value for key, value in \n",
" causal_structure_1.evaluate_objective(Evaluator, data=data_test,\n",
" inputs=['calwpct', 'mealpct'], outputs=outputs,\n",
" metric=\"r2\").items()\n",
" if key in outputs\n",
"})\n",
"display(\"prediction based on structured model:\")\n",
"display({\n",
" key: value for key, value in \n",
" causal_structure_2.evaluate_objective(Evaluator, data=data_test,\n",
" inputs=['calwpct', 'mealpct'], outputs=outputs,\n",
" metric=\"r2\").items()\n",
" if key in outputs\n",
"})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that in this scenario the structured model outperforms the simple model, but both models still achieve decent scores."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Example 2: predicting based on CalWorks only\n",
"Let's try a prediction based on only the percentage in the CalWorks program."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'prediction based on simple model:'"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"{'readscr': -0.033142940725477965,\n",
" 'mathscr': 0.11518897669104289,\n",
" 'testscr': 0.029182778885538996}"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"'prediction based on structured model:'"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"{'readscr': 0.32351923326782916,\n",
" 'testscr': 0.3688571477752969,\n",
" 'mathscr': 0.37847951309032846}"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(\"prediction based on simple model:\")\n",
"display({\n",
" key: value for key, value in \n",
" causal_structure_1.evaluate_objective(Evaluator, data=data_test,\n",
" inputs=['calwpct'], outputs=outputs,\n",
" metric=\"r2\").items()\n",
" if key in outputs\n",
"})\n",
"display(\"prediction based on structured model:\")\n",
"display({\n",
" key: value for key, value in \n",
" causal_structure_2.evaluate_objective(Evaluator, data=data_test,\n",
" inputs=['calwpct'], outputs=outputs,\n",
" metric=\"r2\").items()\n",
" if key in outputs\n",
"})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this scenario the structured model outperforms the simple model significantly."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Example 3: predicting based on district average income\n",
"Let's try a prediction based on only the district average income."
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'prediction based on simple model:'"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"{'readscr': 0.2816587802736583,\n",
" 'mathscr': 0.2002535289464773,\n",
" 'testscr': 0.2320481281685418}"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"'prediction based on structured model:'"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"{'readscr': 0.34039792627046295,\n",
" 'testscr': 0.3543316153774535,\n",
" 'mathscr': 0.3358698650229551}"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(\"prediction based on simple model:\")\n",
"display({\n",
" key: value for key, value in \n",
" causal_structure_1.evaluate_objective(Evaluator, data=data_test,\n",
" inputs=['avginc'], outputs=outputs,\n",
" metric=\"r2\").items()\n",
" if key in outputs\n",
"})\n",
"display(\"prediction based on structured model:\")\n",
"display({\n",
" key: value for key, value in \n",
" causal_structure_2.evaluate_objective(Evaluator, data=data_test,\n",
" inputs=['avginc'], outputs=outputs,\n",
" metric=\"r2\").items()\n",
" if key in outputs\n",
"})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this scenario the structured model outperforms the simple model significantly."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To summarize, a causal structure based on realistic dependency assumptions will not necessarily outperform a black-box model with only inputs and outputs. It is however, more robust when making predicions based on incomplete data."
]
}
],
"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.9"
}
},
"nbformat": 4,
"nbformat_minor": 4
}